diff options
| -rw-r--r-- | interface.ccl | 6 | ||||
| -rw-r--r-- | param.ccl | 7 | ||||
| -rw-r--r-- | src/Dissipation_4_3_min_err_coeff.F90 | 100 | ||||
| -rw-r--r-- | src/Dissipation_4_3_min_err_coeff_alt.F90 | 248 | ||||
| -rw-r--r-- | src/Dissipation_6_5_min_err_coeff.F90 | 255 | ||||
| -rw-r--r-- | src/Dissipation_6_5_min_err_coeff_alt.F90 | 129 | ||||
| -rw-r--r-- | src/dissipation.c | 30 | ||||
| -rw-r--r-- | src/dissipation_coeff.F90 | 17 | ||||
| -rw-r--r-- | src/make.code.defn | 1 |
9 files changed, 512 insertions, 281 deletions
diff --git a/interface.ccl b/interface.ccl index 0c85802..e9d1e9e 100644 --- a/interface.ccl +++ b/interface.ccl @@ -51,6 +51,12 @@ CCTK_INT FUNCTION \ SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH) REQUIRES FUNCTION SymmetryTableHandleForGrid +CCTK_INT FUNCTION MultiPatch_GetMap (CCTK_POINTER_TO_CONST IN cctkGH) +USES FUNCTION MultiPatch_GetMap + +CCTK_INT FUNCTION MultiPatch_GetMaps (CCTK_POINTER_TO_CONST IN cctkGH) +USES FUNCTION MultiPatch_GetMaps + CCTK_INT FUNCTION SymmetryHandleOfName (CCTK_STRING IN sym_name) REQUIRES FUNCTION SymmetryHandleOfName @@ -38,11 +38,16 @@ REAL epsdis "Dissipation strength" STEERABLE=always *:* :: "Values typical between 0 and 1" } 0.2 -REAL diss_fraction "Fractional size of the transition region for the full restricted dissipation operator" +REAL diss_fraction[3] "Fractional size of the transition region for the full restricted dissipation operator" { 0:0.5 :: "" } 0.2 +REAL h_scaling[3] "Scaling factor for the local grid spacing in the dissipation operators" +{ + 0:* :: "Positive please" +} 1.0 + STRING vars "List of evolved grid functions that should have dissipation added" STEERABLE=always { .* :: "Must be a valid list of grid functions" diff --git a/src/Dissipation_4_3_min_err_coeff.F90 b/src/Dissipation_4_3_min_err_coeff.F90 index 1a196ac..7d61389 100644 --- a/src/Dissipation_4_3_min_err_coeff.F90 +++ b/src/Dissipation_4_3_min_err_coeff.F90 @@ -2,14 +2,18 @@ #include "cctk.h" #include "cctk_Parameters.h" +#include "cctk_Functions.h" subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, & - delta, epsilon, dfl, rhs) + delta, epsilon, dfl, npatches, patch, rhs) + + use dissipation_coeff implicit none DECLARE_CCTK_PARAMETERS + DECLARE_CCTK_FUNCTIONS CCTK_INT, dimension(3), intent(in) :: lsh, gsh, lbnd CCTK_INT :: ni, nj, nk @@ -18,28 +22,52 @@ subroutine dissipation_4_3_opt (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 - CCTK_REAL, dimension(:,:), allocatable :: a, d, b, h + CCTK_REAL, dimension(:,:), allocatable :: atmp, d, b, h + 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) ) - a = zero; d = zero; b = zero; h = 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 + + if ( .not. savedx(patch) ) then + allocate ( atmp(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni) ) + atmp = zero; d = zero; b = zero; h = zero - call set_dmatrix ( d, bb(1:2) ) + call set_dmatrix ( d, bb(1:2) ) - call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl ) + call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl(1) ) - call set_hmatrix ( h, bb(1:2) ) + call set_hmatrix ( h, bb(1:2) ) + + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + + allocate ( xcoeff(patch)%coeff(ni,ni) ) + + xcoeff(patch)%coeff = atmp - a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + savedx(patch) = .true. + + deallocate ( atmp, d, b, h ) + + end if + + a => xcoeff(patch)%coeff if ( scale_with_h > 0 ) then idel = epsilon / ( 16 * delta(1) ) @@ -125,20 +153,31 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, & a(i+2,i) * var(i+2,:,:) ) * idel end do - deallocate ( a, d, b, h ) - if ( zero_derivs_y == 0 ) then - allocate ( a(nj,nj), d(nj,nj), b(nj,nj), h(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) ) + atmp = zero; d = zero; b = zero; h = zero - call set_dmatrix ( d, bb(3:4) ) + call set_dmatrix ( d, bb(3:4) ) - 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) ) - a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + allocate ( ycoeff(patch)%coeff(nj,nj) ) + + ycoeff(patch)%coeff = atmp + + savedy(patch) = .true. + + deallocate ( atmp, d, b, h ) + + end if + + a => ycoeff(patch)%coeff + if ( scale_with_h > 0 ) then idel = epsilon / ( 16 * delta(2) ) else @@ -224,21 +263,33 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, & a(j+2,j) * var(:,j+2,:) ) * idel end do - deallocate ( a, d, b, h ) end if if ( zero_derivs_z == 0 ) then - allocate ( a(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk) ) - a = zero; d = zero; b = zero; h = zero; + if ( .not. savedz(patch) ) then + allocate ( atmp(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk) ) + atmp = zero; d = zero; b = zero; h = zero - call set_dmatrix ( d, bb(5:6) ) + call set_dmatrix ( d, bb(5:6) ) - 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) ) - a = - transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + allocate ( zcoeff(patch)%coeff(nk,nk) ) + + zcoeff(patch)%coeff = atmp + + savedz(patch) = .true. + + deallocate ( atmp, d, b, h ) + + end if + + a => zcoeff(patch)%coeff + if ( scale_with_h > 0 ) then idel = epsilon / ( 16 * delta(3) ) else @@ -324,7 +375,6 @@ subroutine dissipation_4_3_opt (var, lsh, gsh, lbnd, bb, gsize, & a(k+2,k) * var(:,:,k+2) ) * idel end do - deallocate ( a, d, b, h ) end if contains 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 diff --git a/src/Dissipation_6_5_min_err_coeff.F90 b/src/Dissipation_6_5_min_err_coeff.F90 index 63e0fbe..991d1ff 100644 --- a/src/Dissipation_6_5_min_err_coeff.F90 +++ b/src/Dissipation_6_5_min_err_coeff.F90 @@ -6,7 +6,9 @@ subroutine dissipation_6_5_opt (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_6_5_opt (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 + 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 + 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 - 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 - 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) + + end if + + allocate ( xcoeff(patch)%coeff(ni,ni) ) + + xcoeff(patch)%coeff = atmp + + savedx(patch) = .true. - a(il:ni,center+1:ni) = tmp(il:ni,center+1:ni) + deallocate ( atmp, d, b, h, tmp ) end if + a => xcoeff(patch)%coeff + if ( scale_with_h > 0 ) then idel = epsilon / ( 64 * delta(1) ) else @@ -230,56 +255,67 @@ subroutine dissipation_6_5_opt (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) ) - - if ( any ( bb(3:4) == 0 ) ) then + call set_hmatrix ( h, bb(3:4) ) - call set_dmatrix ( d, bb(3:4) ) + if ( any ( bb(3:4) == 0 ) ) then - a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + call set_dmatrix ( d, bb(3:4) ) - else - - if ( lsh(2) /= gsh(2) ) then - call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) - end if + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - center = gsh(2) / 2 - jr = center + gsize(2) + else - call set_dmatrix_half ( d(1:jr,1:jr), 1 ) + if ( lsh(2) /= gsh(2) ) then + call CCTK_WARN(0, 'Internal error. Bounding box inconsistent with grid' ) + end if - 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) ) ) ) ) + center = gsh(2) / 2 + jr = center + gsize(2) - a(1:jr,1:center) = tmp(1:jr,1:center) + call set_dmatrix_half ( d(1:jr,1:jr), 1 ) - jl = center + 1 - 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), & + tmp(1:jr,1:jr) = -transpose ( & + matmul ( h(1:jr,1:jr), & 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) + transpose ( d(1:jr,1:jr) ), & + matmul ( b(1:jr,1:jr), & + d(1:jr,1:jr) ) ) ) ) + + atmp(1:jr,1:center) = tmp(1:jr,1:center) + + jl = center + 1 - 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) ) ) ) ) + + atmp(jl:nj,center+1:nj) = tmp(jl:nj,center+1:nj) - end if + 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 + if ( scale_with_h > 0 ) then idel = epsilon / ( 64 * delta(2) ) else @@ -431,57 +467,69 @@ subroutine dissipation_6_5_opt (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 + gsize(3) + center = gsh(3) / 2 + kr = center + 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 - 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 - 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) - end if + 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 + if ( scale_with_h > 0 ) then idel = epsilon / ( 64 * delta(3) ) else @@ -633,7 +681,6 @@ subroutine dissipation_6_5_opt (var, lsh, gsh, lbnd, bb, gsize, & end do - deallocate ( a, d, b, h, tmp ) end if contains diff --git a/src/Dissipation_6_5_min_err_coeff_alt.F90 b/src/Dissipation_6_5_min_err_coeff_alt.F90 index ac9b19a..e67b2a4 100644 --- a/src/Dissipation_6_5_min_err_coeff_alt.F90 +++ b/src/Dissipation_6_5_min_err_coeff_alt.F90 @@ -6,7 +6,9 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, & - delta, epsdisl, dfl, rhs) + delta, epsilon, dfl, npatches, patch, rhs) + + use dissipation_coeff implicit none @@ -20,30 +22,53 @@ subroutine dissipation_6_5_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) :: epsdisl, 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 + 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 + + if ( .not. savedx(patch) ) then + allocate ( atmp(ni,ni), d(ni,ni), b(ni,ni), h(ni,ni) ) + atmp = zero; d = zero; b = zero; h = zero + + call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl(1) ) + + call set_hmatrix ( h, bb(1:2) ) - call set_bmatrix ( b, bb(1:2), lsh(1), gsh(1), lbnd(1), delta(1), dfl ) + call set_dmatrix ( d, bb(1:2) ) - call set_hmatrix ( h, bb(1:2) ) + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) - call set_dmatrix ( d, bb(1:2) ) + allocate ( xcoeff(patch)%coeff(ni,ni) ) + + xcoeff(patch)%coeff = atmp + + savedx(patch) = .true. + + deallocate ( atmp, d, b, h ) + + end if - a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + a => xcoeff(patch)%coeff - idel = epsdisl / ( 256.0_wp * delta(1) ) + idel = epsilon / ( 256.0_wp * delta(1) ) if ( bb(1) == 0 ) then @@ -206,28 +231,39 @@ subroutine dissipation_6_5_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) ) + atmp = zero; d = zero; b = zero; h = 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_dmatrix ( d, bb(3:4) ) + call set_hmatrix ( h, bb(3:4) ) - a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + call set_dmatrix ( d, bb(3:4) ) - idel = epsdisl / ( 256.0_wp * delta(2) ) + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + + allocate ( ycoeff(patch)%coeff(nj,nj) ) + + ycoeff(patch)%coeff = atmp + + savedy(patch) = .true. + + deallocate ( atmp, d, b, h ) + + end if + + a => ycoeff(patch)%coeff + + idel = epsilon / ( 256.0_wp * delta(2) ) - if ( bb(3) == 0 ) then - - jl = 1 + gsize(2) + if ( bb(3) == 0 ) then - else + jl = 1 + gsize(2) + else + rhs(:,1,:) = rhs(:,1,:) + & ( a(1,1) * var(:,1,:) + a(2,1) * var(:,2,:) + & a(3,1) * var(:,3,:) + a(4,1) * var(:,4,:) + & @@ -365,11 +401,11 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, & a(nj,nj) * var(:,nj,:) ) * idel jr = nj - 7 - + end if - + do j = jl, jr - + rhs(:,j,:) = rhs(:,j,:) + & ( a(j-4,j) * var(:,j-4,:) + & a(j-3,j) * var(:,j-3,:) + & @@ -383,22 +419,34 @@ subroutine dissipation_6_5_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 - - call set_bmatrix ( b, bb(5:6), lsh(3), gsh(3), lbnd(3), delta(3), dfl ) - - call set_hmatrix ( h, bb(5:6) ) - - call set_dmatrix ( d, bb(5:6) ) + if ( .not. savedz(patch) ) then + allocate ( atmp(nk,nk), d(nk,nk), b(nk,nk), h(nk,nk) ) + atmp = zero; d = zero; b = zero; h = zero - a = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + 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_dmatrix ( d, bb(5:6) ) + + atmp = -transpose ( matmul ( h, matmul ( transpose(d), matmul ( b, d ) ) ) ) + + allocate ( zcoeff(patch)%coeff(nk,nk) ) + + zcoeff(patch)%coeff = atmp + + savedz(patch) = .true. + + deallocate ( atmp, d, b, h ) + + end if - idel = epsdisl / ( 256.0_wp * delta(3) ) + a => zcoeff(patch)%coeff + + idel = epsilon / ( 256.0_wp * delta(3) ) if ( bb(5) == 0 ) then @@ -561,7 +609,6 @@ subroutine dissipation_6_5_alt (var, lsh, gsh, lbnd, bb, gsize, & end do - deallocate ( a, d, b, h, tmp ) end if contains @@ -605,7 +652,7 @@ contains implicit none - CCTK_REAL, dimension(:,:), intent(out) :: b + CCTK_REAL, dimension(:,:), intent(inout) :: b CCTK_INT, dimension(2), intent(in) :: bb CCTK_INT, intent(in) :: lsh, gsh, lbnd CCTK_REAL, intent(in) :: h, dfl diff --git a/src/dissipation.c b/src/dissipation.c index 48f7ba4..8f85d32 100644 --- a/src/dissipation.c +++ b/src/dissipation.c @@ -6,6 +6,7 @@ #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" +#include "cctk_Functions.h" #include "util_ErrorCodes.h" #include "util_Table.h" @@ -19,6 +20,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_4_3_opt) (const CCTK_REAL *var, const CCTK_REAL *dx, const CCTK_REAL *epsdis, const CCTK_REAL *diss_fraction, + const CCTK_INT *npatches, + const CCTK_INT *patch, CCTK_REAL *rhs); void CCTK_FCALL CCTK_FNAME(dissipation_4_3_alt) (const CCTK_REAL *var, const CCTK_INT *lsh, @@ -29,6 +32,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_4_3_alt) (const CCTK_REAL *var, const CCTK_REAL *dx, const CCTK_REAL *epsdis, const CCTK_REAL *diss_fraction, + const CCTK_INT *npatches, + const CCTK_INT *patch, CCTK_REAL *rhs); void CCTK_FCALL CCTK_FNAME(dissipation_6_5_opt) (const CCTK_REAL *var, const CCTK_INT *ni, @@ -39,6 +44,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_6_5_opt) (const CCTK_REAL *var, const CCTK_REAL *dx, const CCTK_REAL *epsdis, const CCTK_REAL *diss_fraction, + const CCTK_INT *npatches, + const CCTK_INT *patch, CCTK_REAL *rhs); void CCTK_FCALL CCTK_FNAME(dissipation_6_5_alt) (const CCTK_REAL *var, const CCTK_INT *ni, @@ -49,6 +56,8 @@ void CCTK_FCALL CCTK_FNAME(dissipation_6_5_alt) (const CCTK_REAL *var, const CCTK_REAL *dx, const CCTK_REAL *epsdis, const CCTK_REAL *diss_fraction, + const CCTK_INT *npatches, + const CCTK_INT *patch, CCTK_REAL *rhs); void CCTK_FCALL CCTK_FNAME(dissipation_2_1) (const CCTK_REAL *var, const CCTK_INT *ni, @@ -153,11 +162,12 @@ apply (int const varindex, char const * const optstring, void * const arg) int ierr; CCTK_INT symtable, pen_sym_handle; CCTK_INT symbnd[6]; + CCTK_INT npatches, patch; assert (varindex >= 0); for (d=0; d<3; ++d) { - dx[d] = CCTK_DELTA_SPACE(d); + dx[d] = CCTK_DELTA_SPACE(d)*h_scaling[d]; gsize[d] = cctk_nghostzones[d]; } @@ -280,28 +290,30 @@ apply (int const varindex, char const * const optstring, void * const arg) assert(0); } } else { + npatches = MultiPatch_GetMaps ( cctkGH ); + patch = MultiPatch_GetMap ( cctkGH ) + 1; switch(order) { case 4: { if ( CCTK_Equals(dissipation_type,"Mattson-Svard-Nordstrom") ) { CCTK_FNAME(dissipation_4_3_opt) - (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, - bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr); + (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize, + dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr); } else { CCTK_FNAME(dissipation_4_3_alt) - (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, - bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr); + (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize, + dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr); } break; } case 6: { if ( CCTK_Equals(dissipation_type,"Mattson-Svard-Nordstrom") ) { CCTK_FNAME(dissipation_6_5_opt) - (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, - bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr); + (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize, + dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr); } else { CCTK_FNAME(dissipation_6_5_alt) - (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, - bbox, gsize, dx, &epsdis, &diss_fraction, rhsptr); + (varptr, cctk_lsh, cctk_gsh, cctk_lbnd, bbox, gsize, + dx, &epsdis, diss_fraction, &npatches, &patch, rhsptr); } break; } diff --git a/src/dissipation_coeff.F90 b/src/dissipation_coeff.F90 new file mode 100644 index 0000000..dfd3d94 --- /dev/null +++ b/src/dissipation_coeff.F90 @@ -0,0 +1,17 @@ +! Arrays to store dissipation operator coefficients for the full restricted norm case. +! $Header$ + +#include "cctk.h" + +module dissipation_coeff + + type coefficients + CCTK_REAL, dimension(:,:), pointer :: coeff + end type coefficients + + type(coefficients), dimension(:), pointer :: xcoeff, ycoeff, zcoeff + + logical :: first = .true. + logical, dimension(:), allocatable :: savedx, savedy, savedz + +end module dissipation_coeff diff --git a/src/make.code.defn b/src/make.code.defn index 1519773..0371d01 100644 --- a/src/make.code.defn +++ b/src/make.code.defn @@ -31,6 +31,7 @@ SRCS = call_derivs.c \ Dissipation_4_3_min_err_coeff_alt.F90 \ Dissipation_6_5_min_err_coeff.F90 \ Dissipation_6_5_min_err_coeff_alt.F90 \ + dissipation_coeff.F90 \ get_coeffs.c \ Coefficients_2_1.F90 \ Coefficients_4_2.F90 \ |