Merge branch 'master' into reconnection
This commit is contained in:
commit
5470bbd0da
@ -876,9 +876,12 @@ 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"
|
||||
case default
|
||||
update_errors => update_errors_l2
|
||||
name_norm = "L2-norm"
|
||||
@ -4178,16 +4181,139 @@ 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:
|
||||
! ---------------------------
|
||||
!
|
||||
! 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;
|
||||
!
|
||||
!===============================================================================
|
||||
!
|
||||
@ -4195,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(:))
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
@ -4258,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
|
||||
#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)
|
||||
@ -4323,7 +4531,7 @@ module evolution
|
||||
|
||||
!-------------------------------------------------------------------------------
|
||||
!
|
||||
end subroutine update_errors_max
|
||||
end subroutine update_errors_lmax
|
||||
|
||||
!===============================================================================
|
||||
!
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user