aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authordiener <diener@f69c4107-0314-4c4f-9ad4-17e986b73f4a>2008-07-14 22:37:12 +0000
committerdiener <diener@f69c4107-0314-4c4f-9ad4-17e986b73f4a>2008-07-14 22:37:12 +0000
commit6af885f4ff90d496a972d2d8177479512dcb4cbe (patch)
treebef0ba5ed06f76fdb8d855f3bf5f8bcbe4348872
parente3b9cd6b30215fb09a45910eafca7cd9121ceffa (diff)
Added functionality to set up public grid functions (sbp_dx, sbp_dy, sbp_dz)
containing the delta's used in the Kreiss-Oliger type dissipation. This is activated with a parameter (use_variable_deltas) and a routine (SBP_DeltaInitial) is scheduled at BASEGRID to set the default values initially. The grid functions can be modified later by a scheduled user function if so desired. git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/SummationByParts/trunk@116 f69c4107-0314-4c4f-9ad4-17e986b73f4a
Diffstat
-rw-r--r--interface.ccl5
-rw-r--r--param.ccl4
-rw-r--r--schedule.ccl15
-rw-r--r--src/DeltaInitial.F9024
-rw-r--r--src/Dissipation_2_1_alt.F90224
-rw-r--r--src/Dissipation_4_2_alt.F90297
-rw-r--r--src/Dissipation_6_3_alt.F90438
-rw-r--r--src/Dissipation_8_4_alt.F90595
-rw-r--r--src/dissipation.c100
-rw-r--r--src/make.code.defn1
10 files changed, 1690 insertions, 13 deletions
diff --git a/interface.ccl b/interface.ccl
index bc6c63b..15fc072 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -136,3 +136,8 @@ CCTK_REAL normmask TYPE=gf TAGS='tensortypealias="scalar" Prolongation="None" ch
{
nmask
} "Mask for the norm calculation"
+
+CCTK_REAL deltas TYPE=gf TIMELEVELS=1 TAGS='tensortypealias="U" Prolongation="None"'
+{
+ sbp_dx, sbp_dy, sbp_dz
+} "Dissipation deltas"
diff --git a/param.ccl b/param.ccl
index eba0e40..583a629 100644
--- a/param.ccl
+++ b/param.ccl
@@ -38,6 +38,10 @@ BOOLEAN use_dissipation "Should we add dissipation"
{
} "no"
+BOOLEAN use_variable_deltas "Use extra grid functions to allow for variable delta's in the dissipation operators"
+{
+} "no"
+
BOOLEAN use_shiftout "Should we use the boundary_shift_out parameters from CoordBase to shift the stencils of derivatives and dissipation"
{
} "no"
diff --git a/schedule.ccl b/schedule.ccl
index e4c492d..3fd0600 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -4,10 +4,23 @@
# Mask for norms
STORAGE: normmask
+if (use_variable_deltas)
+{
+ STORAGE: deltas
+}
+
SCHEDULE SBP_SetNormMask AT basegrid
{
LANG: Fortran
-} "Setup the mask for the calulcation of the norm"
+} "Setup the mask for the calculation of the norm"
+
+if (use_variable_deltas)
+{
+ SCHEDULE SBP_DeltaInitial AT basegrid
+ {
+ LANG: Fortran
+ } "Initialize dissipation deltas"
+}
if (check_grid_sizes)
{
diff --git a/src/DeltaInitial.F90 b/src/DeltaInitial.F90
new file mode 100644
index 0000000..97b43b3
--- /dev/null
+++ b/src/DeltaInitial.F90
@@ -0,0 +1,24 @@
+! Check if there are enough grid cells and ghost zones to compute the finite
+! differences.
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+#include "cctk_Functions.h"
+
+subroutine SBP_DeltaInitial(CCTK_ARGUMENTS)
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_ARGUMENTS
+ DECLARE_CCTK_FUNCTIONS
+
+!$OMP PARALLEL WORKSHARE
+ sbp_dx = CCTK_DELTA_SPACE(1)*h_scaling(1)
+ sbp_dy = CCTK_DELTA_SPACE(2)*h_scaling(2)
+ sbp_dz = CCTK_DELTA_SPACE(3)*h_scaling(3)
+!$OMP END PARALLEL WORKSHARE
+
+end subroutine SBP_DeltaInitial
diff --git a/src/Dissipation_2_1_alt.F90 b/src/Dissipation_2_1_alt.F90
index 3501cdf..0c0cda5 100644
--- a/src/Dissipation_2_1_alt.F90
+++ b/src/Dissipation_2_1_alt.F90
@@ -33,12 +33,15 @@ subroutine dissipation_2_1_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
il = 1 + gsize(1)
else
ol = offset(1)
+
+!$OMP PARALLEL WORKSHARE
rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
a(3,1) * var(3+ol,:,:) ) * idel
rhs(2+ol,:,:) = rhs(2+ol,:,:) + &
( a(1,2) * var(1+ol,:,:) + a(2,2) * var(2+ol,:,:) + &
a(3,2) * var(3+ol,:,:) + a(4,2) * var(4+ol,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
il = 3 + ol
end if
@@ -46,21 +49,27 @@ subroutine dissipation_2_1_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
ir = ni - gsize(1)
else
or = ni - offset(2)
+
+!$OMP PARALLEL WORKSHARE
rhs(or-1,:,:) = rhs(or-1,:,:) + &
( a(1,2) * var(or,:,:) + a(2,2) * var(or-1,:,:) + &
a(3,2) * var(or-2,:,:) + a(4,2) * var(or-3,:,:) ) * idel
rhs(or,:,:) = rhs(or,:,:) + &
( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
a(3,1) * var(or-2,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
ir = or - 2
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
( -6.0_wp * var(il:ir,:,:) + &
4.0_wp * ( var(il-1:ir-1,:,:) + &
var(il+1:ir+1,:,:) ) - &
( var(il-2:ir-2,:,:) + &
var(il+2:ir+2,:,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
if ( zero_derivs_y == 0 ) then
call set_coeff ( a )
@@ -71,12 +80,15 @@ subroutine dissipation_2_1_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jl = 1 + gsize(2)
else
ol = offset(3)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
a(3,1) * var(:,3+ol,:) ) * idel
rhs(:,2+ol,:) = rhs(:,2+ol,:) + &
( a(1,2) * var(:,1+ol,:) + a(2,2) * var(:,2+ol,:) + &
a(3,2) * var(:,3+ol,:) + a(4,2) * var(:,4+ol,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jl = 3 + ol
end if
@@ -84,21 +96,27 @@ subroutine dissipation_2_1_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jr = nj - gsize(2)
else
or = nj - offset(4)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,or-1,:) = rhs(:,or-1,:) + &
( a(1,2) * var(:,or,:) + a(2,2) * var(:,or-1,:) + &
a(3,2) * var(:,or-2,:) + a(4,2) * var(:,or-3,:) ) * idel
rhs(:,or,:) = rhs(:,or,:) + &
( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
a(3,1) * var(:,or-2,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jr = or - 2
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
( -6.0_wp * var(:,jl:jr,:) + &
4.0_wp * ( var(:,jl-1:jr-1,:) + &
var(:,jl+1:jr+1,:) ) - &
( var(:,jl-2:jr-2,:) + &
var(:,jl+2:jr+2,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
if ( zero_derivs_z == 0 ) then
@@ -110,12 +128,15 @@ subroutine dissipation_2_1_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kl = 1 + gsize(3)
else
ol = offset(5)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
a(3,1) * var(:,:,3+ol) ) * idel
rhs(:,:,2+ol) = rhs(:,:,2+ol) + &
( a(1,2) * var(:,:,1+ol) + a(2,2) * var(:,:,2+ol) + &
a(3,2) * var(:,:,3+ol) + a(4,2) * var(:,:,4+ol) ) * idel
+!$OMP END PARALLEL WORKSHARE
kl = 3 + ol
end if
@@ -123,21 +144,27 @@ subroutine dissipation_2_1_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kr = nk - gsize(3)
else
or = nk - offset(6)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,or-1) = rhs(:,:,or-1) + &
( a(1,2) * var(:,:,or) + a(2,2) * var(:,:,or-1) + &
a(3,2) * var(:,:,or-2) + a(4,2) * var(:,:,or-3) ) * idel
rhs(:,:,or) = rhs(:,:,or) + &
( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
a(3,1) * var(:,:,or-2) ) * idel
+!$OMP END PARALLEL WORKSHARE
kr = or - 2
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
( -6.0_wp * var(:,:,kl:kr) + &
4.0_wp * ( var(:,:,kl-1:kr-1) + &
var(:,:,kl+1:kr+1) ) - &
( var(:,:,kl-2:kr-2) + &
var(:,:,kl+2:kr+2) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
contains
@@ -162,3 +189,200 @@ contains
end subroutine set_coeff
end subroutine dissipation_2_1_alt
+
+subroutine dissipation_2_1_delta (var, ni, nj, nk, bb, gsize, offset, &
+ dx, dy, dz, epsilon, rhs)
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: ni, nj, nk
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
+ CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
+ CCTK_INT, dimension(6), intent(in) :: bb
+ CCTK_INT, dimension(3), intent(in) :: gsize
+ CCTK_INT, dimension(6), intent(in) :: offset
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: dx, dy, dz
+ CCTK_REAL, intent(in) :: epsilon
+
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+ CCTK_REAL, dimension(4,2) :: a
+ CCTK_REAL :: idel
+
+ CCTK_INT :: il, ir, jl, jr, kl, kr, ol, or
+
+ call set_coeff ( a )
+
+ idel = epsilon / 16
+
+ if ( bb(1) == 0 ) then
+ il = 1 + gsize(1)
+ else
+ ol = offset(1)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
+ ( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
+ a(3,1) * var(3+ol,:,:) ) * idel / dx(1+ol,:,:)
+ rhs(2+ol,:,:) = rhs(2+ol,:,:) + &
+ ( a(1,2) * var(1+ol,:,:) + a(2,2) * var(2+ol,:,:) + &
+ a(3,2) * var(3+ol,:,:) + a(4,2) * var(4+ol,:,:) ) * &
+ idel / dx(2+ol,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ il = 3 + ol
+ end if
+ if ( bb(2) == 0 ) then
+ ir = ni - gsize(1)
+ else
+ or = ni - offset(2)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(or-1,:,:) = rhs(or-1,:,:) + &
+ ( a(1,2) * var(or,:,:) + a(2,2) * var(or-1,:,:) + &
+ a(3,2) * var(or-2,:,:) + a(4,2) * var(or-3,:,:) ) * &
+ idel / dx(or-1,:,:)
+ rhs(or,:,:) = rhs(or,:,:) + &
+ ( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
+ a(3,1) * var(or-2,:,:) ) * idel / dx(or,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ ir = or - 2
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+ ( -6.0_wp * var(il:ir,:,:) + &
+ 4.0_wp * ( var(il-1:ir-1,:,:) + &
+ var(il+1:ir+1,:,:) ) - &
+ ( var(il-2:ir-2,:,:) + &
+ var(il+2:ir+2,:,:) ) ) * &
+ idel / dx(il:ir,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ if ( zero_derivs_y == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 16
+
+ if ( bb(3) == 0 ) then
+ jl = 1 + gsize(2)
+ else
+ ol = offset(3)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
+ ( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
+ a(3,1) * var(:,3+ol,:) ) * idel / dy(:,1+ol,:)
+ rhs(:,2+ol,:) = rhs(:,2+ol,:) + &
+ ( a(1,2) * var(:,1+ol,:) + a(2,2) * var(:,2+ol,:) + &
+ a(3,2) * var(:,3+ol,:) + a(4,2) * var(:,4+ol,:) ) * &
+ idel /dy(:,2+ol,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jl = 3 + ol
+ end if
+ if ( bb(4) == 0 ) then
+ jr = nj - gsize(2)
+ else
+ or = nj - offset(4)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,or-1,:) = rhs(:,or-1,:) + &
+ ( a(1,2) * var(:,or,:) + a(2,2) * var(:,or-1,:) + &
+ a(3,2) * var(:,or-2,:) + a(4,2) * var(:,or-3,:) ) * &
+ idel / dy(:,or-1,:)
+ rhs(:,or,:) = rhs(:,or,:) + &
+ ( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
+ a(3,1) * var(:,or-2,:) ) * idel / dy(:,or,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jr = or - 2
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+ ( -6.0_wp * var(:,jl:jr,:) + &
+ 4.0_wp * ( var(:,jl-1:jr-1,:) + &
+ var(:,jl+1:jr+1,:) ) - &
+ ( var(:,jl-2:jr-2,:) + &
+ var(:,jl+2:jr+2,:) ) ) * &
+ idel / dy(:,jl:jr,:)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+ if ( zero_derivs_z == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 16
+
+ if ( bb(5) == 0 ) then
+ kl = 1 + gsize(3)
+ else
+ ol = offset(5)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
+ ( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
+ a(3,1) * var(:,:,3+ol) ) * idel / dz (:,:,1+ol)
+ rhs(:,:,2+ol) = rhs(:,:,2+ol) + &
+ ( a(1,2) * var(:,:,1+ol) + a(2,2) * var(:,:,2+ol) + &
+ a(3,2) * var(:,:,3+ol) + a(4,2) * var(:,:,4+ol) ) * &
+ idel / dz(:,:,2+ol)
+!$OMP END PARALLEL WORKSHARE
+
+ kl = 3 + ol
+ end if
+ if ( bb(6) == 0 ) then
+ kr = nk - gsize(3)
+ else
+ or = nk - offset(6)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,or-1) = rhs(:,:,or-1) + &
+ ( a(1,2) * var(:,:,or) + a(2,2) * var(:,:,or-1) + &
+ a(3,2) * var(:,:,or-2) + a(4,2) * var(:,:,or-3) ) * &
+ idel / dz(:,:,or-1)
+ rhs(:,:,or) = rhs(:,:,or) + &
+ ( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
+ a(3,1) * var(:,:,or-2) ) * idel / dz(:,:,or)
+!$OMP END PARALLEL WORKSHARE
+
+ kr = or - 2
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+ ( -6.0_wp * var(:,:,kl:kr) + &
+ 4.0_wp * ( var(:,:,kl-1:kr-1) + &
+ var(:,:,kl+1:kr+1) ) - &
+ ( var(:,:,kl-2:kr-2) + &
+ var(:,:,kl+2:kr+2) ) ) * &
+ idel / dz(:,:,kl:kr)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+contains
+
+ subroutine set_coeff ( a )
+
+ implicit none
+
+ CCTK_REAL, dimension(4,2), intent(OUT) :: a
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+
+ a(1,1) = -2.0_wp
+ a(2,1) = 4.0_wp
+ a(3,1) = -2.0_wp
+ a(4,1) = 0.0_wp
+ a(1,2) = 2.0_wp
+ a(2,2) = -5.0_wp
+ a(3,2) = 4.0_wp
+ a(4,2) = -1.0_wp
+
+ end subroutine set_coeff
+
+end subroutine dissipation_2_1_delta
diff --git a/src/Dissipation_4_2_alt.F90 b/src/Dissipation_4_2_alt.F90
index 577ed39..7a615d6 100644
--- a/src/Dissipation_4_2_alt.F90
+++ b/src/Dissipation_4_2_alt.F90
@@ -33,6 +33,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
il = 1 + gsize(1)
else
ol = offset(1)
+!$OMP PARALLEL WORKSHARE
rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
a(3,1) * var(3+ol,:,:) + a(4,1) * var(4+ol,:,:) ) * idel
@@ -49,6 +50,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(3,4) * var(3+ol,:,:) + a(4,4) * var(4+ol,:,:) + &
a(5,4) * var(5+ol,:,:) + a(6,4) * var(6+ol,:,:) + &
a(7,4) * var(7+ol,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
il = 5 + ol
end if
@@ -56,6 +58,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
ir = ni - gsize(1)
else
or = ni - offset(2)
+!$OMP PARALLEL WORKSHARE
rhs(or-3,:,:) = rhs(or-3,:,:) + &
( a(1,4) * var(or,:,:) + a(2,4) * var(or-1,:,:) + &
a(3,4) * var(or-2,:,:) + a(4,4) * var(or-3,:,:) + &
@@ -72,9 +75,11 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
rhs(or,:,:) = rhs(or,:,:) + &
( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
a(3,1) * var(or-2,:,:) + a(4,1) * var(or-3,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
ir = or - 4
end if
+!$OMP PARALLEL WORKSHARE
rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
( -20.0_wp * var(il:ir,:,:) + &
15.0_wp * ( var(il-1:ir-1,:,:) + &
@@ -83,6 +88,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(il+2:ir+2,:,:) ) + &
( var(il-3:ir-3,:,:) + &
var(il+3:ir+3,:,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
if ( zero_derivs_y == 0 ) then
call set_coeff ( a )
@@ -93,6 +99,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jl = 1 + gsize(2)
else
ol = offset(3)
+!$OMP PARALLEL WORKSHARE
rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
a(3,1) * var(:,3+ol,:) + a(4,1) * var(:,4+ol,:) ) * idel
@@ -109,6 +116,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(3,4) * var(:,3+ol,:) + a(4,4) * var(:,4+ol,:) + &
a(5,4) * var(:,5+ol,:) + a(6,4) * var(:,6+ol,:) + &
a(7,4) * var(:,7+ol,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jl = 5 + ol
end if
@@ -116,6 +124,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jr = nj - gsize(2)
else
or = nj - offset(4)
+!$OMP PARALLEL WORKSHARE
rhs(:,or-3,:) = rhs(:,or-3,:) + &
( a(1,4) * var(:,or,:) + a(2,4) * var(:,or-1,:) + &
a(3,4) * var(:,or-2,:) + a(4,4) * var(:,or-3,:) + &
@@ -132,9 +141,11 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
rhs(:,or,:) = rhs(:,or,:) + &
( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
a(3,1) * var(:,or-2,:) + a(4,1) * var(:,or-3,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jr = or - 4
end if
+!$OMP PARALLEL WORKSHARE
rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
( -20.0_wp * var(:,jl:jr,:) + &
15.0_wp * ( var(:,jl-1:jr-1,:) + &
@@ -143,6 +154,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(:,jl+2:jr+2,:) ) + &
( var(:,jl-3:jr-3,:) + &
var(:,jl+3:jr+3,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
if ( zero_derivs_z == 0 ) then
@@ -154,6 +166,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kl = 1 + gsize(3)
else
ol = offset(5)
+!$OMP PARALLEL WORKSHARE
rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
a(3,1) * var(:,:,3+ol) + a(4,1) * var(:,:,4+ol) ) * idel
@@ -170,6 +183,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(3,4) * var(:,:,3+ol) + a(4,4) * var(:,:,4+ol) + &
a(5,4) * var(:,:,5+ol) + a(6,4) * var(:,:,6+ol) + &
a(7,4) * var(:,:,7+ol) ) * idel
+!$OMP END PARALLEL WORKSHARE
kl = 5 + ol
end if
@@ -177,6 +191,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kr = nk - gsize(3)
else
or = nk - offset(6)
+!$OMP PARALLEL WORKSHARE
rhs(:,:,or-3) = rhs(:,:,or-3) + &
( a(1,4) * var(:,:,or) + a(2,4) * var(:,:,or-1) + &
a(3,4) * var(:,:,or-2) + a(4,4) * var(:,:,or-3) + &
@@ -193,9 +208,11 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
rhs(:,:,or) = rhs(:,:,or) + &
( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
a(3,1) * var(:,:,or-2) + a(4,1) * var(:,:,or-3) ) * idel
+!$OMP END PARALLEL WORKSHARE
kr = or - 4
end if
+!$OMP PARALLEL WORKSHARE
rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
( -20.0_wp * var(:,:,kl:kr) + &
15.0_wp * ( var(:,:,kl-1:kr-1) + &
@@ -204,6 +221,7 @@ subroutine dissipation_4_2_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(:,:,kl+2:kr+2) ) + &
( var(:,:,kl-3:kr-3) + &
var(:,:,kl+3:kr+3) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
contains
@@ -248,3 +266,282 @@ contains
end subroutine set_coeff
end subroutine dissipation_4_2_alt
+
+subroutine dissipation_4_2_delta (var, ni, nj, nk, bb, gsize, offset, &
+ dx, dy, dz, epsilon, rhs)
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: ni, nj, nk
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
+ CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
+ CCTK_INT, dimension(6), intent(in) :: bb
+ CCTK_INT, dimension(3), intent(in) :: gsize
+ CCTK_INT, dimension(6), intent(in) :: offset
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: dx, dy, dz
+ CCTK_REAL, intent(in) :: epsilon
+
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+ CCTK_REAL, dimension(7,4) :: a
+ CCTK_REAL :: idel
+
+ CCTK_INT :: il, ir, jl, jr, kl, kr, ol, or
+
+ call set_coeff ( a )
+
+ idel = epsilon / 64
+
+ if ( bb(1) == 0 ) then
+ il = 1 + gsize(1)
+ else
+ ol = offset(1)
+!$OMP PARALLEL WORKSHARE
+ rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
+ ( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
+ a(3,1) * var(3+ol,:,:) + a(4,1) * var(4+ol,:,:) ) * &
+ idel / dx(1+ol,:,:)
+ rhs(2+ol,:,:) = rhs(2+ol,:,:) + &
+ ( a(1,2) * var(1+ol,:,:) + a(2,2) * var(2+ol,:,:) + &
+ a(3,2) * var(3+ol,:,:) + a(4,2) * var(4+ol,:,:) + &
+ a(5,2) * var(5+ol,:,:) ) * idel / dx(2+ol,:,:)
+ rhs(3+ol,:,:) = rhs(3+ol,:,:) + &
+ ( a(1,3) * var(1+ol,:,:) + a(2,3) * var(2+ol,:,:) + &
+ a(3,3) * var(3+ol,:,:) + a(4,3) * var(4+ol,:,:) + &
+ a(5,3) * var(5+ol,:,:) + a(6,3) * var(6+ol,:,:) ) * &
+ idel / dx(3+ol,:,:)
+ rhs(4+ol,:,:) = rhs(4+ol,:,:) + &
+ ( a(1,4) * var(1+ol,:,:) + a(2,4) * var(2+ol,:,:) + &
+ a(3,4) * var(3+ol,:,:) + a(4,4) * var(4+ol,:,:) + &
+ a(5,4) * var(5+ol,:,:) + a(6,4) * var(6+ol,:,:) + &
+ a(7,4) * var(7+ol,:,:) ) * idel / dx(4+ol,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ il = 5 + ol
+ end if
+ if ( bb(2) == 0 ) then
+ ir = ni - gsize(1)
+ else
+ or = ni - offset(2)
+!$OMP PARALLEL WORKSHARE
+ rhs(or-3,:,:) = rhs(or-3,:,:) + &
+ ( a(1,4) * var(or,:,:) + a(2,4) * var(or-1,:,:) + &
+ a(3,4) * var(or-2,:,:) + a(4,4) * var(or-3,:,:) + &
+ a(5,4) * var(or-4,:,:) + a(6,4) * var(or-5,:,:) + &
+ a(7,4) * var(or-6,:,:) ) * idel / dx(or-3,:,:)
+ rhs(or-2,:,:) = rhs(or-2,:,:) + &
+ ( a(1,3) * var(or,:,:) + a(2,3) * var(or-1,:,:) + &
+ a(3,3) * var(or-2,:,:) + a(4,3) * var(or-3,:,:) + &
+ a(5,3) * var(or-4,:,:) + a(6,3) * var(or-5,:,:) ) * &
+ idel / dx(or-2,:,:)
+ rhs(or-1,:,:) = rhs(or-1,:,:) + &
+ ( a(1,2) * var(or,:,:) + a(2,2) * var(or-1,:,:) + &
+ a(3,2) * var(or-2,:,:) + a(4,2) * var(or-3,:,:) + &
+ a(5,2) * var(or-4,:,:) ) * idel / dx(or-1,:,:)
+ rhs(or,:,:) = rhs(or,:,:) + &
+ ( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
+ a(3,1) * var(or-2,:,:) + a(4,1) * var(or-3,:,:) ) * &
+ idel / dx(or,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ ir = or - 4
+ end if
+!$OMP PARALLEL WORKSHARE
+ rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+ ( -20.0_wp * var(il:ir,:,:) + &
+ 15.0_wp * ( var(il-1:ir-1,:,:) + &
+ var(il+1:ir+1,:,:) ) - &
+ 6.0_wp * ( var(il-2:ir-2,:,:) + &
+ var(il+2:ir+2,:,:) ) + &
+ ( var(il-3:ir-3,:,:) + &
+ var(il+3:ir+3,:,:) ) ) * idel / dx(il:ir,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ if ( zero_derivs_y == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 64
+
+ if ( bb(3) == 0 ) then
+ jl = 1 + gsize(2)
+ else
+ ol = offset(3)
+!$OMP PARALLEL WORKSHARE
+ rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
+ ( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
+ a(3,1) * var(:,3+ol,:) + a(4,1) * var(:,4+ol,:) ) * &
+ idel / dy(:,1+ol,:)
+ rhs(:,2+ol,:) = rhs(:,2+ol,:) + &
+ ( a(1,2) * var(:,1+ol,:) + a(2,2) * var(:,2+ol,:) + &
+ a(3,2) * var(:,3+ol,:) + a(4,2) * var(:,4+ol,:) + &
+ a(5,2) * var(:,5+ol,:) ) * idel / dy(:,2+ol,:)
+ rhs(:,3+ol,:) = rhs(:,3+ol,:) + &
+ ( a(1,3) * var(:,1+ol,:) + a(2,3) * var(:,2+ol,:) + &
+ a(3,3) * var(:,3+ol,:) + a(4,3) * var(:,4+ol,:) + &
+ a(5,3) * var(:,5+ol,:) + a(6,3) * var(:,6+ol,:) ) * &
+ idel / dy(:,3+ol,:)
+ rhs(:,4+ol,:) = rhs(:,4+ol,:) + &
+ ( a(1,4) * var(:,1+ol,:) + a(2,4) * var(:,2+ol,:) + &
+ a(3,4) * var(:,3+ol,:) + a(4,4) * var(:,4+ol,:) + &
+ a(5,4) * var(:,5+ol,:) + a(6,4) * var(:,6+ol,:) + &
+ a(7,4) * var(:,7+ol,:) ) * idel / dy(:,4+ol,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jl = 5 + ol
+ end if
+ if ( bb(4) == 0 ) then
+ jr = nj - gsize(2)
+ else
+ or = nj - offset(4)
+!$OMP PARALLEL WORKSHARE
+ rhs(:,or-3,:) = rhs(:,or-3,:) + &
+ ( a(1,4) * var(:,or,:) + a(2,4) * var(:,or-1,:) + &
+ a(3,4) * var(:,or-2,:) + a(4,4) * var(:,or-3,:) + &
+ a(5,4) * var(:,or-4,:) + a(6,4) * var(:,or-5,:) + &
+ a(7,4) * var(:,or-6,:) ) * idel / dy(:,or-3,:)
+ rhs(:,or-2,:) = rhs(:,or-2,:) + &
+ ( a(1,3) * var(:,or,:) + a(2,3) * var(:,or-1,:) + &
+ a(3,3) * var(:,or-2,:) + a(4,3) * var(:,or-3,:) + &
+ a(5,3) * var(:,or-4,:) + a(6,3) * var(:,or-5,:) ) * &
+ idel / dy(:,or-2,:)
+ rhs(:,or-1,:) = rhs(:,or-1,:) + &
+ ( a(1,2) * var(:,or,:) + a(2,2) * var(:,or-1,:) + &
+ a(3,2) * var(:,or-2,:) + a(4,2) * var(:,or-3,:) + &
+ a(5,2) * var(:,or-4,:) ) * idel / dy(:,or-1,:)
+ rhs(:,or,:) = rhs(:,or,:) + &
+ ( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
+ a(3,1) * var(:,or-2,:) + a(4,1) * var(:,or-3,:) ) * &
+ idel / dy(:,or,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jr = or - 4
+ end if
+!$OMP PARALLEL WORKSHARE
+ rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+ ( -20.0_wp * var(:,jl:jr,:) + &
+ 15.0_wp * ( var(:,jl-1:jr-1,:) + &
+ var(:,jl+1:jr+1,:) ) - &
+ 6.0_wp * ( var(:,jl-2:jr-2,:) + &
+ var(:,jl+2:jr+2,:) ) + &
+ ( var(:,jl-3:jr-3,:) + &
+ var(:,jl+3:jr+3,:) ) ) * &
+ idel / dy(:,jl:jr,:)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+ if ( zero_derivs_z == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 64
+
+ if ( bb(5) == 0 ) then
+ kl = 1 + gsize(3)
+ else
+ ol = offset(5)
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
+ ( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
+ a(3,1) * var(:,:,3+ol) + a(4,1) * var(:,:,4+ol) ) * &
+ idel / dz(:,:,1+ol)
+ rhs(:,:,2+ol) = rhs(:,:,2+ol) + &
+ ( a(1,2) * var(:,:,1+ol) + a(2,2) * var(:,:,2+ol) + &
+ a(3,2) * var(:,:,3+ol) + a(4,2) * var(:,:,4+ol) + &
+ a(5,2) * var(:,:,5+ol) ) * idel / dz(:,:,2+ol)
+ rhs(:,:,3+ol) = rhs(:,:,3+ol) + &
+ ( a(1,3) * var(:,:,1+ol) + a(2,3) * var(:,:,2+ol) + &
+ a(3,3) * var(:,:,3+ol) + a(4,3) * var(:,:,4+ol) + &
+ a(5,3) * var(:,:,5+ol) + a(6,3) * var(:,:,6+ol) ) * &
+ idel / dz(:,:,3+ol)
+ rhs(:,:,4+ol) = rhs(:,:,4+ol) + &
+ ( a(1,4) * var(:,:,1+ol) + a(2,4) * var(:,:,2+ol) + &
+ a(3,4) * var(:,:,3+ol) + a(4,4) * var(:,:,4+ol) + &
+ a(5,4) * var(:,:,5+ol) + a(6,4) * var(:,:,6+ol) + &
+ a(7,4) * var(:,:,7+ol) ) * idel / dz(:,:,4+ol)
+!$OMP END PARALLEL WORKSHARE
+
+ kl = 5 + ol
+ end if
+ if ( bb(6) == 0 ) then
+ kr = nk - gsize(3)
+ else
+ or = nk - offset(6)
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,or-3) = rhs(:,:,or-3) + &
+ ( a(1,4) * var(:,:,or) + a(2,4) * var(:,:,or-1) + &
+ a(3,4) * var(:,:,or-2) + a(4,4) * var(:,:,or-3) + &
+ a(5,4) * var(:,:,or-4) + a(6,4) * var(:,:,or-5) + &
+ a(7,4) * var(:,:,or-6) ) * idel / dz(:,:,or-3)
+ rhs(:,:,or-2) = rhs(:,:,or-2) + &
+ ( a(1,3) * var(:,:,or) + a(2,3) * var(:,:,or-1) + &
+ a(3,3) * var(:,:,or-2) + a(4,3) * var(:,:,or-3) + &
+ a(5,3) * var(:,:,or-4) + a(6,3) * var(:,:,or-5) ) * &
+ idel / dz(:,:,or-2)
+ rhs(:,:,or-1) = rhs(:,:,or-1) + &
+ ( a(1,2) * var(:,:,or) + a(2,2) * var(:,:,or-1) + &
+ a(3,2) * var(:,:,or-2) + a(4,2) * var(:,:,or-3) + &
+ a(5,2) * var(:,:,or-4) ) * idel / dz(:,:,or-1)
+ rhs(:,:,or) = rhs(:,:,or) + &
+ ( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
+ a(3,1) * var(:,:,or-2) + a(4,1) * var(:,:,or-3) ) * &
+ idel / dz(:,:,or)
+!$OMP END PARALLEL WORKSHARE
+
+ kr = or - 4
+ end if
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+ ( -20.0_wp * var(:,:,kl:kr) + &
+ 15.0_wp * ( var(:,:,kl-1:kr-1) + &
+ var(:,:,kl+1:kr+1) ) - &
+ 6.0_wp * ( var(:,:,kl-2:kr-2) + &
+ var(:,:,kl+2:kr+2) ) + &
+ ( var(:,:,kl-3:kr-3) + &
+ var(:,:,kl+3:kr+3) ) ) * &
+ idel / dz(:,:,kl:kr)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+contains
+
+ subroutine set_coeff ( a )
+
+ implicit none
+
+ CCTK_REAL, dimension(7,4), intent(OUT) :: a
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+
+ a(1,1) = -2.8235294117647058823529411764705882352941176470588_wp
+ a(2,1) = 8.4705882352941176470588235294117647058823529411765_wp
+ a(3,1) = -8.4705882352941176470588235294117647058823529411765_wp
+ a(4,1) = 2.8235294117647058823529411764705882352941176470588_wp
+ a(5,1) = 0.0_wp
+ a(6,1) = 0.0_wp
+ a(7,1) = 0.0_wp
+ a(1,2) = 2.4406779661016949152542372881355932203389830508475_wp
+ a(2,2) = -8.1355932203389830508474576271186440677966101694915_wp
+ a(3,2) = 9.7627118644067796610169491525423728813559322033898_wp
+ a(4,2) = -4.8813559322033898305084745762711864406779661016949_wp
+ a(5,2) = 0.81355932203389830508474576271186440677966101694915_wp
+ a(6,2) = 0.0_wp
+ a(7,2) = 0.0_wp
+ a(1,3) = -3.3488372093023255813953488372093023255813953488372_wp
+ a(2,3) = 13.395348837209302325581395348837209302325581395349_wp
+ a(3,3) = -21.209302325581395348837209302325581395348837209302_wp
+ a(4,3) = 16.744186046511627906976744186046511627906976744186_wp
+ a(5,3) = -6.6976744186046511627906976744186046511627906976744_wp
+ a(6,3) = 1.1162790697674418604651162790697674418604651162791_wp
+ a(7,3) = 0.0_wp
+ a(1,4) = 0.97959183673469387755102040816326530612244897959184_wp
+ a(2,4) = -5.8775510204081632653061224489795918367346938775510_wp
+ a(3,4) = 14.693877551020408163265306122448979591836734693878_wp
+ a(4,4) = -19.591836734693877551020408163265306122448979591837_wp
+ a(5,4) = 14.693877551020408163265306122448979591836734693878_wp
+ a(6,4) = -5.8775510204081632653061224489795918367346938775510_wp
+ a(7,4) = 0.97959183673469387755102040816326530612244897959184_wp
+
+ end subroutine set_coeff
+
+end subroutine dissipation_4_2_delta
diff --git a/src/Dissipation_6_3_alt.F90 b/src/Dissipation_6_3_alt.F90
index 1a4e00b..310d245 100644
--- a/src/Dissipation_6_3_alt.F90
+++ b/src/Dissipation_6_3_alt.F90
@@ -33,6 +33,8 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
il = 1 + gsize(1)
else
ol = offset(1)
+
+!$OMP PARALLEL WORKSHARE
rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
a(3,1) * var(3+ol,:,:) + a(4,1) * var(4+ol,:,:) + &
@@ -63,6 +65,7 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(6,6) * var(6+ol,:,:) + a(7,6) * var(7+ol,:,:) + &
a(8,6) * var(8+ol,:,:) + a(9,6) * var(9+ol,:,:) + &
a(10,6) * var(10+ol,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
il = 7 + ol
end if
@@ -70,6 +73,8 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
ir = ni - gsize(1)
else
or = ni - offset(2)
+
+!$OMP PARALLEL WORKSHARE
rhs(or-5,:,:) = rhs(or-5,:,:) + &
( a(2,6) * var(or-1,:,:) + a(3,6) * var(or-2,:,:) + &
a(4,6) * var(or-3,:,:) + a(5,6) * var(or-4,:,:) + &
@@ -100,9 +105,12 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
a(3,1) * var(or-2,:,:) + a(4,1) * var(or-3,:,:) + &
a(5,1) * var(or-4,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
ir = or - 6
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
( -70.0_wp * var(il:ir,:,:) + &
56.0_wp * ( var(il-1:ir-1,:,:) + &
@@ -113,6 +121,7 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(il+3:ir+3,:,:) ) - &
( var(il-4:ir-4,:,:) + &
var(il+4:ir+4,:,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
if ( zero_derivs_y == 0 ) then
call set_coeff ( a )
@@ -123,6 +132,8 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jl = 1 + gsize(2)
else
ol = offset(3)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
a(3,1) * var(:,3+ol,:) + a(4,1) * var(:,4+ol,:) + &
@@ -153,6 +164,7 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(6,6) * var(:,6+ol,:) + a(7,6) * var(:,7+ol,:) + &
a(8,6) * var(:,8+ol,:) + a(9,6) * var(:,9+ol,:) + &
a(10,6) * var(:,10+ol,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jl = 7 + ol
end if
@@ -160,6 +172,8 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jr = nj - gsize(2)
else
or = nj - offset(4)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,or-5,:) = rhs(:,or-5,:) + &
( a(2,6) * var(:,or-1,:) + a(3,6) * var(:,or-2,:) + &
a(4,6) * var(:,or-3,:) + a(5,6) * var(:,or-4,:) + &
@@ -190,9 +204,12 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
a(3,1) * var(:,or-2,:) + a(4,1) * var(:,or-3,:) + &
a(5,1) * var(:,or-4,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jr = or - 6
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
( -70.0_wp * var(:,jl:jr,:) + &
56.0_wp * ( var(:,jl-1:jr-1,:) + &
@@ -203,6 +220,7 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(:,jl+3:jr+3,:) ) - &
( var(:,jl-4:jr-4,:) + &
var(:,jl+4:jr+4,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
if ( zero_derivs_z == 0 ) then
@@ -214,6 +232,8 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kl = 1 + gsize(3)
else
ol = offset(5)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
a(3,1) * var(:,:,3+ol) + a(4,1) * var(:,:,4+ol) + &
@@ -244,6 +264,7 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(6,6) * var(:,:,6+ol) + a(7,6) * var(:,:,7+ol) + &
a(8,6) * var(:,:,8+ol) + a(9,6) * var(:,:,9+ol) + &
a(10,6) * var(:,:,10+ol) ) * idel
+!$OMP END PARALLEL WORKSHARE
kl = 7 + ol
end if
@@ -251,6 +272,8 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kr = nk - gsize(3)
else
or = nk - offset(6)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,or-5) = rhs(:,:,or-5) + &
( a(2,6) * var(:,:,or-1) + a(3,6) * var(:,:,or-2) + &
a(4,6) * var(:,:,or-3) + a(5,6) * var(:,:,or-4) + &
@@ -281,8 +304,12 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
a(3,1) * var(:,:,or-2) + a(4,1) * var(:,:,or-3) + &
a(5,1) * var(:,:,or-4) ) * idel
+!$OMP END PARALLEL WORKSHARE
+
kr = or - 6
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
( -70.0_wp * var(:,:,kl:kr) + &
56.0_wp * ( var(:,:,kl-1:kr-1) + &
@@ -293,6 +320,7 @@ subroutine dissipation_6_3_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(:,:,kl+3:kr+3) ) - &
( var(:,:,kl-4:kr-4) + &
var(:,:,kl+4:kr+4) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
contains
@@ -369,3 +397,413 @@ contains
end subroutine set_coeff
end subroutine dissipation_6_3_alt
+
+subroutine dissipation_6_3_delta (var, ni, nj, nk, bb, gsize, offset, &
+ dx, dy, dz, epsilon, rhs)
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: ni, nj, nk
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
+ CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
+ CCTK_INT, dimension(6), intent(in) :: bb
+ CCTK_INT, dimension(3), intent(in) :: gsize
+ CCTK_INT, dimension(6), intent(in) :: offset
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: dx, dy, dz
+ CCTK_REAL, intent(in) :: epsilon
+
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+ CCTK_REAL, dimension(10,6) :: a
+ CCTK_REAL :: idel
+
+ CCTK_INT :: il, ir, jl, jr, kl, kr, ol, or
+
+ call set_coeff ( a )
+
+ idel = epsilon / 256
+
+ if ( bb(1) == 0 ) then
+ il = 1 + gsize(1)
+ else
+ ol = offset(1)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
+ ( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
+ a(3,1) * var(3+ol,:,:) + a(4,1) * var(4+ol,:,:) + &
+ a(5,1) * var(5+ol,:,:) ) * idel / dx(1+ol,:,:)
+ rhs(2+ol,:,:) = rhs(2+ol,:,:) + &
+ ( a(1,2) * var(1+ol,:,:) + a(2,2) * var(2+ol,:,:) + &
+ a(3,2) * var(3+ol,:,:) + a(4,2) * var(4+ol,:,:) + &
+ a(5,2) * var(5+ol,:,:) + a(6,2) * var(6+ol,:,:) ) * &
+ idel / dx(2+ol,:,:)
+ rhs(3+ol,:,:) = rhs(3+ol,:,:) + &
+ ( a(1,3) * var(1+ol,:,:) + a(2,3) * var(2+ol,:,:) + &
+ a(3,3) * var(3+ol,:,:) + a(4,3) * var(4+ol,:,:) + &
+ a(5,3) * var(5+ol,:,:) + a(6,3) * var(6+ol,:,:) + &
+ a(7,3) * var(7+ol,:,:) ) * idel / dx(3+ol,:,:)
+ rhs(4+ol,:,:) = rhs(4+ol,:,:) + &
+ ( a(1,4) * var(1+ol,:,:) + a(2,4) * var(2+ol,:,:) + &
+ a(3,4) * var(3+ol,:,:) + a(4,4) * var(4+ol,:,:) + &
+ a(5,4) * var(5+ol,:,:) + a(6,4) * var(6+ol,:,:) + &
+ a(7,4) * var(7+ol,:,:) + a(8,4) * var(8+ol,:,:) ) * &
+ idel / dx(4+ol,:,:)
+ rhs(5+ol,:,:) = rhs(5+ol,:,:) + &
+ ( a(1,5) * var(1+ol,:,:) + a(2,5) * var(2+ol,:,:) + &
+ a(3,5) * var(3+ol,:,:) + a(4,5) * var(4+ol,:,:) + &
+ a(5,5) * var(5+ol,:,:) + a(6,5) * var(6+ol,:,:) + &
+ a(7,5) * var(7+ol,:,:) + a(8,5) * var(8+ol,:,:) + &
+ a(9,5) * var(9+ol,:,:) ) * idel / dx(5+ol,:,:)
+ rhs(6+ol,:,:) = rhs(6+ol,:,:) + &
+ ( a(2,6) * var(2+ol,:,:) + a(3,6) * var(3+ol,:,:) + &
+ a(4,6) * var(4+ol,:,:) + a(5,6) * var(5+ol,:,:) + &
+ a(6,6) * var(6+ol,:,:) + a(7,6) * var(7+ol,:,:) + &
+ a(8,6) * var(8+ol,:,:) + a(9,6) * var(9+ol,:,:) + &
+ a(10,6) * var(10+ol,:,:) ) * idel / dx(6+ol,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ il = 7 + ol
+ end if
+ if ( bb(2) == 0 ) then
+ ir = ni - gsize(1)
+ else
+ or = ni - offset(2)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(or-5,:,:) = rhs(or-5,:,:) + &
+ ( a(2,6) * var(or-1,:,:) + a(3,6) * var(or-2,:,:) + &
+ a(4,6) * var(or-3,:,:) + a(5,6) * var(or-4,:,:) + &
+ a(6,6) * var(or-5,:,:) + a(7,6) * var(or-6,:,:) + &
+ a(8,6) * var(or-7,:,:) + a(9,6) * var(or-8,:,:) + &
+ a(10,6) * var(or-9,:,:) ) * idel / dx(or-5,:,:)
+ rhs(or-4,:,:) = rhs(or-4,:,:) + &
+ ( a(1,5) * var(or,:,:) + a(2,5) * var(or-1,:,:) + &
+ a(3,5) * var(or-2,:,:) + a(4,5) * var(or-3,:,:) + &
+ a(5,5) * var(or-4,:,:) + a(6,5) * var(or-5,:,:) + &
+ a(7,5) * var(or-6,:,:) + a(8,5) * var(or-7,:,:) + &
+ a(9,5) * var(or-8,:,:) ) * idel / dx(or-4,:,:)
+ rhs(or-3,:,:) = rhs(or-3,:,:) + &
+ ( a(1,4) * var(or,:,:) + a(2,4) * var(or-1,:,:) + &
+ a(3,4) * var(or-2,:,:) + a(4,4) * var(or-3,:,:) + &
+ a(5,4) * var(or-4,:,:) + a(6,4) * var(or-5,:,:) + &
+ a(7,4) * var(or-6,:,:) + a(8,4) * var(or-7,:,:) ) * &
+ idel / dx(or-3,:,:)
+ rhs(or-2,:,:) = rhs(or-2,:,:) + &
+ ( a(1,3) * var(or,:,:) + a(2,3) * var(or-1,:,:) + &
+ a(3,3) * var(or-2,:,:) + a(4,3) * var(or-3,:,:) + &
+ a(5,3) * var(or-4,:,:) + a(6,3) * var(or-5,:,:) + &
+ a(7,3) * var(or-6,:,:) ) * idel / dx(or-2,:,:)
+ rhs(or-1,:,:) = rhs(or-1,:,:) + &
+ ( a(1,2) * var(or,:,:) + a(2,2) * var(or-1,:,:) + &
+ a(3,2) * var(or-2,:,:) + a(4,2) * var(or-3,:,:) + &
+ a(5,2) * var(or-4,:,:) + a(6,2) * var(or-5,:,:) ) * &
+ idel / dx(or-1,:,:)
+ rhs(or,:,:) = rhs(or,:,:) + &
+ ( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
+ a(3,1) * var(or-2,:,:) + a(4,1) * var(or-3,:,:) + &
+ a(5,1) * var(or-4,:,:) ) * idel / dx(or,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ ir = or - 6
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+ ( -70.0_wp * var(il:ir,:,:) + &
+ 56.0_wp * ( var(il-1:ir-1,:,:) + &
+ var(il+1:ir+1,:,:) ) - &
+ 28.0_wp * ( var(il-2:ir-2,:,:) + &
+ var(il+2:ir+2,:,:) ) + &
+ 8.0_wp * ( var(il-3:ir-3,:,:) + &
+ var(il+3:ir+3,:,:) ) - &
+ ( var(il-4:ir-4,:,:) + &
+ var(il+4:ir+4,:,:) ) ) * idel / dx(il:ir,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ if ( zero_derivs_y == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 256
+
+ if ( bb(3) == 0 ) then
+ jl = 1 + gsize(2)
+ else
+ ol = offset(3)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
+ ( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
+ a(3,1) * var(:,3+ol,:) + a(4,1) * var(:,4+ol,:) + &
+ a(5,1) * var(:,5+ol,:) ) * idel / dy(:,1+ol,:)
+ rhs(:,2+ol,:) = rhs(:,2+ol,:) + &
+ ( a(1,2) * var(:,1+ol,:) + a(2,2) * var(:,2+ol,:) + &
+ a(3,2) * var(:,3+ol,:) + a(4,2) * var(:,4+ol,:) + &
+ a(5,2) * var(:,5+ol,:) + a(6,2) * var(:,6+ol,:) ) * &
+ idel / dy(:,2+ol,:)
+ rhs(:,3+ol,:) = rhs(:,3+ol,:) + &
+ ( a(1,3) * var(:,1+ol,:) + a(2,3) * var(:,2+ol,:) + &
+ a(3,3) * var(:,3+ol,:) + a(4,3) * var(:,4+ol,:) + &
+ a(5,3) * var(:,5+ol,:) + a(6,3) * var(:,6+ol,:) + &
+ a(7,3) * var(:,7+ol,:) ) * idel / dy(:,3+ol,:)
+ rhs(:,4+ol,:) = rhs(:,4+ol,:) + &
+ ( a(1,4) * var(:,1+ol,:) + a(2,4) * var(:,2+ol,:) + &
+ a(3,4) * var(:,3+ol,:) + a(4,4) * var(:,4+ol,:) + &
+ a(5,4) * var(:,5+ol,:) + a(6,4) * var(:,6+ol,:) + &
+ a(7,4) * var(:,7+ol,:) + a(8,4) * var(:,8+ol,:) ) * &
+ idel / dy(:,4+ol,:)
+ rhs(:,5+ol,:) = rhs(:,5+ol,:) + &
+ ( a(1,5) * var(:,1+ol,:) + a(2,5) * var(:,2+ol,:) + &
+ a(3,5) * var(:,3+ol,:) + a(4,5) * var(:,4+ol,:) + &
+ a(5,5) * var(:,5+ol,:) + a(6,5) * var(:,6+ol,:) + &
+ a(7,5) * var(:,7+ol,:) + a(8,5) * var(:,8+ol,:) + &
+ a(9,5) * var(:,9+ol,:) ) * idel / dy(:,5+ol,:)
+ rhs(:,6+ol,:) = rhs(:,6+ol,:) + &
+ ( a(2,6) * var(:,2+ol,:) + a(3,6) * var(:,3+ol,:) + &
+ a(4,6) * var(:,4+ol,:) + a(5,6) * var(:,5+ol,:) + &
+ a(6,6) * var(:,6+ol,:) + a(7,6) * var(:,7+ol,:) + &
+ a(8,6) * var(:,8+ol,:) + a(9,6) * var(:,9+ol,:) + &
+ a(10,6) * var(:,10+ol,:) ) * idel / dy(:,6+ol,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jl = 7 + ol
+ end if
+ if ( bb(4) == 0 ) then
+ jr = nj - gsize(2)
+ else
+ or = nj - offset(4)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,or-5,:) = rhs(:,or-5,:) + &
+ ( a(2,6) * var(:,or-1,:) + a(3,6) * var(:,or-2,:) + &
+ a(4,6) * var(:,or-3,:) + a(5,6) * var(:,or-4,:) + &
+ a(6,6) * var(:,or-5,:) + a(7,6) * var(:,or-6,:) + &
+ a(8,6) * var(:,or-7,:) + a(9,6) * var(:,or-8,:) + &
+ a(10,6) * var(:,or-9,:) ) * idel / dy(:,or-5,:)
+ rhs(:,or-4,:) = rhs(:,or-4,:) + &
+ ( a(1,5) * var(:,or,:) + a(2,5) * var(:,or-1,:) + &
+ a(3,5) * var(:,or-2,:) + a(4,5) * var(:,or-3,:) + &
+ a(5,5) * var(:,or-4,:) + a(6,5) * var(:,or-5,:) + &
+ a(7,5) * var(:,or-6,:) + a(8,5) * var(:,or-7,:) + &
+ a(9,5) * var(:,or-8,:) ) * idel / dy(:,or-4,:)
+ rhs(:,or-3,:) = rhs(:,or-3,:) + &
+ ( a(1,4) * var(:,or,:) + a(2,4) * var(:,or-1,:) + &
+ a(3,4) * var(:,or-2,:) + a(4,4) * var(:,or-3,:) + &
+ a(5,4) * var(:,or-4,:) + a(6,4) * var(:,or-5,:) + &
+ a(7,4) * var(:,or-6,:) + a(8,4) * var(:,or-7,:) ) * &
+ idel / dy(:,or-3,:)
+ rhs(:,or-2,:) = rhs(:,or-2,:) + &
+ ( a(1,3) * var(:,or,:) + a(2,3) * var(:,or-1,:) + &
+ a(3,3) * var(:,or-2,:) + a(4,3) * var(:,or-3,:) + &
+ a(5,3) * var(:,or-4,:) + a(6,3) * var(:,or-5,:) + &
+ a(7,3) * var(:,or-6,:) ) * idel / dy(:,or-2,:)
+ rhs(:,or-1,:) = rhs(:,or-1,:) + &
+ ( a(1,2) * var(:,or,:) + a(2,2) * var(:,or-1,:) + &
+ a(3,2) * var(:,or-2,:) + a(4,2) * var(:,or-3,:) + &
+ a(5,2) * var(:,or-4,:) + a(6,2) * var(:,or-5,:) ) * &
+ idel / dy(:,or-1,:)
+ rhs(:,or,:) = rhs(:,or,:) + &
+ ( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
+ a(3,1) * var(:,or-2,:) + a(4,1) * var(:,or-3,:) + &
+ a(5,1) * var(:,or-4,:) ) * idel / dy(:,or,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jr = or - 6
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+ ( -70.0_wp * var(:,jl:jr,:) + &
+ 56.0_wp * ( var(:,jl-1:jr-1,:) + &
+ var(:,jl+1:jr+1,:) ) - &
+ 28.0_wp * ( var(:,jl-2:jr-2,:) + &
+ var(:,jl+2:jr+2,:) ) + &
+ 8.0_wp * ( var(:,jl-3:jr-3,:) + &
+ var(:,jl+3:jr+3,:) ) - &
+ ( var(:,jl-4:jr-4,:) + &
+ var(:,jl+4:jr+4,:) ) ) * &
+ idel / dy(:,jl:jr,:)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+ if ( zero_derivs_z == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 256
+
+ if ( bb(5) == 0 ) then
+ kl = 1 + gsize(3)
+ else
+ ol = offset(5)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
+ ( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
+ a(3,1) * var(:,:,3+ol) + a(4,1) * var(:,:,4+ol) + &
+ a(5,1) * var(:,:,5+ol) ) * idel / dz(:,:,1+ol)
+ rhs(:,:,2+ol) = rhs(:,:,2+ol) + &
+ ( a(1,2) * var(:,:,1+ol) + a(2,2) * var(:,:,2+ol) + &
+ a(3,2) * var(:,:,3+ol) + a(4,2) * var(:,:,4+ol) + &
+ a(5,2) * var(:,:,5+ol) + a(6,2) * var(:,:,6+ol) ) * &
+ idel / dz(:,:,2+ol)
+ rhs(:,:,3+ol) = rhs(:,:,3+ol) + &
+ ( a(1,3) * var(:,:,1+ol) + a(2,3) * var(:,:,2+ol) + &
+ a(3,3) * var(:,:,3+ol) + a(4,3) * var(:,:,4+ol) + &
+ a(5,3) * var(:,:,5+ol) + a(6,3) * var(:,:,6+ol) + &
+ a(7,3) * var(:,:,7+ol) ) * idel / dz(:,:,3+ol)
+ rhs(:,:,4+ol) = rhs(:,:,4+ol) + &
+ ( a(1,4) * var(:,:,1+ol) + a(2,4) * var(:,:,2+ol) + &
+ a(3,4) * var(:,:,3+ol) + a(4,4) * var(:,:,4+ol) + &
+ a(5,4) * var(:,:,5+ol) + a(6,4) * var(:,:,6+ol) + &
+ a(7,4) * var(:,:,7+ol) + a(8,4) * var(:,:,8+ol) ) * &
+ idel / dz(:,:,4+ol)
+ rhs(:,:,5+ol) = rhs(:,:,5+ol) + &
+ ( a(1,5) * var(:,:,1+ol) + a(2,5) * var(:,:,2+ol) + &
+ a(3,5) * var(:,:,3+ol) + a(4,5) * var(:,:,4+ol) + &
+ a(5,5) * var(:,:,5+ol) + a(6,5) * var(:,:,6+ol) + &
+ a(7,5) * var(:,:,7+ol) + a(8,5) * var(:,:,8+ol) + &
+ a(9,5) * var(:,:,9+ol) ) * idel / dz(:,:,5+ol)
+ rhs(:,:,6+ol) = rhs(:,:,6+ol) + &
+ ( a(2,6) * var(:,:,2+ol) + a(3,6) * var(:,:,3+ol) + &
+ a(4,6) * var(:,:,4+ol) + a(5,6) * var(:,:,5+ol) + &
+ a(6,6) * var(:,:,6+ol) + a(7,6) * var(:,:,7+ol) + &
+ a(8,6) * var(:,:,8+ol) + a(9,6) * var(:,:,9+ol) + &
+ a(10,6) * var(:,:,10+ol) ) * idel / dz(:,:,6+ol)
+!$OMP END PARALLEL WORKSHARE
+
+ kl = 7 + ol
+ end if
+ if ( bb(6) == 0 ) then
+ kr = nk - gsize(3)
+ else
+ or = nk - offset(6)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,or-5) = rhs(:,:,or-5) + &
+ ( a(2,6) * var(:,:,or-1) + a(3,6) * var(:,:,or-2) + &
+ a(4,6) * var(:,:,or-3) + a(5,6) * var(:,:,or-4) + &
+ a(6,6) * var(:,:,or-5) + a(7,6) * var(:,:,or-6) + &
+ a(8,6) * var(:,:,or-7) + a(9,6) * var(:,:,or-8) + &
+ a(10,6) * var(:,:,or-9) ) * idel / dz(:,:,or-5)
+ rhs(:,:,or-4) = rhs(:,:,or-4) + &
+ ( a(1,5) * var(:,:,or) + a(2,5) * var(:,:,or-1) + &
+ a(3,5) * var(:,:,or-2) + a(4,5) * var(:,:,or-3) + &
+ a(5,5) * var(:,:,or-4) + a(6,5) * var(:,:,or-5) + &
+ a(7,5) * var(:,:,or-6) + a(8,5) * var(:,:,or-7) + &
+ a(9,5) * var(:,:,or-8) ) * idel / dz(:,:,or-4)
+ rhs(:,:,or-3) = rhs(:,:,or-3) + &
+ ( a(1,4) * var(:,:,or) + a(2,4) * var(:,:,or-1) + &
+ a(3,4) * var(:,:,or-2) + a(4,4) * var(:,:,or-3) + &
+ a(5,4) * var(:,:,or-4) + a(6,4) * var(:,:,or-5) + &
+ a(7,4) * var(:,:,or-6) + a(8,4) * var(:,:,or-7) ) * &
+ idel / dz(:,:,or-3)
+ rhs(:,:,or-2) = rhs(:,:,or-2) + &
+ ( a(1,3) * var(:,:,or) + a(2,3) * var(:,:,or-1) + &
+ a(3,3) * var(:,:,or-2) + a(4,3) * var(:,:,or-3) + &
+ a(5,3) * var(:,:,or-4) + a(6,3) * var(:,:,or-5) + &
+ a(7,3) * var(:,:,or-6) ) * idel / dz(:,:,or-2)
+ rhs(:,:,or-1) = rhs(:,:,or-1) + &
+ ( a(1,2) * var(:,:,or) + a(2,2) * var(:,:,or-1) + &
+ a(3,2) * var(:,:,or-2) + a(4,2) * var(:,:,or-3) + &
+ a(5,2) * var(:,:,or-4) + a(6,2) * var(:,:,or-5) ) * &
+ idel / dz(:,:,or-1)
+ rhs(:,:,or) = rhs(:,:,or) + &
+ ( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
+ a(3,1) * var(:,:,or-2) + a(4,1) * var(:,:,or-3) + &
+ a(5,1) * var(:,:,or-4) ) * idel / dz(:,:,or)
+!$OMP END PARALLEL WORKSHARE
+
+ kr = or - 6
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+ ( -70.0_wp * var(:,:,kl:kr) + &
+ 56.0_wp * ( var(:,:,kl-1:kr-1) + &
+ var(:,:,kl+1:kr+1) ) - &
+ 28.0_wp * ( var(:,:,kl-2:kr-2) + &
+ var(:,:,kl+2:kr+2) ) + &
+ 8.0_wp * ( var(:,:,kl-3:kr-3) + &
+ var(:,:,kl+3:kr+3) ) - &
+ ( var(:,:,kl-4:kr-4) + &
+ var(:,:,kl+4:kr+4) ) ) * &
+ idel / dz(:,:,kl:kr)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+contains
+
+ subroutine set_coeff ( a )
+
+ implicit none
+
+ CCTK_REAL, dimension(10,6), intent(OUT) :: a
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+
+ a(1,1) = -3.1650670378782328375705179866656897941241116565316_wp
+ a(2,1) = 12.660268151512931350282071946662759176496446626126_wp
+ a(3,1) = -18.990402227269397025423107919994138764744669939190_wp
+ a(4,1) = 12.660268151512931350282071946662759176496446626126_wp
+ a(5,1) = -3.1650670378782328375705179866656897941241116565316_wp
+ a(6,1) = 0.0_wp
+ a(7,1) = 0.0_wp
+ a(8,1) = 0.0_wp
+ a(9,1) = 0.0_wp
+ a(10,1) = 0.0_wp
+ a(1,2) = 2.8768833763422958461666527928077915591442603845834_wp
+ a(2,2) = -12.226754349454757346208274369433114126363106634479_wp
+ a(3,2) = 20.138183634396070923166569549654540914009822692084_wp
+ a(4,2) = -15.822858569882627153916590360442853575293432115209_wp
+ a(5,2) = 5.7537667526845916923333055856155831182885207691667_wp
+ a(6,2) = -0.71922084408557396154166319820194788978606509614584_wp
+ a(7,2) = 0.0_wp
+ a(8,2) = 0.0_wp
+ a(9,2) = 0.0_wp
+ a(10,2) = 0.0_wp
+ a(1,3) = -9.5610475839173736628550350424197713021025451862781_wp
+ a(2,3) = 44.618222058281077093323496864625599409811877535965_wp
+ a(3,3) = -84.455920324603467355219476208041313168572482478790_wp
+ a(4,3) = 82.862412393950571744743637034304684618222058281077_wp
+ a(5,3) = -44.618222058281077093323496864625599409811877535965_wp
+ a(6,3) = 12.748063445223164883806713389893028402803393581704_wp
+ a(7,3) = -1.5935079306528956104758391737366285503504241977130_wp
+ a(8,3) = 0.0_wp
+ a(9,3) = 0.0_wp
+ a(10,3) = 0.0_wp
+ a(1,4) = 3.2244821795111028176898675125956335137152453816010_wp
+ a(2,4) = -17.734651987311065497294271319275984325433849598806_wp
+ a(3,4) = 41.918268333644336629968277663743235678298189960814_wp
+ a(4,4) = -55.622317596566523605150214592274678111587982832618_wp
+ a(5,4) = 45.142750513155439447658145176338869192013435342415_wp
+ a(6,4) = -22.571375256577719723829072588169434596006717671207_wp
+ a(7,4) = 6.4489643590222056353797350251912670274304907632021_wp
+ a(8,4) = -0.80612054487777570442246687814890837842881134540026_wp
+ a(9,4) = 0.0_wp
+ a(10,4) = 0.0_wp
+ a(1,5) = -1.0968642884346832550463374381109559476958232829758_wp
+ a(2,5) = 8.7749143074774660403706995048876475815665862638060_wp
+ a(3,5) = -30.712200076171131141297448267106766535483051923321_wp
+ a(4,5) = 61.424400152342262282594896534213533070966103846642_wp
+ a(5,5) = -76.780500190427827853243620667766916338707629808303_wp
+ a(6,5) = 61.424400152342262282594896534213533070966103846642_wp
+ a(7,5) = -30.712200076171131141297448267106766535483051923321_wp
+ a(8,5) = 8.7749143074774660403706995048876475815665862638060_wp
+ a(9,5) = -1.0968642884346832550463374381109559476958232829758_wp
+ a(10,5) = 0.0_wp
+ a(1,6) = 0.0_wp
+ a(2,6) = -0.98627885208100271683294901942878016483641925983425_wp
+ a(3,6) = 7.8902308166480217346635921554302413186913540786740_wp
+ a(4,6) = -27.615807858268076071322572544005844615419739275359_wp
+ a(5,6) = 55.231615716536152142645145088011689230839478550718_wp
+ a(6,6) = -69.039519645670190178306431360014611538549348188398_wp
+ a(7,6) = 55.231615716536152142645145088011689230839478550718_wp
+ a(8,6) = -27.615807858268076071322572544005844615419739275359_wp
+ a(9,6) = 7.8902308166480217346635921554302413186913540786740_wp
+ a(10,6) = -0.98627885208100271683294901942878016483641925983425_wp
+
+ end subroutine set_coeff
+
+end subroutine dissipation_6_3_delta
diff --git a/src/Dissipation_8_4_alt.F90 b/src/Dissipation_8_4_alt.F90
index 572b1ad..993600f 100644
--- a/src/Dissipation_8_4_alt.F90
+++ b/src/Dissipation_8_4_alt.F90
@@ -33,6 +33,8 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
il = 1 + gsize(1)
else
ol = offset(1)
+
+!$OMP PARALLEL WORKSHARE
rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
a(3,1) * var(3+ol,:,:) + a(4,1) * var(4+ol,:,:) + &
@@ -80,6 +82,7 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(9,8) * var(9+ol,:,:) + a(10,8) * var(10+ol,:,:) + &
a(11,8) * var(11+ol,:,:) + a(12,8) * var(12+ol,:,:) + &
a(13,8) * var(13+ol,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
il = 9 + ol
end if
@@ -87,6 +90,8 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
ir = ni - gsize(1)
else
or = ni - offset(2)
+
+!$OMP PARALLEL WORKSHARE
rhs(or-7,:,:) = rhs(or-7,:,:) + &
( a(3,8) * var(or-2,:,:) + a(4,8) * var(or-3,:,:) + &
a(5,8) * var(or-4,:,:) + a(6,8) * var(or-5,:,:) + &
@@ -134,9 +139,12 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
a(3,1) * var(or-2,:,:) + a(4,1) * var(or-3,:,:) + &
a(5,1) * var(or-4,:,:) + a(6,1) * var(or-5,:,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
ir = or - 8
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
( -252.0_wp * var(il:ir,:,:) + &
210.0_wp * ( var(il-1:ir-1,:,:) + &
@@ -149,6 +157,7 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(il+4:ir+4,:,:) ) + &
( var(il-5:ir-5,:,:) + &
var(il+5:ir+5,:,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
if ( zero_derivs_y == 0 ) then
call set_coeff ( a )
@@ -159,6 +168,8 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jl = 1 + gsize(2)
else
ol = offset(3);
+
+!$OMP PARALLEL WORKSHARE
rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
a(3,1) * var(:,3+ol,:) + a(4,1) * var(:,4+ol,:) + &
@@ -206,6 +217,7 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(9,8) * var(:,9+ol,:) + a(10,8) * var(:,10+ol,:) + &
a(11,8) * var(:,11+ol,:) + a(12,8) * var(:,12+ol,:) + &
a(13,8) * var(:,13+ol,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jl = 9 + ol
end if
@@ -213,6 +225,8 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
jr = nj - gsize(2)
else
or = nj - offset(4)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,or-7,:) = rhs(:,or-7,:) + &
( a(3,8) * var(:,or-2,:) + a(4,8) * var(:,or-3,:) + &
a(5,8) * var(:,or-4,:) + a(6,8) * var(:,or-5,:) + &
@@ -260,9 +274,12 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
a(3,1) * var(:,or-2,:) + a(4,1) * var(:,or-3,:) + &
a(5,1) * var(:,or-4,:) + a(6,1) * var(:,or-5,:) ) * idel
+!$OMP END PARALLEL WORKSHARE
jr = or - 8
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
( -252.0_wp * var(:,jl:jr,:) + &
210.0_wp * ( var(:,jl-1:jr-1,:) + &
@@ -275,6 +292,7 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(:,jl+4:jr+4,:) ) + &
( var(:,jl-5:jr-5,:) + &
var(:,jl+5:jr+5,:) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
if ( zero_derivs_z == 0 ) then
@@ -286,6 +304,8 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kl = 1 + gsize(3)
else
ol = offset(5);
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
a(3,1) * var(:,:,3+ol) + a(4,1) * var(:,:,4+ol) + &
@@ -333,6 +353,7 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
a(9,8) * var(:,:,9+ol) + a(10,8) * var(:,:,10+ol) + &
a(11,8) * var(:,:,11+ol) + a(12,8) * var(:,:,12+ol) + &
a(13,8) * var(:,:,13+ol ) ) * idel
+!$OMP END PARALLEL WORKSHARE
kl = 9 + ol
end if
@@ -340,6 +361,8 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
kr = nk - gsize(3)
else
or = nk - offset(6)
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,or-7) = rhs(:,:,or-7) + &
( a(3,8) * var(:,:,or-2) + a(4,8) * var(:,:,or-3) + &
a(5,8) * var(:,:,or-4) + a(6,8) * var(:,:,or-5) + &
@@ -387,9 +410,12 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
a(3,1) * var(:,:,or-2) + a(4,1) * var(:,:,or-3) + &
a(5,1) * var(:,:,or-4) + a(6,1) * var(:,:,or-5) ) * idel
+!$OMP END PARALLEL WORKSHARE
kr = or - 8
end if
+
+!$OMP PARALLEL WORKSHARE
rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
( -252.0_wp * var(:,:,kl:kr) + &
210.0_wp * ( var(:,:,kl-1:kr-1) + &
@@ -402,6 +428,7 @@ subroutine dissipation_8_4_alt (var, ni, nj, nk, bb, gsize, offset, delta, epsil
var(:,:,kl+4:kr+4) ) + &
( var(:,:,kl-5:kr-5) + &
var(:,:,kl+5:kr+5) ) ) * idel
+!$OMP END PARALLEL WORKSHARE
end if
contains
@@ -522,3 +549,571 @@ contains
end subroutine set_coeff
end subroutine dissipation_8_4_alt
+
+subroutine dissipation_8_4_delta (var, ni, nj, nk, bb, gsize, offset, &
+ dx, dy, dz, epsilon, rhs)
+
+ implicit none
+
+ DECLARE_CCTK_PARAMETERS
+
+ integer :: ni, nj, nk
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: var
+ CCTK_REAL, dimension(ni,nj,nk), intent(inout) :: rhs
+ CCTK_INT, dimension(6), intent(in) :: bb
+ CCTK_INT, dimension(3), intent(in) :: gsize
+ CCTK_INT, dimension(6), intent(in) :: offset
+ CCTK_REAL, dimension(ni,nj,nk), intent(in) :: dx, dy, dz
+ CCTK_REAL, intent(in) :: epsilon
+
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+ CCTK_REAL, dimension(13,8) :: a
+ CCTK_REAL :: idel
+
+ CCTK_INT :: il, ir, jl, jr, kl, kr, ol, or
+
+ call set_coeff ( a )
+
+ idel = epsilon / 1024
+
+ if ( bb(1) == 0 ) then
+ il = 1 + gsize(1)
+ else
+ ol = offset(1)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(1+ol,:,:) = rhs(1+ol,:,:) + &
+ ( a(1,1) * var(1+ol,:,:) + a(2,1) * var(2+ol,:,:) + &
+ a(3,1) * var(3+ol,:,:) + a(4,1) * var(4+ol,:,:) + &
+ a(5,1) * var(5+ol,:,:) + a(6,1) * var(6+ol,:,:) ) * &
+ idel / dx(1+ol,:,:)
+ rhs(2+ol,:,:) = rhs(2+ol,:,:) + &
+ ( a(1,2) * var(1+ol,:,:) + a(2,2) * var(2+ol,:,:) + &
+ a(3,2) * var(3+ol,:,:) + a(4,2) * var(4+ol,:,:) + &
+ a(5,2) * var(5+ol,:,:) + a(6,2) * var(6+ol,:,:) + &
+ a(7,2) * var(7+ol,:,:) ) * idel / dx(2+ol,:,:)
+ rhs(3+ol,:,:) = rhs(3+ol,:,:) + &
+ ( a(1,3) * var(1+ol,:,:) + a(2,3) * var(2+ol,:,:) + &
+ a(3,3) * var(3+ol,:,:) + a(4,3) * var(4+ol,:,:) + &
+ a(5,3) * var(5+ol,:,:) + a(6,3) * var(6+ol,:,:) + &
+ a(7,3) * var(7+ol,:,:) + a(8,3) * var(8+ol,:,:) ) * &
+ idel / dx(3+ol,:,:)
+ rhs(4+ol,:,:) = rhs(4+ol,:,:) + &
+ ( a(1,4) * var(1+ol,:,:) + a(2,4) * var(2+ol,:,:) + &
+ a(3,4) * var(3+ol,:,:) + a(4,4) * var(4+ol,:,:) + &
+ a(5,4) * var(5+ol,:,:) + a(6,4) * var(6+ol,:,:) + &
+ a(7,4) * var(7+ol,:,:) + a(8,4) * var(8+ol,:,:) + &
+ a(9,4) * var(9+ol,:,:) ) * idel / dx(4+ol,:,:)
+ rhs(5+ol,:,:) = rhs(5+ol,:,:) + &
+ ( a(1,5) * var(1+ol,:,:) + a(2,5) * var(2+ol,:,:) + &
+ a(3,5) * var(3+ol,:,:) + a(4,5) * var(4+ol,:,:) + &
+ a(5,5) * var(5+ol,:,:) + a(6,5) * var(6+ol,:,:) + &
+ a(7,5) * var(7+ol,:,:) + a(8,5) * var(8+ol,:,:) + &
+ a(9,5) * var(9+ol,:,:) + a(10,5) * var(10+ol,:,:) ) * &
+ idel / dx(5+ol,:,:)
+ rhs(6+ol,:,:) = rhs(6+ol,:,:) + &
+ ( a(1,6) * var(1+ol,:,:) + a(2,6) * var(2+ol,:,:) + &
+ a(3,6) * var(3+ol,:,:) + a(4,6) * var(4+ol,:,:) + &
+ a(5,6) * var(5+ol,:,:) + a(6,6) * var(6+ol,:,:) + &
+ a(7,6) * var(7+ol,:,:) + a(8,6) * var(8+ol,:,:) + &
+ a(9,6) * var(9+ol,:,:) + a(10,6) * var(10+ol,:,:) + &
+ a(11,6) * var(11+ol,:,:) ) * idel / dx(6+ol,:,:)
+ rhs(7+ol,:,:) = rhs(7+ol,:,:) + &
+ ( a(2,7) * var(2+ol,:,:) + a(3,7) * var(3+ol,:,:) + &
+ a(4,7) * var(4+ol,:,:) + a(5,7) * var(5+ol,:,:) + &
+ a(6,7) * var(6+ol,:,:) + a(7,7) * var(7+ol,:,:) + &
+ a(8,7) * var(8+ol,:,:) + a(9,7) * var(9+ol,:,:) + &
+ a(10,7) * var(10+ol,:,:) + a(11,7) * var(11+ol,:,:) + &
+ a(12,7) * var(12+ol,:,:) ) * idel / dx(7+ol,:,:)
+ rhs(8+ol,:,:) = rhs(8+ol,:,:) + &
+ ( a(3,8) * var(3+ol,:,:) + a(4,8) * var(4+ol,:,:) + &
+ a(5,8) * var(5+ol,:,:) + a(6,8) * var(6+ol,:,:) + &
+ a(7,8) * var(7+ol,:,:) + a(8,8) * var(8+ol,:,:) + &
+ a(9,8) * var(9+ol,:,:) + a(10,8) * var(10+ol,:,:) + &
+ a(11,8) * var(11+ol,:,:) + a(12,8) * var(12+ol,:,:) + &
+ a(13,8) * var(13+ol,:,:) ) * idel / dx(8+ol,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ il = 9 + ol
+ end if
+ if ( bb(2) == 0 ) then
+ ir = ni - gsize(1)
+ else
+ or = ni - offset(2)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(or-7,:,:) = rhs(or-7,:,:) + &
+ ( a(3,8) * var(or-2,:,:) + a(4,8) * var(or-3,:,:) + &
+ a(5,8) * var(or-4,:,:) + a(6,8) * var(or-5,:,:) + &
+ a(7,8) * var(or-6,:,:) + a(8,8) * var(or-7,:,:) + &
+ a(9,8) * var(or-8,:,:) + a(10,8) * var(or-9,:,:) + &
+ a(11,8) * var(or-10,:,:) + a(12,8) * var(or-11,:,:) + &
+ a(13,8) * var(or-12,:,:) ) * idel / dx(or-7,:,:)
+ rhs(or-6,:,:) = rhs(or-6,:,:) + &
+ ( a(2,7) * var(or-1,:,:) + a(3,7) * var(or-2,:,:) + &
+ a(4,7) * var(or-3,:,:) + a(5,7) * var(or-4,:,:) + &
+ a(6,7) * var(or-5,:,:) + a(7,7) * var(or-6,:,:) + &
+ a(8,7) * var(or-7,:,:) + a(9,7) * var(or-8,:,:) + &
+ a(10,7) * var(or-9,:,:) + a(11,7) * var(or-10,:,:) + &
+ a(12,7) * var(or-11,:,:) ) * idel / dx(or-6,:,:)
+ rhs(or-5,:,:) = rhs(or-5,:,:) + &
+ ( a(1,6) * var(or,:,:) + a(2,6) * var(or-1,:,:) + &
+ a(3,6) * var(or-2,:,:) + a(4,6) * var(or-3,:,:) + &
+ a(5,6) * var(or-4,:,:) + a(6,6) * var(or-5,:,:) + &
+ a(7,6) * var(or-6,:,:) + a(8,6) * var(or-7,:,:) + &
+ a(9,6) * var(or-8,:,:) + a(10,6) * var(or-9,:,:) + &
+ a(11,6) * var(or-10,:,:) ) * idel / dx(or-5,:,:)
+ rhs(or-4,:,:) = rhs(or-4,:,:) + &
+ ( a(1,5) * var(or,:,:) + a(2,5) * var(or-1,:,:) + &
+ a(3,5) * var(or-2,:,:) + a(4,5) * var(or-3,:,:) + &
+ a(5,5) * var(or-4,:,:) + a(6,5) * var(or-5,:,:) + &
+ a(7,5) * var(or-6,:,:) + a(8,5) * var(or-7,:,:) + &
+ a(9,5) * var(or-8,:,:) + a(10,5) * var(or-9,:,:) ) * &
+ idel / dx(or-4,:,:)
+ rhs(or-3,:,:) = rhs(or-3,:,:) + &
+ ( a(1,4) * var(or,:,:) + a(2,4) * var(or-1,:,:) + &
+ a(3,4) * var(or-2,:,:) + a(4,4) * var(or-3,:,:) + &
+ a(5,4) * var(or-4,:,:) + a(6,4) * var(or-5,:,:) + &
+ a(7,4) * var(or-6,:,:) + a(8,4) * var(or-7,:,:) + &
+ a(9,4) * var(or-8,:,:) ) * idel / dx(or-3,:,:)
+ rhs(or-2,:,:) = rhs(or-2,:,:) + &
+ ( a(1,3) * var(or,:,:) + a(2,3) * var(or-1,:,:) + &
+ a(3,3) * var(or-2,:,:) + a(4,3) * var(or-3,:,:) + &
+ a(5,3) * var(or-4,:,:) + a(6,3) * var(or-5,:,:) + &
+ a(7,3) * var(or-6,:,:) + a(8,3) * var(or-7,:,:) ) * &
+ idel / dx(or-2,:,:)
+ rhs(or-1,:,:) = rhs(or-1,:,:) + &
+ ( a(1,2) * var(or,:,:) + a(2,2) * var(or-1,:,:) + &
+ a(3,2) * var(or-2,:,:) + a(4,2) * var(or-3,:,:) + &
+ a(5,2) * var(or-4,:,:) + a(6,2) * var(or-5,:,:) + &
+ a(7,2) * var(or-6,:,:) ) * idel / dx(or-1,:,:)
+ rhs(or,:,:) = rhs(or,:,:) + &
+ ( a(1,1) * var(or,:,:) + a(2,1) * var(or-1,:,:) + &
+ a(3,1) * var(or-2,:,:) + a(4,1) * var(or-3,:,:) + &
+ a(5,1) * var(or-4,:,:) + a(6,1) * var(or-5,:,:) ) * &
+ idel / dx(or,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ ir = or - 8
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(il:ir,:,:) = rhs(il:ir,:,:) + &
+ ( -252.0_wp * var(il:ir,:,:) + &
+ 210.0_wp * ( var(il-1:ir-1,:,:) + &
+ var(il+1:ir+1,:,:) ) - &
+ 120.0_wp * ( var(il-2:ir-2,:,:) + &
+ var(il+2:ir+2,:,:) ) + &
+ 45.0_wp * ( var(il-3:ir-3,:,:) + &
+ var(il+3:ir+3,:,:) ) - &
+ 10.0_wp * ( var(il-4:ir-4,:,:) + &
+ var(il+4:ir+4,:,:) ) + &
+ ( var(il-5:ir-5,:,:) + &
+ var(il+5:ir+5,:,:) ) ) * idel / dx(il:ir,:,:)
+!$OMP END PARALLEL WORKSHARE
+
+ if ( zero_derivs_y == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 1024
+
+ if ( bb(3) == 0 ) then
+ jl = 1 + gsize(2)
+ else
+ ol = offset(3);
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,1+ol,:) = rhs(:,1+ol,:) + &
+ ( a(1,1) * var(:,1+ol,:) + a(2,1) * var(:,2+ol,:) + &
+ a(3,1) * var(:,3+ol,:) + a(4,1) * var(:,4+ol,:) + &
+ a(5,1) * var(:,5+ol,:) + a(6,1) * var(:,6+ol,:) ) * &
+ idel / dy(:,1+ol,:)
+ rhs(:,2+ol,:) = rhs(:,2+ol,:) + &
+ ( a(1,2) * var(:,1+ol,:) + a(2,2) * var(:,2+ol,:) + &
+ a(3,2) * var(:,3+ol,:) + a(4,2) * var(:,4+ol,:) + &
+ a(5,2) * var(:,5+ol,:) + a(6,2) * var(:,6+ol,:) + &
+ a(7,2) * var(:,7+ol,:) ) * idel / dy(:,2+ol,:)
+ rhs(:,3+ol,:) = rhs(:,3+ol,:) + &
+ ( a(1,3) * var(:,1+ol,:) + a(2,3) * var(:,2+ol,:) + &
+ a(3,3) * var(:,3+ol,:) + a(4,3) * var(:,4+ol,:) + &
+ a(5,3) * var(:,5+ol,:) + a(6,3) * var(:,6+ol,:) + &
+ a(7,3) * var(:,7+ol,:) + a(8,3) * var(:,8+ol,:) ) * &
+ idel / dy(:,3+ol,:)
+ rhs(:,4+ol,:) = rhs(:,4+ol,:) + &
+ ( a(1,4) * var(:,1+ol,:) + a(2,4) * var(:,2+ol,:) + &
+ a(3,4) * var(:,3+ol,:) + a(4,4) * var(:,4+ol,:) + &
+ a(5,4) * var(:,5+ol,:) + a(6,4) * var(:,6+ol,:) + &
+ a(7,4) * var(:,7+ol,:) + a(8,4) * var(:,8+ol,:) + &
+ a(9,4) * var(:,9+ol,:) ) * idel / dy(:,4+ol,:)
+ rhs(:,5+ol,:) = rhs(:,5+ol,:) + &
+ ( a(1,5) * var(:,1+ol,:) + a(2,5) * var(:,2+ol,:) + &
+ a(3,5) * var(:,3+ol,:) + a(4,5) * var(:,4+ol,:) + &
+ a(5,5) * var(:,5+ol,:) + a(6,5) * var(:,6+ol,:) + &
+ a(7,5) * var(:,7+ol,:) + a(8,5) * var(:,8+ol,:) + &
+ a(9,5) * var(:,9+ol,:) + a(10,5) * var(:,10+ol,:) ) * &
+ idel / dy(:,5+ol,:)
+ rhs(:,6+ol,:) = rhs(:,6+ol,:) + &
+ ( a(1,6) * var(:,1+ol,:) + a(2,6) * var(:,2+ol,:) + &
+ a(3,6) * var(:,3+ol,:) + a(4,6) * var(:,4+ol,:) + &
+ a(5,6) * var(:,5+ol,:) + a(6,6) * var(:,6+ol,:) + &
+ a(7,6) * var(:,7+ol,:) + a(8,6) * var(:,8+ol,:) + &
+ a(9,6) * var(:,9+ol,:) + a(10,6) * var(:,10+ol,:) + &
+ a(11,6) * var(:,11+ol,:) ) * idel / dy(:,6+ol,:)
+ rhs(:,7+ol,:) = rhs(:,7+ol,:) + &
+ ( a(2,7) * var(:,2+ol,:) + a(3,7) * var(:,3+ol,:) + &
+ a(4,7) * var(:,4+ol,:) + a(5,7) * var(:,5+ol,:) + &
+ a(6,7) * var(:,6+ol,:) + a(7,7) * var(:,7+ol,:) + &
+ a(8,7) * var(:,8+ol,:) + a(9,7) * var(:,9+ol,:) + &
+ a(10,7) * var(:,10+ol,:) + a(11,7) * var(:,11+ol,:) + &
+ a(12,7) * var(:,12+ol,:) ) * idel / dy(:,7+ol,:)
+ rhs(:,8+ol,:) = rhs(:,8+ol,:) + &
+ ( a(3,8) * var(:,3+ol,:) + a(4,8) * var(:,4+ol,:) + &
+ a(5,8) * var(:,5+ol,:) + a(6,8) * var(:,6+ol,:) + &
+ a(7,8) * var(:,7+ol,:) + a(8,8) * var(:,8+ol,:) + &
+ a(9,8) * var(:,9+ol,:) + a(10,8) * var(:,10+ol,:) + &
+ a(11,8) * var(:,11+ol,:) + a(12,8) * var(:,12+ol,:) + &
+ a(13,8) * var(:,13+ol,:) ) * idel / dy(:,8+ol,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jl = 9 + ol
+ end if
+ if ( bb(4) == 0 ) then
+ jr = nj - gsize(2)
+ else
+ or = nj - offset(4)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,or-7,:) = rhs(:,or-7,:) + &
+ ( a(3,8) * var(:,or-2,:) + a(4,8) * var(:,or-3,:) + &
+ a(5,8) * var(:,or-4,:) + a(6,8) * var(:,or-5,:) + &
+ a(7,8) * var(:,or-6,:) + a(8,8) * var(:,or-7,:) + &
+ a(9,8) * var(:,or-8,:) + a(10,8) * var(:,or-9,:) + &
+ a(11,8) * var(:,or-10,:) + a(12,8) * var(:,or-11,:) + &
+ a(13,8) * var(:,or-12,:) ) * idel / dy(:,or-7,:)
+ rhs(:,or-6,:) = rhs(:,or-6,:) + &
+ ( a(2,7) * var(:,or-1,:) + a(3,7) * var(:,or-2,:) + &
+ a(4,7) * var(:,or-3,:) + a(5,7) * var(:,or-4,:) + &
+ a(6,7) * var(:,or-5,:) + a(7,7) * var(:,or-6,:) + &
+ a(8,7) * var(:,or-7,:) + a(9,7) * var(:,or-8,:) + &
+ a(10,7) * var(:,or-9,:) + a(11,7) * var(:,or-10,:) + &
+ a(12,7) * var(:,or-11,:) ) * idel / dy(:,or-6,:)
+ rhs(:,or-5,:) = rhs(:,or-5,:) + &
+ ( a(1,6) * var(:,or,:) + a(2,6) * var(:,or-1,:) + &
+ a(3,6) * var(:,or-2,:) + a(4,6) * var(:,or-3,:) + &
+ a(5,6) * var(:,or-4,:) + a(6,6) * var(:,or-5,:) + &
+ a(7,6) * var(:,or-6,:) + a(8,6) * var(:,or-7,:) + &
+ a(9,6) * var(:,or-8,:) + a(10,6) * var(:,or-9,:) + &
+ a(11,6) * var(:,or-10,:) ) * idel / dy(:,or-5,:)
+ rhs(:,or-4,:) = rhs(:,or-4,:) + &
+ ( a(1,5) * var(:,or,:) + a(2,5) * var(:,or-1,:) + &
+ a(3,5) * var(:,or-2,:) + a(4,5) * var(:,or-3,:) + &
+ a(5,5) * var(:,or-4,:) + a(6,5) * var(:,or-5,:) + &
+ a(7,5) * var(:,or-6,:) + a(8,5) * var(:,or-7,:) + &
+ a(9,5) * var(:,or-8,:) + a(10,5) * var(:,or-9,:) ) * &
+ idel / dy(:,or-4,:)
+ rhs(:,or-3,:) = rhs(:,or-3,:) + &
+ ( a(1,4) * var(:,or,:) + a(2,4) * var(:,or-1,:) + &
+ a(3,4) * var(:,or-2,:) + a(4,4) * var(:,or-3,:) + &
+ a(5,4) * var(:,or-4,:) + a(6,4) * var(:,or-5,:) + &
+ a(7,4) * var(:,or-6,:) + a(8,4) * var(:,or-7,:) + &
+ a(9,4) * var(:,or-8,:) ) * idel / dy(:,or-3,:)
+ rhs(:,or-2,:) = rhs(:,or-2,:) + &
+ ( a(1,3) * var(:,or,:) + a(2,3) * var(:,or-1,:) + &
+ a(3,3) * var(:,or-2,:) + a(4,3) * var(:,or-3,:) + &
+ a(5,3) * var(:,or-4,:) + a(6,3) * var(:,or-5,:) + &
+ a(7,3) * var(:,or-6,:) + a(8,3) * var(:,or-7,:) ) * &
+ idel / dy(:,or-2,:)
+ rhs(:,or-1,:) = rhs(:,or-1,:) + &
+ ( a(1,2) * var(:,or,:) + a(2,2) * var(:,or-1,:) + &
+ a(3,2) * var(:,or-2,:) + a(4,2) * var(:,or-3,:) + &
+ a(5,2) * var(:,or-4,:) + a(6,2) * var(:,or-5,:) + &
+ a(7,2) * var(:,or-6,:) ) * idel / dy(:,or-1,:)
+ rhs(:,or,:) = rhs(:,or,:) + &
+ ( a(1,1) * var(:,or,:) + a(2,1) * var(:,or-1,:) + &
+ a(3,1) * var(:,or-2,:) + a(4,1) * var(:,or-3,:) + &
+ a(5,1) * var(:,or-4,:) + a(6,1) * var(:,or-5,:) ) * &
+ idel / dy(:,or,:)
+!$OMP END PARALLEL WORKSHARE
+
+ jr = or - 8
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,jl:jr,:) = rhs(:,jl:jr,:) + &
+ ( -252.0_wp * var(:,jl:jr,:) + &
+ 210.0_wp * ( var(:,jl-1:jr-1,:) + &
+ var(:,jl+1:jr+1,:) ) - &
+ 120.0_wp * ( var(:,jl-2:jr-2,:) + &
+ var(:,jl+2:jr+2,:) ) + &
+ 45.0_wp * ( var(:,jl-3:jr-3,:) + &
+ var(:,jl+3:jr+3,:) ) - &
+ 10.0_wp * ( var(:,jl-4:jr-4,:) + &
+ var(:,jl+4:jr+4,:) ) + &
+ ( var(:,jl-5:jr-5,:) + &
+ var(:,jl+5:jr+5,:) ) ) * &
+ idel / dy(:,jl:jr,:)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+ if ( zero_derivs_z == 0 ) then
+ call set_coeff ( a )
+
+ idel = epsilon / 1024
+
+ if ( bb(5) == 0 ) then
+ kl = 1 + gsize(3)
+ else
+ ol = offset(5);
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,1+ol) = rhs(:,:,1+ol) + &
+ ( a(1,1) * var(:,:,1+ol) + a(2,1) * var(:,:,2+ol) + &
+ a(3,1) * var(:,:,3+ol) + a(4,1) * var(:,:,4+ol) + &
+ a(5,1) * var(:,:,5+ol) + a(6,1) * var(:,:,6+ol ) ) * &
+ idel / dz(:,:,1+ol)
+ rhs(:,:,2+ol) = rhs(:,:,2+ol) + &
+ ( a(1,2) * var(:,:,1+ol) + a(2,2) * var(:,:,2+ol) + &
+ a(3,2) * var(:,:,3+ol) + a(4,2) * var(:,:,4+ol) + &
+ a(5,2) * var(:,:,5+ol) + a(6,2) * var(:,:,6+ol) + &
+ a(7,2) * var(:,:,7+ol ) ) * idel / dz(:,:,2+ol)
+ rhs(:,:,3+ol) = rhs(:,:,3+ol) + &
+ ( a(1,3) * var(:,:,1+ol) + a(2,3) * var(:,:,2+ol) + &
+ a(3,3) * var(:,:,3+ol) + a(4,3) * var(:,:,4+ol) + &
+ a(5,3) * var(:,:,5+ol) + a(6,3) * var(:,:,6+ol) + &
+ a(7,3) * var(:,:,7+ol) + a(8,3) * var(:,:,8+ol ) ) * &
+ idel / dz(:,:,3+ol)
+ rhs(:,:,4+ol) = rhs(:,:,4+ol) + &
+ ( a(1,4) * var(:,:,1+ol) + a(2,4) * var(:,:,2+ol) + &
+ a(3,4) * var(:,:,3+ol) + a(4,4) * var(:,:,4+ol) + &
+ a(5,4) * var(:,:,5+ol) + a(6,4) * var(:,:,6+ol) + &
+ a(7,4) * var(:,:,7+ol) + a(8,4) * var(:,:,8+ol) + &
+ a(9,4) * var(:,:,9+ol ) ) * idel / dz(:,:,4+ol)
+ rhs(:,:,5+ol) = rhs(:,:,5+ol) + &
+ ( a(1,5) * var(:,:,1+ol) + a(2,5) * var(:,:,2+ol) + &
+ a(3,5) * var(:,:,3+ol) + a(4,5) * var(:,:,4+ol) + &
+ a(5,5) * var(:,:,5+ol) + a(6,5) * var(:,:,6+ol) + &
+ a(7,5) * var(:,:,7+ol) + a(8,5) * var(:,:,8+ol) + &
+ a(9,5) * var(:,:,9+ol) + a(10,5) * var(:,:,10+ol ) ) * &
+ idel / dz(:,:,5+ol)
+ rhs(:,:,6+ol) = rhs(:,:,6+ol) + &
+ ( a(1,6) * var(:,:,1+ol) + a(2,6) * var(:,:,2+ol) + &
+ a(3,6) * var(:,:,3+ol) + a(4,6) * var(:,:,4+ol) + &
+ a(5,6) * var(:,:,5+ol) + a(6,6) * var(:,:,6+ol) + &
+ a(7,6) * var(:,:,7+ol) + a(8,6) * var(:,:,8+ol) + &
+ a(9,6) * var(:,:,9+ol) + a(10,6) * var(:,:,10+ol) + &
+ a(11,6) * var(:,:,11+ol ) ) * idel / dz(:,:,6+ol)
+ rhs(:,:,7+ol) = rhs(:,:,7+ol) + &
+ ( a(2,7) * var(:,:,2+ol) + a(3,7) * var(:,:,3+ol) + &
+ a(4,7) * var(:,:,4+ol) + a(5,7) * var(:,:,5+ol) + &
+ a(6,7) * var(:,:,6+ol) + a(7,7) * var(:,:,7+ol) + &
+ a(8,7) * var(:,:,8+ol) + a(9,7) * var(:,:,9+ol) + &
+ a(10,7) * var(:,:,10+ol) + a(11,7) * var(:,:,11+ol) + &
+ a(12,7) * var(:,:,12+ol ) ) * idel / dz(:,:,7+ol)
+ rhs(:,:,8+ol) = rhs(:,:,8+ol) + &
+ ( a(3,8) * var(:,:,3+ol) + a(4,8) * var(:,:,4+ol) + &
+ a(5,8) * var(:,:,5+ol) + a(6,8) * var(:,:,6+ol) + &
+ a(7,8) * var(:,:,7+ol) + a(8,8) * var(:,:,8+ol) + &
+ a(9,8) * var(:,:,9+ol) + a(10,8) * var(:,:,10+ol) + &
+ a(11,8) * var(:,:,11+ol) + a(12,8) * var(:,:,12+ol) + &
+ a(13,8) * var(:,:,13+ol ) ) * idel / dz(:,:,8+ol)
+!$OMP END PARALLEL WORKSHARE
+
+ kl = 9 + ol
+ end if
+ if ( bb(6) == 0 ) then
+ kr = nk - gsize(3)
+ else
+ or = nk - offset(6)
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,or-7) = rhs(:,:,or-7) + &
+ ( a(3,8) * var(:,:,or-2) + a(4,8) * var(:,:,or-3) + &
+ a(5,8) * var(:,:,or-4) + a(6,8) * var(:,:,or-5) + &
+ a(7,8) * var(:,:,or-6) + a(8,8) * var(:,:,or-7) + &
+ a(9,8) * var(:,:,or-8) + a(10,8) * var(:,:,or-9) + &
+ a(11,8) * var(:,:,or-10) + a(12,8) * var(:,:,or-11) + &
+ a(13,8) * var(:,:,or-12) ) * idel / dz(:,:,or-7)
+ rhs(:,:,or-6) = rhs(:,:,or-6) + &
+ ( a(2,7) * var(:,:,or-1) + a(3,7) * var(:,:,or-2) + &
+ a(4,7) * var(:,:,or-3) + a(5,7) * var(:,:,or-4) + &
+ a(6,7) * var(:,:,or-5) + a(7,7) * var(:,:,or-6) + &
+ a(8,7) * var(:,:,or-7) + a(9,7) * var(:,:,or-8) + &
+ a(10,7) * var(:,:,or-9) + a(11,7) * var(:,:,or-10) + &
+ a(12,7) * var(:,:,or-11) ) * idel / dz(:,:,or-6)
+ rhs(:,:,or-5) = rhs(:,:,or-5) + &
+ ( a(1,6) * var(:,:,or) + a(2,6) * var(:,:,or-1) + &
+ a(3,6) * var(:,:,or-2) + a(4,6) * var(:,:,or-3) + &
+ a(5,6) * var(:,:,or-4) + a(6,6) * var(:,:,or-5) + &
+ a(7,6) * var(:,:,or-6) + a(8,6) * var(:,:,or-7) + &
+ a(9,6) * var(:,:,or-8) + a(10,6) * var(:,:,or-9) + &
+ a(11,6) * var(:,:,or-10) ) * idel / dz(:,:,or-5)
+ rhs(:,:,or-4) = rhs(:,:,or-4) + &
+ ( a(1,5) * var(:,:,or) + a(2,5) * var(:,:,or-1) + &
+ a(3,5) * var(:,:,or-2) + a(4,5) * var(:,:,or-3) + &
+ a(5,5) * var(:,:,or-4) + a(6,5) * var(:,:,or-5) + &
+ a(7,5) * var(:,:,or-6) + a(8,5) * var(:,:,or-7) + &
+ a(9,5) * var(:,:,or-8) + a(10,5) * var(:,:,or-9) ) * &
+ idel / dz(:,:,or-4)
+ rhs(:,:,or-3) = rhs(:,:,or-3) + &
+ ( a(1,4) * var(:,:,or) + a(2,4) * var(:,:,or-1) + &
+ a(3,4) * var(:,:,or-2) + a(4,4) * var(:,:,or-3) + &
+ a(5,4) * var(:,:,or-4) + a(6,4) * var(:,:,or-5) + &
+ a(7,4) * var(:,:,or-6) + a(8,4) * var(:,:,or-7) + &
+ a(9,4) * var(:,:,or-8) ) * idel / dz(:,:,or-3)
+ rhs(:,:,or-2) = rhs(:,:,or-2) + &
+ ( a(1,3) * var(:,:,or) + a(2,3) * var(:,:,or-1) + &
+ a(3,3) * var(:,:,or-2) + a(4,3) * var(:,:,or-3) + &
+ a(5,3) * var(:,:,or-4) + a(6,3) * var(:,:,or-5) + &
+ a(7,3) * var(:,:,or-6) + a(8,3) * var(:,:,or-7) ) * &
+ idel / dz(:,:,or-2)
+ rhs(:,:,or-1) = rhs(:,:,or-1) + &
+ ( a(1,2) * var(:,:,or) + a(2,2) * var(:,:,or-1) + &
+ a(3,2) * var(:,:,or-2) + a(4,2) * var(:,:,or-3) + &
+ a(5,2) * var(:,:,or-4) + a(6,2) * var(:,:,or-5) + &
+ a(7,2) * var(:,:,or-6) ) * idel / dz(:,:,or-1)
+ rhs(:,:,or) = rhs(:,:,or) + &
+ ( a(1,1) * var(:,:,or) + a(2,1) * var(:,:,or-1) + &
+ a(3,1) * var(:,:,or-2) + a(4,1) * var(:,:,or-3) + &
+ a(5,1) * var(:,:,or-4) + a(6,1) * var(:,:,or-5) ) * &
+ idel / dz(:,:,or)
+!$OMP END PARALLEL WORKSHARE
+
+ kr = or - 8
+ end if
+
+!$OMP PARALLEL WORKSHARE
+ rhs(:,:,kl:kr) = rhs(:,:,kl:kr) + &
+ ( -252.0_wp * var(:,:,kl:kr) + &
+ 210.0_wp * ( var(:,:,kl-1:kr-1) + &
+ var(:,:,kl+1:kr+1) ) - &
+ 120.0_wp * ( var(:,:,kl-2:kr-2) + &
+ var(:,:,kl+2:kr+2) ) + &
+ 45.0_wp * ( var(:,:,kl-3:kr-3) + &
+ var(:,:,kl+3:kr+3) ) - &
+ 10.0_wp * ( var(:,:,kl-4:kr-4) + &
+ var(:,:,kl+4:kr+4) ) + &
+ ( var(:,:,kl-5:kr-5) + &
+ var(:,:,kl+5:kr+5) ) ) * &
+ idel / dz(:,:,kl:kr)
+!$OMP END PARALLEL WORKSHARE
+ end if
+
+contains
+
+ subroutine set_coeff ( a )
+
+ implicit none
+
+ CCTK_REAL, dimension(13,8), intent(OUT) :: a
+ CCTK_REAL :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+
+ a(1,1) = -3.3910872088637970174997113084967416241083103770745_wp
+ a(2,1) = 16.955436044318985087498556542483708120541551885372_wp
+ a(3,1) = -33.910872088637970174997113084967416241083103770745_wp
+ a(4,1) = 33.910872088637970174997113084967416241083103770745_wp
+ a(5,1) = -16.955436044318985087498556542483708120541551885372_wp
+ a(6,1) = 3.3910872088637970174997113084967416241083103770745_wp
+ a(7,1) = 0.0_wp
+ a(8,1) = 0.0_wp
+ a(9,1) = 0.0_wp
+ a(10,1) = 0.0_wp
+ a(11,1) = 0.0_wp
+ a(12,1) = 0.0_wp
+ a(13,1) = 0.0_wp
+ a(1,2) = 3.2771399440263630592058029074141137010783820566473_wp
+ a(2,2) = -17.041127708937087907870175118553391245607586694566_wp
+ a(3,2) = 36.048539384289993651263831981555250711862202623121_wp
+ a(4,2) = -39.325679328316356710469634888969364412940584679768_wp
+ a(5,2) = 22.939979608184541414440620351898795907548674396531_wp
+ a(6,2) = -6.5542798880527261184116058148282274021567641132947_wp
+ a(7,2) = 0.65542798880527261184116058148282274021567641132947_wp
+ a(8,2) = 0.0_wp
+ a(9,2) = 0.0_wp
+ a(10,2) = 0.0_wp
+ a(11,2) = 0.0_wp
+ a(12,2) = 0.0_wp
+ a(13,2) = 0.0_wp
+ a(1,3) = -38.842059631038967294446317614758441308222147295410_wp
+ a(2,3) = 213.63132797071432011945474688117142719522181012475_wp
+ a(3,3) = -489.40995135109098791002360194595636048359905592216_wp
+ a(4,3) = 602.05192428110399306391792302875584027744328307885_wp
+ a(5,3) = -427.26265594142864023890949376234285439044362024951_wp
+ a(6,3) = 174.78926833967535282500842926641298588699966282934_wp
+ a(7,3) = -38.842059631038967294446317614758441308222147295410_wp
+ a(8,3) = 3.8842059631038967294446317614758441308222147295410_wp
+ a(9,3) = 0.0_wp
+ a(10,3) = 0.0_wp
+ a(11,3) = 0.0_wp
+ a(12,3) = 0.0_wp
+ a(13,3) = 0.0_wp
+ a(1,4) = 5.5613835719414344378807953109542612676331552744485_wp
+ a(2,4) = -33.368301431648606627284771865725567605798931646691_wp
+ a(3,4) = 86.201445365092233787152327319791049648313906753952_wp
+ a(4,4) = -125.68726872587641829610597402756630464850930920254_wp
+ a(5,4) = 114.00836322479940597655630387456235598647968312619_wp
+ a(6,4) = -66.736602863297213254569543731451135211597863293382_wp
+ a(7,4) = 25.026226073736454970463578899294175704349198735018_wp
+ a(8,4) = -5.5613835719414344378807953109542612676331552744485_wp
+ a(9,4) = 0.55613835719414344378807953109542612676331552744485_wp
+ a(10,4) = 0.0_wp
+ a(11,4) = 0.0_wp
+ a(12,4) = 0.0_wp
+ a(13,4) = 0.0_wp
+ a(1,5) = -12.115101476661536355654081268132755978592914829047_wp
+ a(2,5) = 84.805710336630754489578568876929291850150403803330_wp
+ a(3,5) = -266.53223248655379982438978789892063152904412623904_wp
+ a(4,5) = 496.71916054312299058181733199344299512230950799093_wp
+ a(5,5) = -608.17809412840912505383487966026435012536432441817_wp
+ a(6,5) = 508.83426201978452693747141326157575110090242281998_wp
+ a(7,5) = -290.76243543987687253569795043518614348622995589713_wp
+ a(8,5) = 109.03591328995382720088673141319480380733623346142_wp
+ a(9,5) = -24.230202953323072711308162536265511957185829658094_wp
+ a(10,5) = 2.4230202953323072711308162536265511957185829658094_wp
+ a(11,5) = 0.0_wp
+ a(12,5) = 0.0_wp
+ a(13,5) = 0.0_wp
+ a(1,6) = 0.78217600900123184961735065035840034142603567514089_wp
+ a(2,6) = -7.8217600900123184961735065035840034142603567514089_wp
+ a(3,6) = 35.197920405055433232780779266128015364171605381340_wp
+ a(4,6) = -93.861121080147821954082078043008040971124281016906_wp
+ a(5,6) = 164.25696189025868841964363657526407169946749177959_wp
+ a(6,6) = -197.10835426831042610357236389031688603936099013550_wp
+ a(7,6) = 164.25696189025868841964363657526407169946749177959_wp
+ a(8,6) = -93.861121080147821954082078043008040971124281016906_wp
+ a(9,6) = 35.197920405055433232780779266128015364171605381340_wp
+ a(10,6) = -7.8217600900123184961735065035840034142603567514089_wp
+ a(11,6) = 0.78217600900123184961735065035840034142603567514089_wp
+ a(12,6) = 0.0_wp
+ a(13,6) = 0.0_wp
+ a(1,7) = 0.0_wp
+ a(2,7) = 1.0830767761393601764536458481012280421614377748694_wp
+ a(3,7) = -10.830767761393601764536458481012280421614377748694_wp
+ a(4,7) = 48.738454926271207940414063164555261897264699869122_wp
+ a(5,7) = -129.96921313672322117443750177214736505937253298433_wp
+ a(6,7) = 227.44612298926563705526562810125788885390193272257_wp
+ a(7,7) = -272.93534758711876446631875372150946662468231926708_wp
+ a(8,7) = 227.44612298926563705526562810125788885390193272257_wp
+ a(9,7) = -129.96921313672322117443750177214736505937253298433_wp
+ a(10,7) = 48.738454926271207940414063164555261897264699869122_wp
+ a(11,7) = -10.830767761393601764536458481012280421614377748694_wp
+ a(12,7) = 1.0830767761393601764536458481012280421614377748694_wp
+ a(13,7) = 0.0_wp
+ a(1,8) = 0.0_wp
+ a(2,8) = 0.0_wp
+ a(3,8) = 0.99075245444434671889501396229410272246695863420506_wp
+ a(4,8) = -9.9075245444434671889501396229410272246695863420506_wp
+ a(5,8) = 44.583860449995602350275628303234622511013138539228_wp
+ a(6,8) = -118.89029453332160626740167547529232669603503610461_wp
+ a(7,8) = 208.05801543331281096795293208176157171806131318306_wp
+ a(8,8) = -249.66961851997537316154351849811388606167357581967_wp
+ a(9,8) = 208.05801543331281096795293208176157171806131318306_wp
+ a(10,8) = -118.89029453332160626740167547529232669603503610461_wp
+ a(11,8) = 44.583860449995602350275628303234622511013138539228_wp
+ a(12,8) = -9.9075245444434671889501396229410272246695863420506_wp
+ a(13,8) = 0.99075245444434671889501396229410272246695863420506_wp
+
+ end subroutine set_coeff
+
+end subroutine dissipation_8_4_delta
diff --git a/src/dissipation.c b/src/dissipation.c
index 548f522..690fce7 100644
--- a/src/dissipation.c
+++ b/src/dissipation.c
@@ -85,6 +85,18 @@ void CCTK_FCALL CCTK_FNAME(dissipation_2_1_alt) (const CCTK_REAL *var,
const CCTK_REAL *dx,
const CCTK_REAL *epsdis,
CCTK_REAL *rhs);
+void CCTK_FCALL CCTK_FNAME(dissipation_2_1_delta) (const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ const CCTK_INT *bbox,
+ const CCTK_INT *gsize,
+ const CCTK_INT *offset,
+ const CCTK_REAL *dx,
+ const CCTK_REAL *dy,
+ const CCTK_REAL *dz,
+ const CCTK_REAL *epsdis,
+ CCTK_REAL *rhs);
void CCTK_FCALL CCTK_FNAME(dissipation_4_2) (const CCTK_REAL *var,
const CCTK_INT *ni,
const CCTK_INT *nj,
@@ -105,6 +117,18 @@ void CCTK_FCALL CCTK_FNAME(dissipation_4_2_alt) (const CCTK_REAL *var,
const CCTK_REAL *dx,
const CCTK_REAL *epsdis,
CCTK_REAL *rhs);
+void CCTK_FCALL CCTK_FNAME(dissipation_4_2_delta) (const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ const CCTK_INT *bbox,
+ const CCTK_INT *gsize,
+ const CCTK_INT *offset,
+ const CCTK_REAL *dx,
+ const CCTK_REAL *dy,
+ const CCTK_REAL *dz,
+ const CCTK_REAL *epsdis,
+ CCTK_REAL *rhs);
void CCTK_FCALL CCTK_FNAME(dissipation_6_3) (const CCTK_REAL *var,
const CCTK_INT *ni,
const CCTK_INT *nj,
@@ -125,6 +149,18 @@ void CCTK_FCALL CCTK_FNAME(dissipation_6_3_alt) (const CCTK_REAL *var,
const CCTK_REAL *dx,
const CCTK_REAL *epsdis,
CCTK_REAL *rhs);
+void CCTK_FCALL CCTK_FNAME(dissipation_6_3_delta) (const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ const CCTK_INT *bbox,
+ const CCTK_INT *gsize,
+ const CCTK_INT *offset,
+ const CCTK_REAL *dx,
+ const CCTK_REAL *dy,
+ const CCTK_REAL *dz,
+ const CCTK_REAL *epsdis,
+ CCTK_REAL *rhs);
void CCTK_FCALL CCTK_FNAME(dissipation_8_4) (const CCTK_REAL *var,
const CCTK_INT *ni,
const CCTK_INT *nj,
@@ -145,6 +181,18 @@ void CCTK_FCALL CCTK_FNAME(dissipation_8_4_alt) (const CCTK_REAL *var,
const CCTK_REAL *dx,
const CCTK_REAL *epsdis,
CCTK_REAL *rhs);
+void CCTK_FCALL CCTK_FNAME(dissipation_8_4_delta) (const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ const CCTK_INT *bbox,
+ const CCTK_INT *gsize,
+ const CCTK_INT *offset,
+ const CCTK_REAL *dx,
+ const CCTK_REAL *dy,
+ const CCTK_REAL *dz,
+ const CCTK_REAL *epsdis,
+ CCTK_REAL *rhs);
void get_shiftout ( const CCTK_POINTER_TO_CONST cctkGH_, CCTK_INT *offset);
void CCTK_FCALL CCTK_FNAME(SBP_Poisoning) ( const CCTK_INT *ni,
@@ -264,9 +312,16 @@ apply (int const varindex, char const * const optstring, void * const arg)
(varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
bbox, gsize, offset, dx, &epsdis, rhsptr);
} else {
- CCTK_FNAME(dissipation_2_1_alt)
- (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
- bbox, gsize, offset, dx, &epsdis, rhsptr);
+ if (!use_variable_deltas) {
+ CCTK_FNAME(dissipation_2_1_alt)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset, dx, &epsdis, rhsptr);
+ } else {
+ CCTK_FNAME(dissipation_2_1_delta)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset,
+ sbp_dx, sbp_dy, sbp_dz, &epsdis, rhsptr);
+ }
}
break;
}
@@ -276,9 +331,16 @@ apply (int const varindex, char const * const optstring, void * const arg)
(varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
bbox, gsize, offset, dx, &epsdis, rhsptr);
} else {
- CCTK_FNAME(dissipation_4_2_alt)
- (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
- bbox, gsize, offset, dx, &epsdis, rhsptr);
+ if (!use_variable_deltas) {
+ CCTK_FNAME(dissipation_4_2_alt)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset, dx, &epsdis, rhsptr);
+ } else {
+ CCTK_FNAME(dissipation_4_2_delta)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset,
+ sbp_dx, sbp_dy, sbp_dz, &epsdis, rhsptr);
+ }
}
break;
}
@@ -288,9 +350,16 @@ apply (int const varindex, char const * const optstring, void * const arg)
(varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
bbox, gsize, offset, dx, &epsdis, rhsptr);
} else {
- CCTK_FNAME(dissipation_6_3_alt)
- (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
- bbox, gsize, offset, dx, &epsdis, rhsptr);
+ if (!use_variable_deltas) {
+ CCTK_FNAME(dissipation_6_3_alt)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset, dx, &epsdis, rhsptr);
+ } else {
+ CCTK_FNAME(dissipation_6_3_delta)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset,
+ sbp_dx, sbp_dy, sbp_dz, &epsdis, rhsptr);
+ }
}
break;
}
@@ -300,9 +369,16 @@ apply (int const varindex, char const * const optstring, void * const arg)
(varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
bbox, gsize, offset, dx, &epsdis, rhsptr);
} else {
- CCTK_FNAME(dissipation_8_4_alt)
- (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
- bbox, gsize, offset, dx, &epsdis, rhsptr);
+ if (!use_variable_deltas) {
+ CCTK_FNAME(dissipation_8_4_alt)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset, dx, &epsdis, rhsptr);
+ } else {
+ CCTK_FNAME(dissipation_8_4_delta)
+ (varptr, &cctk_lsh[0], &cctk_lsh[1], &cctk_lsh[2],
+ bbox, gsize, offset,
+ sbp_dx, sbp_dy, sbp_dz, &epsdis, rhsptr);
+ }
}
break;
}
diff --git a/src/make.code.defn b/src/make.code.defn
index 150ab07..8eaaf16 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -41,6 +41,7 @@ SRCS = call_derivs.c \
Dissipation_6_5_min_err_coeff.F90 \
Dissipation_6_5_min_err_coeff_alt.F90 \
dissipation_coeff.F90 \
+ DeltaInitial.F90 \
get_coeffs.c \
get_up_coeffs.c \
get_coeffs2.c \
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-08 10:09:07 +0000