diff options
author | diener <diener@f69c4107-0314-4c4f-9ad4-17e986b73f4a> | 2006-03-10 00:52:17 +0000 |
---|---|---|
committer | diener <diener@f69c4107-0314-4c4f-9ad4-17e986b73f4a> | 2006-03-10 00:52:17 +0000 |
commit | 3aaf16cfccc2ad7b57344ddc37fd45e69f6b7ee5 (patch) | |
tree | 6a46bbbd766659a9f52707e18d839eec460dc7bc /src/Dissipation_4_3_min_err_coeff.F90 | |
parent | 8fdb80b563f1834c4fdef6d87cc12fe81262ccb3 (diff) |
Two part update of the restricted full norm dissipation operators. The first
part is an optimisation. Previously the coefficients were constructed every
time the dissipation was applied. This turned out to have a significant
overhead, especially for large domains. Now the coefficients are constructed
only once for each patch on each processor and stored. The thorn now uses
the aliased functions MultiPatch_GetMap and MultiPatch_GetMaps (provided by
thorn MultiPatch) to figure out how many patches are present and which patch
is currently being worked on. These dissipation operators will not work unless
the MultiPatch thorn is used. Right now this should not be a problem (but can
be changed if it is), unless somebody uses dissipation operators outside of
the multipatch infrastructure.
The second part is an update aimed at making the restricted full norm
dissipation operators better behaved on grids with very different number of
grid points in different directions. For that reason the parameter:
diss_fraction
has been made vector valued (of length 3), so that the width of the transition
region from the boundary to the interior operators can be controlled
separately in the 3 dimensions. Unfortunately this means that any parameter
files, that explicitly sets this parameter, will have to be modified. I don't
think too many people are using the restricted full dissipation operators yet,
so it shouldn't cause too many problems. In additiona new parameter:
h_scaling
has been introduced. This parameter is also vector valued of length 3 and
allows the local grid spacing to be rescaled to more closely resemble the
physical grid spacing in any dimension.
Note, these changes only affect the non-diagonal dissipation operators.
git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@61 f69c4107-0314-4c4f-9ad4-17e986b73f4a
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 |