diff options
Diffstat (limited to 'src/Dissipation_4_3_min_err_coeff_alt.F90')
| -rw-r--r-- | src/Dissipation_4_3_min_err_coeff_alt.F90 | 248 |
1 files changed, 147 insertions, 101 deletions
diff --git a/src/Dissipation_4_3_min_err_coeff_alt.F90 b/src/Dissipation_4_3_min_err_coeff_alt.F90 index ef619d0..3310704 100644 --- a/src/Dissipation_4_3_min_err_coeff_alt.F90 +++ b/src/Dissipation_4_3_min_err_coeff_alt.F90 @@ -6,7 +6,9 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, & - delta, epsilon, dfl, rhs) + delta, epsilon, dfl, npatches, patch, rhs) + + use dissipation_coeff implicit none @@ -20,65 +22,88 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, & CCTK_INT, dimension(6), intent(in) :: bb CCTK_INT, dimension(3), intent(in) :: gsize CCTK_REAL, dimension(3), intent(in) :: delta - CCTK_REAL, intent(in) :: epsilon, dfl + CCTK_REAL, intent(in) :: epsilon + CCTK_REAL, dimension(3), intent(in) :: dfl + CCTK_INT, intent(in) :: npatches, patch CCTK_REAL :: zero = 0.0 integer, parameter :: wp = kind(zero) integer :: i, j, k, center - CCTK_REAL, dimension(:,:), allocatable :: a, d, b, h, tmp + CCTK_REAL, dimension(:,:), allocatable :: atmp, d, b, h, tmp + CCTK_REAL, dimension(:,:), pointer :: a CCTK_REAL :: idel CCTK_INT :: il, ir, jl, jr, kl, kr ni = lsh(1); nj = lsh(2); nk = lsh(3) - allocate ( a(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni), tmp(ni,ni) ) - a = zero; d = zero; b = zero; h = zero; tmp = zero + if ( first ) then + allocate ( savedx(npatches), savedy(npatches), savedz(npatches) ) + savedx = .false.; savedy = .false.; savedz = .false. + allocate ( xcoeff(npatches), ycoeff(npatches), zcoeff(npatches) ) + first = .false. + end if - call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl ) + if ( .not. savedx(patch) ) then + allocate ( atmp(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni), tmp(ni,ni) ) + atmp = zero; d = zero; b = zero; h = zero; tmp = zero - call set_hmatrix ( h, bb(1:2) ) + call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl(1) ) - if ( any ( bb(1:2) == 0 ) ) then + call set_hmatrix ( h, bb(1:2) ) - call set_dmatrix ( d, bb(1:2) ) + if ( any ( bb(1:2) == 0 ) ) then - a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + call set_dmatrix ( d, bb(1:2) ) - else + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - if ( lsh(1) /= gsh(1) ) then - call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) - end if - - center = gsh(1) / 2 - ir = center + max(3,gsize(1)) + else - call set_dmatrix_half ( d(1:ir,1:ir), 1 ) + if ( lsh(1) /= gsh(1) ) then + call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) + end if - tmp(1:ir,1:ir) = -transpose ( & - matmul ( h(1:ir,1:ir), & - matmul ( & - transpose ( d(1:ir,1:ir) ), & - matmul ( b(1:ir,1:ir), & - d(1:ir,1:ir) ) ) ) ) + center = gsh(1) / 2 + ir = center + max(3,gsize(1)) - a(1:ir,1:center) = tmp(1:ir,1:center) + call set_dmatrix_half ( d(1:ir,1:ir), 1 ) - il = center + 1 - max(3,gsize(1)) - d = zero - call set_dmatrix_half ( d(il:ni,il:ni), 2 ) - tmp(il:ni,il:ni) = -transpose ( & - matmul ( h(il:ni,il:ni), & + tmp(1:ir,1:ir) = -transpose ( & + matmul ( h(1:ir,1:ir), & matmul ( & - transpose ( d(il:ni,il:ni) ), & - matmul ( b(il:ni,il:ni), & - d(il:ni,il:ni) ) ) ) ) + transpose ( d(1:ir,1:ir) ), & + matmul ( b(1:ir,1:ir), & + d(1:ir,1:ir) ) ) ) ) + + atmp(1:ir,1:center) = tmp(1:ir,1:center) + + il = center + 1 - max(3,gsize(1)) + d = zero + call set_dmatrix_half ( d(il:ni,il:ni), 2 ) + tmp(il:ni,il:ni) = -transpose ( & + matmul ( h(il:ni,il:ni), & + matmul ( & + transpose ( d(il:ni,il:ni) ), & + matmul ( b(il:ni,il:ni), & + d(il:ni,il:ni) ) ) ) ) + + atmp(il:ni,center+1:ni) = tmp(il:ni,center+1:ni) - a(il:ni,center+1:ni) = tmp(il:ni,center+1:ni) + end if + + allocate ( xcoeff(patch)%coeff(ni,ni) ) + + xcoeff(patch)%coeff = atmp + + savedx(patch) = .true. + + deallocate ( atmp, d, b, h, tmp ) end if + a => xcoeff(patch)%coeff + idel = epsilon / ( 64 * delta(1) ) if ( bb(1) == 0 ) then @@ -180,59 +205,69 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, & end do - deallocate ( a, d, b, h, tmp ) - if ( zero_derivs_y == 0 ) then - allocate ( a(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj), tmp(nj,nj) ) - a = zero; d = zero; b = zero; h = zero; + if ( .not. savedy(patch) ) then + allocate ( atmp(nj,nj), d(nj,nj), b(nj,nj), h(nj,nj), tmp(nj,nj) ) + atmp = zero; d = zero; b = zero; h = zero; tmp = zero - call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl ) + call set_bmatrix ( b, bb(3:4), lsh(2), gsh(2), lbnd(2), delta(2), dfl(2) ) - call set_hmatrix ( h, bb(3:4) ) + call set_hmatrix ( h, bb(3:4) ) - if ( any ( bb(3:4) == 0 ) ) then + if ( any ( bb(3:4) == 0 ) ) then - call set_dmatrix ( d, bb(3:4) ) + call set_dmatrix ( d, bb(3:4) ) - a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - else + else - if ( lsh(2) /= gsh(2) ) then - call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) - end if + if ( lsh(2) /= gsh(2) ) then + call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) + end if - center = gsh(2) / 2 - jr = center + max(3,gsize(2)) + center = gsh(2) / 2 + jr = center + max(3,gsize(2)) - call set_dmatrix_half ( d(1:jr,1:jr), 1 ) + call set_dmatrix_half ( d(1:jr,1:jr), 1 ) - tmp(1:jr,1:jr) = -transpose ( & - matmul ( h(1:jr,1:jr), & - matmul ( & - transpose ( d(1:jr,1:jr) ), & - matmul ( b(1:jr,1:jr), & - d(1:jr,1:jr) ) ) ) ) + tmp(1:jr,1:jr) = -transpose ( & + matmul ( h(1:jr,1:jr), & + matmul ( & + transpose ( d(1:jr,1:jr) ), & + matmul ( b(1:jr,1:jr), & + d(1:jr,1:jr) ) ) ) ) - a(1:jr,1:center) = tmp(1:jr,1:center) + atmp(1:jr,1:center) = tmp(1:jr,1:center) - jl = center + 1 - max(3,gsize(2)) - d = zero - call set_dmatrix_half ( d(jl:nj,jl:nj), 2 ) - tmp(jl:nj,jl:nj) = -transpose ( & - matmul ( h(jl:nj,jl:nj), & - matmul ( & - transpose ( d(jl:nj,jl:nj) ), & - matmul ( b(jl:nj,jl:nj), & - d(jl:nj,jl:nj) ) ) ) ) + jl = center + 1 - max(3,gsize(2)) + d = zero + call set_dmatrix_half ( d(jl:nj,jl:nj), 2 ) + tmp(jl:nj,jl:nj) = -transpose ( & + matmul ( h(jl:nj,jl:nj), & + matmul ( & + transpose ( d(jl:nj,jl:nj) ), & + matmul ( b(jl:nj,jl:nj), & + d(jl:nj,jl:nj) ) ) ) ) - a(jl:nj,center+1:nj) = tmp(jl:nj,center+1:nj) + atmp(jl:nj,center+1:nj) = tmp(jl:nj,center+1:nj) + allocate ( ycoeff(patch)%coeff(nj,nj) ) + + ycoeff(patch)%coeff = atmp + + savedy(patch) = .true. + + deallocate ( atmp, d, b, h, tmp ) + + end if + end if + a => ycoeff(patch)%coeff + idel = epsilon / ( 64 * delta(2) ) - if ( bb(3) == 0 ) then jl = 1 + gsize(2) @@ -332,57 +367,69 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, & end do - deallocate ( a, d, b, h, tmp ) end if if ( zero_derivs_z == 0 ) then - allocate ( a(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk), tmp(nk,nk) ) - a = zero; d = zero; b = zero; h = zero; tmp = zero + if ( .not. savedz(patch) ) then + allocate ( atmp(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk), tmp(nk,nk) ) + atmp = zero; d = zero; b = zero; h = zero; tmp = zero - call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl ) + call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl(3) ) - call set_hmatrix ( h, bb(5:6) ) + call set_hmatrix ( h, bb(5:6) ) - if ( any ( bb(5:6) == 0 ) ) then + if ( any ( bb(5:6) == 0 ) ) then - call set_dmatrix ( d, bb(5:6) ) + call set_dmatrix ( d, bb(5:6) ) - a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - else + else - if ( lsh(3) /= gsh(3) ) then - call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) - end if + if ( lsh(3) /= gsh(3) ) then + call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) + end if - center = gsh(3) / 2 - kr = center + max(3,gsize(3)) + center = gsh(3) / 2 + kr = center + max(3,gsize(3)) - call set_dmatrix_half ( d(1:kr,1:kr), 1 ) + call set_dmatrix_half ( d(1:kr,1:kr), 1 ) - tmp(1:kr,1:kr) = -transpose ( & - matmul ( h(1:kr,1:kr), & - matmul ( & - transpose ( d(1:kr,1:kr) ), & - matmul ( b(1:kr,1:kr), & - d(1:kr,1:kr) ) ) ) ) + tmp(1:kr,1:kr) = -transpose ( & + matmul ( h(1:kr,1:kr), & + matmul ( & + transpose ( d(1:kr,1:kr) ), & + matmul ( b(1:kr,1:kr), & + d(1:kr,1:kr) ) ) ) ) - a(1:kr,1:center) = tmp(1:kr,1:center) + atmp(1:kr,1:center) = tmp(1:kr,1:center) - kl = center + 1 - max(3,gsize(3)) - d = zero - call set_dmatrix_half ( d(kl:nk,kl:nk), 2 ) - tmp(kl:nk,kl:nk) = -transpose ( & - matmul ( h(kl:nk,kl:nk), & - matmul ( & - transpose ( d(kl:nk,kl:nk) ), & - matmul ( b(kl:nk,kl:nk), & - d(kl:nk,kl:nk) ) ) ) ) + kl = center + 1 - max(3,gsize(3)) + d = zero + call set_dmatrix_half ( d(kl:nk,kl:nk), 2 ) + tmp(kl:nk,kl:nk) = -transpose ( & + matmul ( h(kl:nk,kl:nk), & + matmul ( & + transpose ( d(kl:nk,kl:nk) ), & + matmul ( b(kl:nk,kl:nk), & + d(kl:nk,kl:nk) ) ) ) ) - a(kl:nk,center+1:nk) = tmp(kl:nk,center+1:nk) + atmp(kl:nk,center+1:nk) = tmp(kl:nk,center+1:nk) + allocate ( zcoeff(patch)%coeff(nk,nk) ) + + zcoeff(patch)%coeff = atmp + + savedz(patch) = .true. + + deallocate ( atmp, d, b, h, tmp ) + + end if + end if + a => zcoeff(patch)%coeff + idel = epsilon / ( 64 * delta(3) ) if ( bb(5) == 0 ) then @@ -484,7 +531,6 @@ subroutine dissipation_4_3_alt (var, lsh, gsh, lbnd, bb, gsize, & end do - deallocate ( a, d, b, h, tmp ) end if contains |