From de49eeedaa75a3d63f4c47af63fe2504f8317a9e Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 28 Dec 2023 10:49:24 -0300 Subject: [PATCH 1/4] EVOLUTION: Add subroutine to calculate errors using the L1-norm. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 123 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 123 insertions(+) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 852212a..91e7103 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -879,6 +879,9 @@ module evolution case("max", "maximum", "infinity") update_errors => update_errors_max name_norm = "maximum norm" + case("L1", "l1", "Manhattan", "manhattan") + update_errors => update_errors_l1 + name_norm = "L1-norm" case default update_errors => update_errors_l2 name_norm = "L2-norm" @@ -4178,6 +4181,126 @@ module evolution ! !=============================================================================== ! +! subroutine UPDATE_ERRORS_L1: +! --------------------------- +! +! Description: +! +! The subroutine is designed to update errors with respect to absolute and +! relative tolerances using the L1 norm. The errors are calculated for +! each variable. +! +! Arguments: +! +! nh - the index of the higher order solution; +! nl - the index of the lower order solution; +! +!=============================================================================== +! + subroutine update_errors_l1(nh, nl) + + use blocks , only : block_data, data_blocks + use blocks , only : get_nleafs, get_dblocks + use coordinates, only : nc => ncells, nb, ne + use equations , only : nf, errors + use helpers , only : print_message +#ifdef MPI + use mpitools , only : reduce_sum +#endif /* MPI */ + use workspace , only : resize_workspace, work, work_in_use + + implicit none + + integer, intent(in) :: nh, nl + + character(len=*), parameter :: loc = 'EVOLUTION:update_errors_l1()' + + type(block_data), pointer :: pdata + + integer :: nblks, n, p + + logical , save :: first = .true. + integer , save :: nt = 0 + real(kind=8), save :: fnorm = 1.0d+00 + + real(kind=8), dimension(nf) :: lerrors + real(kind=8), dimension(:,:,:), pointer, save :: uhi, ulo, div + +!$omp threadprivate(first,nt,uhi,ulo,div) +!$ integer :: omp_get_thread_num + +!------------------------------------------------------------------------------- +! + errors(:) = 0.0d+00 + + nblks = get_dblocks() + +!$omp parallel default(shared) private(pdata,lerrors,n,p) + if (first) then +!$ nt = omp_get_thread_num() + n = nc**NDIMS + +#if NDIMS == 3 + uhi(1:nc,1:nc,1:nc) => work( 1: n,nt) + ulo(1:nc,1:nc,1:nc) => work( n+1:2*n,nt) + div(1:nc,1:nc,1:nc) => work(2*n+1:3*n,nt) +#else /* NDIMS == 3 */ + uhi(1:nc,1:nc,1: 1) => work( 1: n,nt) + ulo(1:nc,1:nc,1: 1) => work( n+1:2*n,nt) + div(1:nc,1:nc,1: 1) => work(2*n+1:3*n,nt) +#endif /* NDIMS == 3 */ + +!$omp single + fnorm = 1.0d+00 / n +!$omp end single + + first = .false. + end if + + if (work_in_use(nt)) & + call print_message(loc, "Workspace is being used right now! " // & + "Corruptions can occur!") + + work_in_use(nt) = .true. +!$omp barrier + +!$omp do reduction(+:errors) + do n = 1, nblks + pdata => data_blocks(n)%ptr + + do p = 1, nf +#if NDIMS == 3 + uhi = pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) + ulo = pdata%uu(p,nb:ne,nb:ne,nb:ne,nl) +#else /* NDIMS == 3 */ + uhi = pdata%uu(p,nb:ne,nb:ne, : ,nh) + ulo = pdata%uu(p,nb:ne,nb:ne, : ,nl) +#endif /* NDIMS == 3 */ + + div = atol + rtol * max(abs(uhi), abs(ulo)) + + lerrors(p) = fnorm * sum(abs(uhi - ulo) / div) + end do + + errors = errors + lerrors + end do +!$omp end do nowait + + work_in_use(nt) = .false. +!$omp end parallel + + errors(:) = errors(:) / get_nleafs() + +#ifdef MPI + call reduce_sum(errors) +#endif /* MPI */ + +!------------------------------------------------------------------------------- +! + end subroutine update_errors_l1 +! +!=============================================================================== +! ! subroutine UPDATE_ERRORS_L2: ! --------------------------- ! From 213879c9c44f6cb4a0e7e14313c5f5b989a5af26 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 28 Dec 2023 10:56:11 -0300 Subject: [PATCH 2/4] EVOLUTION: Update subroutine to calculate errors using L2-norm. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 94 +++++++++++++++++++++++++++++++------------ 1 file changed, 68 insertions(+), 26 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 91e7103..7f0f10b 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -4304,13 +4304,16 @@ module evolution ! subroutine UPDATE_ERRORS_L2: ! --------------------------- ! -! Subroutine updates the errors with respect to the absolute and relative -! tolerances using the L2 norm. The errors are calculated per variable. +! Description: ! -! Arguments: +! The subroutine is designed to update errors with respect to absolute and +! relative tolerances using the L2 norm. The errors are calculated for +! each variable. ! -! nh - the index of the higher order solution -! nl - the index of the lower order solution +! Arguments: +! +! nh - the index of the higher order solution; +! nl - the index of the lower order solution; ! !=============================================================================== ! @@ -4318,62 +4321,101 @@ module evolution use blocks , only : block_data, data_blocks use blocks , only : get_nleafs, get_dblocks - use coordinates, only : ncells, nb, ne + use coordinates, only : nc => ncells, nb, ne use equations , only : nf, errors + use helpers , only : print_message #ifdef MPI use mpitools , only : reduce_sum #endif /* MPI */ + use workspace , only : resize_workspace, work, work_in_use implicit none integer, intent(in) :: nh, nl + character(len=*), parameter :: loc = 'EVOLUTION:update_errors_l2()' + type(block_data), pointer :: pdata - integer :: l, n, p - real(kind=8) :: err + integer :: nblks, n, p logical , save :: first = .true. + integer , save :: nt = 0 real(kind=8), save :: fnorm = 1.0d+00 + real(kind=8), dimension(nf) :: lerrors + real(kind=8), dimension(:,:,:), pointer, save :: uhi, ulo, div + +!$omp threadprivate(first,nt,uhi,ulo,div) +!$ integer :: omp_get_thread_num + !------------------------------------------------------------------------------- ! + errors(:) = 0.0d+00 + + nblks = get_dblocks() + +!$omp parallel default(shared) private(pdata,lerrors,n,p) if (first) then - fnorm = 1.0d+00 / ncells**NDIMS +!$ nt = omp_get_thread_num() + n = nc**NDIMS + +#if NDIMS == 3 + uhi(1:nc,1:nc,1:nc) => work( 1: n,nt) + ulo(1:nc,1:nc,1:nc) => work( n+1:2*n,nt) + div(1:nc,1:nc,1:nc) => work(2*n+1:3*n,nt) +#else /* NDIMS == 3 */ + uhi(1:nc,1:nc,1: 1) => work( 1: n,nt) + ulo(1:nc,1:nc,1: 1) => work( n+1:2*n,nt) + div(1:nc,1:nc,1: 1) => work(2*n+1:3*n,nt) +#endif /* NDIMS == 3 */ + +!$omp single + fnorm = 1.0d+00 / n +!$omp end single first = .false. end if - errors(:) = 0.0d+00 + if (work_in_use(nt)) & + call print_message(loc, "Workspace is being used right now! " // & + "Corruptions can occur!") - n = get_dblocks() + work_in_use(nt) = .true. +!$omp barrier + +!$omp do reduction(+:errors) + do n = 1, nblks + pdata => data_blocks(n)%ptr -!$omp parallel do default(shared) private(pdata,p,err) - do l = 1, n - pdata => data_blocks(l)%ptr do p = 1, nf #if NDIMS == 3 - err = sum((abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) & - - pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)) / (atol + rtol * & - max(abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh)), & - abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)))))**2) + uhi = pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) + ulo = pdata%uu(p,nb:ne,nb:ne,nb:ne,nl) #else /* NDIMS == 3 */ - err = sum((abs(pdata%uu(p,nb:ne,nb:ne, : ,nh) & - - pdata%uu(p,nb:ne,nb:ne, : ,nl)) / (atol + rtol * & - max(abs(pdata%uu(p,nb:ne,nb:ne, : ,nh)), & - abs(pdata%uu(p,nb:ne,nb:ne, : ,nl)))))**2) + uhi = pdata%uu(p,nb:ne,nb:ne, : ,nh) + ulo = pdata%uu(p,nb:ne,nb:ne, : ,nl) #endif /* NDIMS == 3 */ -!$omp atomic update - errors(p) = errors(p) + fnorm * err + + div = atol + rtol * max(abs(uhi), abs(ulo)) + + lerrors(p) = fnorm * sum((abs(uhi - ulo) / div)**2) end do + + errors = errors + lerrors end do -!$omp end parallel do +!$omp end do nowait + + work_in_use(nt) = .false. +!$omp end parallel + + errors(:) = errors(:) / get_nleafs() #ifdef MPI call reduce_sum(errors) #endif /* MPI */ - errors(:) = sqrt(errors(:) / get_nleafs()) + errors(:) = sqrt(errors(:)) !------------------------------------------------------------------------------- ! From 60d13dff6cd4f0335906b9bdde65a91e7de59627 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 28 Dec 2023 11:14:34 -0300 Subject: [PATCH 3/4] EVOLUTION: Update subroutine to calculate errors using Lmax-norm. Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 105 +++++++++++++++++++++++++++++------------- 1 file changed, 74 insertions(+), 31 deletions(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 7f0f10b..8db757e 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -876,9 +876,9 @@ module evolution call get_parameter("error_norm", error_norm) select case(trim(error_norm)) - case("max", "maximum", "infinity") - update_errors => update_errors_max - name_norm = "maximum norm" + case("Lmax", "lmax", "max", "maximum", "infinity") + update_errors => update_errors_lmax + name_norm = "Lmax-norm" case("L1", "l1", "Manhattan", "manhattan") update_errors => update_errors_l1 name_norm = "L1-norm" @@ -4423,64 +4423,107 @@ module evolution ! !=============================================================================== ! -! subroutine UPDATE_ERRORS_MAX: -! ---------------------------- +! subroutine UPDATE_ERRORS_LMAX: +! ----------------------------- ! -! Subroutine updates the errors with respect to the absolute and relative -! tolerances using the maximum norm. The errors are calculated per variable. +! Description: ! -! Arguments: +! The subroutine is designed to update errors with respect to absolute and +! relative tolerances using the Lmax norm. The errors are calculated for +! each variable. ! -! nh - the index of the higher order solution -! nl - the index of the lower order solution +! Arguments: +! +! nh - the index of the higher order solution; +! nl - the index of the lower order solution; ! !=============================================================================== ! - subroutine update_errors_max(nh, nl) + subroutine update_errors_lmax(nh, nl) use blocks , only : block_data, data_blocks, get_dblocks - use coordinates, only : nb, ne + use coordinates, only : nc => ncells, nb, ne use equations , only : nf, errors + use helpers , only : print_message #ifdef MPI - use mpitools , only : reduce_maximum + use mpitools , only : reduce_sum #endif /* MPI */ + use workspace , only : resize_workspace, work, work_in_use implicit none integer, intent(in) :: nh, nl - integer :: l, n, p - real(kind=8) :: err + character(len=*), parameter :: loc = 'EVOLUTION:update_errors_lmax()' type(block_data), pointer :: pdata + integer :: nblks, n, p + + logical , save :: first = .true. + integer , save :: nt = 0 + + real(kind=8), dimension(nf) :: lerrors + real(kind=8), dimension(:,:,:), pointer, save :: uhi, ulo, div + +!$omp threadprivate(first,nt,uhi,ulo,div) +!$ integer :: omp_get_thread_num + !------------------------------------------------------------------------------- ! errors(:) = 0.0d+00 - n = get_dblocks() + nblks = get_dblocks() + +!$omp parallel default(shared) private(pdata,lerrors,n,p) + if (first) then +!$ nt = omp_get_thread_num() + n = nc**NDIMS + +#if NDIMS == 3 + uhi(1:nc,1:nc,1:nc) => work( 1: n,nt) + ulo(1:nc,1:nc,1:nc) => work( n+1:2*n,nt) + div(1:nc,1:nc,1:nc) => work(2*n+1:3*n,nt) +#else /* NDIMS == 3 */ + uhi(1:nc,1:nc,1: 1) => work( 1: n,nt) + ulo(1:nc,1:nc,1: 1) => work( n+1:2*n,nt) + div(1:nc,1:nc,1: 1) => work(2*n+1:3*n,nt) +#endif /* NDIMS == 3 */ + + first = .false. + end if + + if (work_in_use(nt)) & + call print_message(loc, "Workspace is being used right now! " // & + "Corruptions can occur!") + + work_in_use(nt) = .true. +!$omp barrier + +!$omp do reduction(max:errors) + do n = 1, nblks + pdata => data_blocks(n)%ptr -!$omp parallel do default(shared) private(pdata,p,err) - do l = 1, n - pdata => data_blocks(l)%ptr do p = 1, nf #if NDIMS == 3 - err = maxval(abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) & - - pdata%uu(p,nb:ne,nb:ne,nb:ne,nl)) / (atol + rtol * & - max(abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nh)), & - abs(pdata%uu(p,nb:ne,nb:ne,nb:ne,nl))))) + uhi = pdata%uu(p,nb:ne,nb:ne,nb:ne,nh) + ulo = pdata%uu(p,nb:ne,nb:ne,nb:ne,nl) #else /* NDIMS == 3 */ - err = maxval(abs(pdata%uu(p,nb:ne,nb:ne, : ,nh) & - - pdata%uu(p,nb:ne,nb:ne, : ,nl)) / (atol + rtol * & - max(abs(pdata%uu(p,nb:ne,nb:ne, : ,nh)), & - abs(pdata%uu(p,nb:ne,nb:ne, : ,nl))))) + uhi = pdata%uu(p,nb:ne,nb:ne, : ,nh) + ulo = pdata%uu(p,nb:ne,nb:ne, : ,nl) #endif /* NDIMS == 3 */ -!$omp atomic update - errors(p) = max(errors(p), err) + + div = atol + rtol * max(abs(uhi), abs(ulo)) + + lerrors(p) = maxval(abs(uhi - ulo) / div) end do + errors = max(errors, lerrors) end do -!$omp end parallel do +!$omp end do nowait + + work_in_use(nt) = .false. +!$omp end parallel #ifdef MPI call reduce_maximum(errors) @@ -4488,7 +4531,7 @@ module evolution !------------------------------------------------------------------------------- ! - end subroutine update_errors_max + end subroutine update_errors_lmax !=============================================================================== ! From 659449669bfb1c3e202d7b87947f54c155a5191e Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Thu, 28 Dec 2023 11:52:22 -0300 Subject: [PATCH 4/4] EVOLUTION: Fix wrong import in update_errors_lmax(). Signed-off-by: Grzegorz Kowal --- sources/evolution.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sources/evolution.F90 b/sources/evolution.F90 index 8db757e..423906f 100644 --- a/sources/evolution.F90 +++ b/sources/evolution.F90 @@ -4446,7 +4446,7 @@ module evolution use equations , only : nf, errors use helpers , only : print_message #ifdef MPI - use mpitools , only : reduce_sum + use mpitools , only : reduce_maximum #endif /* MPI */ use workspace , only : resize_workspace, work, work_in_use