diff options
Diffstat (limited to 'src/Dissipation_4_3_min_err_coeff.F90')
| -rw-r--r-- | src/Dissipation_4_3_min_err_coeff.F90 | 100 |
1 files changed, 75 insertions, 25 deletions
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 |