aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--interface.ccl14
-rw-r--r--param.ccl6
-rw-r--r--src/Derivatives_4_3.F90242
-rw-r--r--src/Get_Coeff.F9032
-rw-r--r--src/call_derivs.c58
-rw-r--r--src/call_derivs_name.c132
-rw-r--r--src/call_norm_mask.c24
-rw-r--r--src/make.code.defn2
-rw-r--r--src/set_norm_mask.F90124
9 files changed, 558 insertions, 76 deletions
diff --git a/interface.ccl b/interface.ccl
index 9eef8aa..c05a3c1 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -25,6 +25,13 @@ SUBROUTINE Diff_gv ( CCTK_POINTER IN cctkGH, \
CCTK_REAL OUT ARRAY dvar )
PROVIDES FUNCTION Diff_gv WITH DiffGv LANGUAGE C
+CCTK_REAL FUNCTION GetScalProdCoeff ()
+PROVIDES FUNCTION GetScalProdCoeff WITH GetCoeff LANGUAGE Fortran
+
+SUBROUTINE SetNormMask ( CCTK_POINTER IN cctkGH, \
+ CCTK_REAL OUT ARRAY mask )
+PROVIDES FUNCTION SetNormMask WITH SetNMask LANGUAGE C
+
CCTK_INT FUNCTION GetDomainSpecification \
(CCTK_INT IN size, \
CCTK_REAL OUT ARRAY physical_min, \
@@ -35,3 +42,10 @@ CCTK_INT FUNCTION GetDomainSpecification \
CCTK_REAL OUT ARRAY exterior_max, \
CCTK_REAL OUT ARRAY spacing)
USES FUNCTION GetDomainSpecification
+
+CCTK_INT FUNCTION \
+ SymmetryTableHandleForGrid (CCTK_POINTER_TO_CONST IN cctkGH)
+REQUIRES FUNCTION SymmetryTableHandleForGrid
+
+CCTK_INT FUNCTION SymmetryHandleOfName (CCTK_STRING IN sym_name)
+REQUIRES FUNCTION SymmetryHandleOfName
diff --git a/param.ccl b/param.ccl
index 4881ebf..285f928 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,6 +1,12 @@
# Parameter definitions for thorn SummationByParts
# $Header$
+KEYWORD norm_type "Type of norm"
+{
+ "Diagonal" :: "Diagonal norm"
+ "Full restricted" :: "Full resricted norm"
+} "Diagonal"
+
INT order "Order of accuracy" STEERABLE=always
{
2:8:2 :: ""
diff --git a/src/Derivatives_4_3.F90 b/src/Derivatives_4_3.F90
new file mode 100644
index 0000000..011a074
--- /dev/null
+++ b/src/Derivatives_4_3.F90
@@ -0,0 +1,242 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine deriv_gf_4_3 ( var, ni, nj, nk, dir, bb, gsize, delta, dvar )
+
+ implicit none
+
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_REAL, parameter :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+ CCTK_INT, intent(IN) :: ni, nj, nk
+ CCTK_REAL, dimension(ni,nj,nk), intent(IN) :: var
+ CCTK_INT, intent(IN) :: dir
+ CCTK_INT, intent(IN) :: bb(2)
+ CCTK_INT, intent(IN) :: gsize
+ CCTK_REAL, intent(IN) :: delta
+ CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: dvar
+
+ CCTK_REAL, dimension(2), save :: a
+ CCTK_REAL, dimension(7,5), save :: q
+ CCTK_REAL :: idel
+ CCTK_REAL :: f0, f1, f2, f3
+
+ CCTK_INT :: il, ir, jl, jr, kl, kr
+
+ logical, save :: first = .true.
+
+ if ( first ) then
+ a(1) = 2.0_wp/3.0_wp; a(2) = -1.0_wp/12.0_wp
+
+ f0 = sqrt(26116897.0_wp)
+ f1 = -56764003702447356523.0_wp + 8154993476273221.0_wp * f0
+ f2 = -55804550303.0_wp + 9650225.0_wp * f0
+ f3 = 3262210757.0_wp + 271861.0_wp * f0
+
+ q(1,1) = -11._wp/6.0_wp; q(2,1) = 3.0_wp
+ q(3,1) = -3.0_wp/2.0_wp; q(4,1) = 1.0_wp/3.0_wp
+ q(5,1) = zero; q(6,1) = zero; q(7,1) = zero
+
+ q(1,2) = -24.0_wp * ( -779042810827742869.0_wp + &
+ 104535124033147.0_wp * f0 ) / f1
+ q(2,2) = -( -176530817412806109689.0_wp + &
+ 29768274816875927.0_wp * f0 ) / ( 6.0_wp * f1 )
+ q(3,2) = 343.0_wp * ( -171079116122226871.0_wp + &
+ 27975630462649.0_wp * f0 ) / f1
+ q(4,2) = -3.0_wp * ( -7475554291248533227.0_wp + &
+ 1648464218793925.0_wp * f0 ) / ( 2.0_wp * f1 )
+ q(5,2) = ( -2383792768180030915.0_wp + &
+ 1179620587812973.0_wp * f0 ) / ( 3.0_wp * f1 )
+ q(6,2) = -1232.0_wp * ( -115724529581315.0_wp + 37280576429.0_wp * f0 ) / f1
+ q(7,2) = zero
+
+ q(1,3) = -12.0_wp * ( -380966843.0_wp + 86315.0_wp * f0 ) / f2
+ q(2,3) = ( 5024933015.0_wp + 2010631.0_wp * f0 ) / ( 3.0_wp * f2 )
+ q(3,3) = -231.0_wp * ( -431968921.0_wp + 86711.0_wp * f0 ) / ( 2.0_wp * f2 )
+ q(4,3) = ( -65931742559.0_wp + 12256337.0_wp * f0 ) / f2
+ q(5,3) = -( -50597298167.0_wp + 9716873.0_wp * f0 ) / ( 6.0_wp * f2 )
+ q(6,3) = -88.0_wp * ( -15453061.0_wp + 2911.0_wp * f0 ) / f2
+ q(7,3) = zero
+
+ q(1,4) = 48.0_wp * ( -56020909845192541.0_wp + &
+ 9790180507043.0_wp * f0 ) / f1
+ q(2,4) = ( -9918249049237586011.0_wp + &
+ 1463702013196501.0_wp * f0 ) / ( 6.0_wp * f1)
+ q(3,4) = -13.0_wp * ( -4130451756851441723.0_wp + &
+ 664278707201077.0_wp * f0 ) / f1
+ q(4,4) = 3.0_wp * ( -26937108467782666617.0_wp + &
+ 5169063172799767.0_wp * f0 ) / ( 2.0_wp * f1 )
+ q(5,4) = -( 6548308508012371315.0_wp + &
+ 3968886380989379.0_wp * f0 ) / ( 3.0_wp * f1 )
+ q(6,4) = 88.0_wp * ( -91337851897923397.0_wp + &
+ 19696768305507.0_wp * f0 ) / f1
+ q(7,4) = 242.0_wp * ( -120683.0_wp + 15.0_wp * f0 ) / f3
+
+ q(1,5) = 264.0_wp * ( -120683.0_wp + 15.0_wp * f0 ) / f3
+ q(2,5) = ( -43118111.0_wp + 23357.0_wp * f0 ) / ( 3.0_wp * f3 )
+ q(3,5) = -47.0_wp * ( -28770085.0_wp + 2259.0_wp * f0 ) / ( 2.0_wp * f3 )
+ q(4,5) = -3.0_wp * ( 1003619433.0_wp + 11777.0_wp * f0 ) / f3
+ q(5,5) = -11.0_wp * ( -384168269.0_wp + 65747.0_wp * f0 ) / ( 6.0_wp * f3 )
+ q(6,5) = 22.0_wp * ( 87290207.0_wp + 10221.0_wp * f0 ) / f3
+ q(7,5) = -66.0_wp * ( 3692405.0_wp + 419.0_wp * f0 ) / f3
+ first = .false.
+ end if
+
+ idel = 1.0_wp / delta
+
+ direction: select case (dir)
+ case (0) direction
+ if ( bb(1) == 0 ) then
+ il = 1 + gsize
+ else
+ dvar(1,:,:) = ( q(1,1) * var(1,:,:) + q(2,1) * var(2,:,:) + &
+ q(3,1) * var(3,:,:) + q(4,1) * var(4,:,:) ) * idel
+ dvar(2,:,:) = ( q(1,2) * var(1,:,:) + q(2,2) * var(2,:,:) + &
+ q(3,2) * var(3,:,:) + q(4,2) * var(4,:,:) + &
+ q(5,2) * var(5,:,:) + q(6,2) * var(6,:,:) ) * idel
+ dvar(3,:,:) = ( q(1,3) * var(1,:,:) + q(2,3) * var(2,:,:) + &
+ q(3,3) * var(3,:,:) + q(4,3) * var(4,:,:) + &
+ q(5,3) * var(5,:,:) + q(6,3) * var(6,:,:) ) * idel
+ dvar(4,:,:) = ( q(1,4) * var(1,:,:) + q(2,4) * var(2,:,:) + &
+ q(3,4) * var(3,:,:) + q(4,4) * var(4,:,:) + &
+ q(5,4) * var(5,:,:) + q(6,4) * var(6,:,:) + &
+ q(7,4) * var(7,:,:) ) * idel
+ dvar(5,:,:) = ( q(1,5) * var(1,:,:) + q(2,5) * var(2,:,:) + &
+ q(3,5) * var(3,:,:) + q(4,5) * var(4,:,:) + &
+ q(5,5) * var(5,:,:) + q(6,5) * var(6,:,:) + &
+ q(7,5) * var(7,:,:) ) * idel
+ il = 6
+ end if
+ if ( bb(2) == 0 ) then
+ ir = ni - gsize
+ else
+ dvar(ni-4,:,:) = - ( q(1,5) * var(ni,:,:) + q(2,5) * var(ni-1,:,:) + &
+ q(3,5) * var(ni-2,:,:) + q(4,5) * var(ni-3,:,:) + &
+ q(5,5) * var(ni-4,:,:) + q(6,5) * var(ni-5,:,:) + &
+ q(7,5) * var(ni-6,:,:) ) * idel
+ dvar(ni-3,:,:) = - ( q(1,4) * var(ni,:,:) + q(2,4) * var(ni-1,:,:) + &
+ q(3,4) * var(ni-2,:,:) + q(4,4) * var(ni-3,:,:) + &
+ q(5,4) * var(ni-4,:,:) + q(6,4) * var(ni-5,:,:) + &
+ q(7,4) * var(ni-6,:,:) ) * idel
+ dvar(ni-2,:,:) = - ( q(1,3) * var(ni,:,:) + q(2,3) * var(ni-1,:,:) + &
+ q(3,3) * var(ni-2,:,:) + q(4,3) * var(ni-3,:,:) + &
+ q(5,3) * var(ni-4,:,:) + &
+ q(6,3) * var(ni-5,:,:) ) * idel
+ dvar(ni-1,:,:) = - ( q(1,2) * var(ni,:,:) + q(2,2) * var(ni-1,:,:) + &
+ q(3,2) * var(ni-2,:,:) + q(4,2) * var(ni-3,:,:) + &
+ q(5,2) * var(ni-4,:,:) + &
+ q(6,2) * var(ni-5,:,:) ) * idel
+ dvar(ni,:,:) = - ( q(1,1) * var(ni,:,:) + q(2,1) * var(ni-1,:,:) + &
+ q(3,1) * var(ni-2,:,:) + &
+ q(4,1) * var(ni-3,:,:) ) * idel
+ ir = ni - 5
+ end if
+ dvar(il:ir,:,:) = ( a(1) * ( var(il+1:ir+1,:,:) - &
+ var(il-1:ir-1,:,:) ) + &
+ a(2) * ( var(il+2:ir+2,:,:) - &
+ var(il-2:ir-2,:,:) ) ) * idel
+ case (1) direction
+ if ( bb(1) == 0 ) then
+ jl = 1 + gsize
+ else
+ dvar(:,1,:) = ( q(1,1) * var(:,1,:) + q(2,1) * var(:,2,:) + &
+ q(3,1) * var(:,3,:) + q(4,1) * var(:,4,:) ) * idel
+ dvar(:,2,:) = ( q(1,2) * var(:,1,:) + q(2,2) * var(:,2,:) + &
+ q(3,2) * var(:,3,:) + q(4,2) * var(:,4,:) + &
+ q(5,2) * var(:,5,:) + q(6,2) * var(:,6,:) ) * idel
+ dvar(:,3,:) = ( q(1,3) * var(:,1,:) + q(2,3) * var(:,2,:) + &
+ q(3,3) * var(:,3,:) + q(4,3) * var(:,4,:) + &
+ q(5,3) * var(:,5,:) + q(6,3) * var(:,6,:) ) * idel
+ dvar(:,4,:) = ( q(1,4) * var(:,1,:) + q(2,4) * var(:,2,:) + &
+ q(3,4) * var(:,3,:) + q(4,4) * var(:,4,:) + &
+ q(5,4) * var(:,5,:) + q(6,4) * var(:,6,:) + &
+ q(7,4) * var(:,7,:) ) * idel
+ dvar(:,5,:) = ( q(1,5) * var(:,1,:) + q(2,5) * var(:,2,:) + &
+ q(3,5) * var(:,3,:) + q(4,5) * var(:,4,:) + &
+ q(5,5) * var(:,5,:) + q(6,5) * var(:,6,:) + &
+ q(7,5) * var(:,7,:) ) * idel
+ jl = 6
+ end if
+ if ( bb(2) == 0 ) then
+ jr = nj - gsize
+ else
+ dvar(:,nj-4,:) = - ( q(1,5) * var(:,nj,:) + q(2,5) * var(:,nj-1,:) + &
+ q(3,5) * var(:,nj-2,:) + q(4,5) * var(:,nj-3,:) + &
+ q(5,5) * var(:,nj-4,:) + q(6,5) * var(:,nj-5,:) + &
+ q(7,5) * var(:,nj-6,:) ) * idel
+ dvar(:,nj-3,:) = - ( q(1,4) * var(:,nj,:) + q(2,4) * var(:,nj-1,:) + &
+ q(3,4) * var(:,nj-2,:) + q(4,4) * var(:,nj-3,:) + &
+ q(5,4) * var(:,nj-4,:) + q(6,4) * var(:,nj-5,:) + &
+ q(7,4) * var(:,nj-6,:) ) * idel
+ dvar(:,nj-2,:) = - ( q(1,3) * var(:,nj,:) + q(2,3) * var(:,nj-1,:) + &
+ q(3,3) * var(:,nj-2,:) + q(4,3) * var(:,nj-3,:) + &
+ q(5,3) * var(:,nj-4,:) + &
+ q(6,3) * var(:,nj-5,:) ) * idel
+ dvar(:,nj-1,:) = - ( q(1,2) * var(:,nj,:) + q(2,2) * var(:,nj-1,:) + &
+ q(3,2) * var(:,nj-2,:) + q(4,2) * var(:,nj-3,:) + &
+ q(5,2) * var(:,nj-4,:) + &
+ q(6,2) * var(:,nj-5,:) ) * idel
+ dvar(:,nj,:) = - ( q(1,1) * var(:,nj,:) + q(2,1) * var(:,nj-1,:) + &
+ q(3,1) * var(:,nj-2,:) + &
+ q(4,1) * var(:,nj-3,:) ) * idel
+ jr = nj - 5
+ end if
+ dvar(:,jl:jr,:) = ( a(1) * ( var(:,jl+1:jr+1,:) - &
+ var(:,jl-1:jr-1,:) ) + &
+ a(2) * ( var(:,jl+2:jr+2,:) - &
+ var(:,jl-2:jr-2,:) ) ) * idel
+ case (2) direction
+ if ( bb(1) == 0 ) then
+ kl = 1 + gsize
+ else
+ dvar(:,:,1) = ( q(1,1) * var(:,:,1) + q(2,1) * var(:,:,2) + &
+ q(3,1) * var(:,:,3) + q(4,1) * var(:,:,4) ) * idel
+ dvar(:,:,2) = ( q(1,2) * var(:,:,1) + q(2,2) * var(:,:,2) + &
+ q(3,2) * var(:,:,3) + q(4,2) * var(:,:,4) + &
+ q(5,2) * var(:,:,5) + q(6,2) * var(:,:,6) ) * idel
+ dvar(:,:,3) = ( q(1,3) * var(:,:,1) + q(2,3) * var(:,:,2) + &
+ q(3,3) * var(:,:,3) + q(4,3) * var(:,:,4) + &
+ q(5,3) * var(:,:,5) + q(6,3) * var(:,:,6) ) * idel
+ dvar(:,:,4) = ( q(1,4) * var(:,:,1) + q(2,4) * var(:,:,2) + &
+ q(3,4) * var(:,:,3) + q(4,4) * var(:,:,4) + &
+ q(5,4) * var(:,:,5) + q(6,4) * var(:,:,6) + &
+ q(7,4) * var(:,:,7) ) * idel
+ dvar(:,:,5) = ( q(1,5) * var(:,:,1) + q(2,5) * var(:,:,2) + &
+ q(3,5) * var(:,:,3) + q(4,5) * var(:,:,4) + &
+ q(5,5) * var(:,:,5) + q(6,5) * var(:,:,6) + &
+ q(7,5) * var(:,:,7) ) * idel
+ kl = 6
+ end if
+ if ( bb(2) == 0 ) then
+ kr = nk - gsize
+ else
+ dvar(:,:,nk-4) = - ( q(1,5) * var(:,:,nk) + q(2,5) * var(:,:,nk-1) + &
+ q(3,5) * var(:,:,nk-2) + q(4,5) * var(:,:,nk-3) + &
+ q(5,5) * var(:,:,nk-4) + q(6,5) * var(:,:,nk-5) + &
+ q(7,5) * var(:,:,nk-6) ) * idel
+ dvar(:,:,nk-3) = - ( q(1,4) * var(:,:,nk) + q(2,4) * var(:,:,nk-1) + &
+ q(3,4) * var(:,:,nk-2) + q(4,4) * var(:,:,nk-3) + &
+ q(5,4) * var(:,:,nk-4) + q(6,4) * var(:,:,nk-5) + &
+ q(7,4) * var(:,:,nk-6) ) * idel
+ dvar(:,:,nk-2) = - ( q(1,3) * var(:,:,nk) + q(2,3) * var(:,:,nk-1) + &
+ q(3,3) * var(:,:,nk-2) + q(4,3) * var(:,:,nk-3) + &
+ q(5,3) * var(:,:,nk-4) + &
+ q(6,3) * var(:,:,nk-5) ) * idel
+ dvar(:,:,nk-1) = - ( q(1,2) * var(:,:,nk) + q(2,2) * var(:,:,nk-1) + &
+ q(3,2) * var(:,:,nk-2) + q(4,2) * var(:,:,nk-3) + &
+ q(5,2) * var(:,:,nk-4) + &
+ q(6,2) * var(:,:,nk-5) ) * idel
+ dvar(:,:,nk) = - ( q(1,1) * var(:,:,nk) + q(2,1) * var(:,:,nk-1) + &
+ q(3,1) * var(:,:,nk-2) + &
+ q(4,1) * var(:,:,nk-3) ) * idel
+ kr = nk - 5
+ end if
+ dvar(:,:,kl:kr) = ( a(1) * ( var(:,:,kl+1:kr+1) - &
+ var(:,:,kl-1:kr-1) ) + &
+ a(2) * ( var(:,:,kl+2:kr+2) - &
+ var(:,:,kl-2:kr-2) ) ) * idel
+ end select direction
+end subroutine deriv_gf_4_3
diff --git a/src/Get_Coeff.F90 b/src/Get_Coeff.F90
new file mode 100644
index 0000000..7b3fcf5
--- /dev/null
+++ b/src/Get_Coeff.F90
@@ -0,0 +1,32 @@
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+CCTK_REAL function GetCoeff ()
+
+ implicit none
+
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_REAL, parameter :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+
+ if ( CCTK_EQUALS(norm_type,'Diagonal') ) then
+ select case (order)
+ case (2)
+ GetCoeff = 1.0_wp / 2.0_wp
+ case (4)
+ GetCoeff = 17.0_wp/48.0_wp
+ case (6)
+ GetCoeff = 13649.0_wp/43200.0_wp
+ case (8)
+ GetCoeff = 1498139.0_wp/5080320.0_wp
+ end select
+ else
+ select case (order)
+ case(4)
+ GetCoeff = 3.0_wp/11.0_wp
+ end select
+ end if
+end function GetCoeff
diff --git a/src/call_derivs.c b/src/call_derivs.c
index e8c051a..f8807c2 100644
--- a/src/call_derivs.c
+++ b/src/call_derivs.c
@@ -33,6 +33,15 @@ void DiffGv ( const CCTK_POINTER cctkGH_, const CCTK_INT dir,
const CCTK_INT *gsize,
const CCTK_REAL *delta,
CCTK_REAL *dvar);
+ void CCTK_FCALL CCTK_FNAME(deriv_gf_4_3)(const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ const CCTK_INT *dir,
+ const CCTK_INT *bb,
+ const CCTK_INT *gsize,
+ const CCTK_REAL *delta,
+ CCTK_REAL *dvar);
void CCTK_FCALL CCTK_FNAME(deriv_gf_6_3)(const CCTK_REAL *var,
const CCTK_INT *ni,
const CCTK_INT *nj,
@@ -78,24 +87,35 @@ void DiffGv ( const CCTK_POINTER cctkGH_, const CCTK_INT dir,
assert(0);
}
- switch(order) {
- case 2: {
- CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
- break;
- }
- case 4: {
- CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
- break;
- }
- case 6: {
- CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
- break;
- }
- case 8: {
- CCTK_FNAME(deriv_gf_8_4)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
- break;
- }
- default:
- assert(0);
+ if ( CCTK_Equals(norm_type,"Diagonal") ) {
+ switch(order) {
+ case 2: {
+ CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ case 4: {
+ CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ case 6: {
+ CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ case 8: {
+ CCTK_FNAME(deriv_gf_8_4)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ default:
+ assert(0);
+ }
+ } else {
+ switch(order) {
+ case 4: {
+ CCTK_FNAME(deriv_gf_4_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ default:
+ assert(0);
+ }
}
}
diff --git a/src/call_derivs_name.c b/src/call_derivs_name.c
index c696e4e..893aa43 100644
--- a/src/call_derivs_name.c
+++ b/src/call_derivs_name.c
@@ -1,82 +1,87 @@
#include "cctk.h"
#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
#include <assert.h>
-void DiffGf ( const CCTK_POINTER cctkGH, const CCTK_INT dir,
+void DiffGf ( const CCTK_POINTER cctkGH_, const CCTK_INT dir,
const char *var_name, const char *dvar_name )
{
+ cGH const * restrict const cctkGH = cctkGH_;
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_ARGUMENTS
+
CCTK_REAL *var, *dvar;
CCTK_INT ni, nj, nk, gsize, ic;
CCTK_REAL delta;
- CCTK_REAL phys_min[3], phys_max[3];
- CCTK_REAL int_min[3], int_max[3];
- CCTK_REAL ext_min[3], ext_max[3];
- CCTK_REAL spacing[3];
CCTK_INT ierr;
CCTK_INT lsh[3], bbox[6], bb[2], nghostzones[3];
- void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(CCTK_REAL *var, const CCTK_INT *ni,
- const CCTK_INT *nj, const CCTK_INT *nk,
+ void CCTK_FCALL CCTK_FNAME(deriv_gf_2_1)(const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
const CCTK_INT *dir,
const CCTK_INT *bb,
const CCTK_INT *gsize,
const CCTK_REAL *delta,
CCTK_REAL *dvar);
- void CCTK_FCALL CCTK_FNAME(deriv_gf_4_2)(CCTK_REAL *var, const CCTK_INT *ni,
- const CCTK_INT *nj, const CCTK_INT *nk,
+ void CCTK_FCALL CCTK_FNAME(deriv_gf_4_2)(const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
const CCTK_INT *dir,
const CCTK_INT *bb,
const CCTK_INT *gsize,
const CCTK_REAL *delta,
CCTK_REAL *dvar);
- void CCTK_FCALL CCTK_FNAME(deriv_gf_6_3)(CCTK_REAL *var, const CCTK_INT *ni,
- const CCTK_INT *nj, const CCTK_INT *nk,
+ void CCTK_FCALL CCTK_FNAME(deriv_gf_4_3)(const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ const CCTK_INT *dir,
+ const CCTK_INT *bb,
+ const CCTK_INT *gsize,
+ const CCTK_REAL *delta,
+ CCTK_REAL *dvar);
+ void CCTK_FCALL CCTK_FNAME(deriv_gf_6_3)(const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ const CCTK_INT *dir,
+ const CCTK_INT *bb,
+ const CCTK_INT *gsize,
+ const CCTK_REAL *delta,
+ CCTK_REAL *dvar);
+ void CCTK_FCALL CCTK_FNAME(deriv_gf_8_4)(const CCTK_REAL *var,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
const CCTK_INT *dir,
const CCTK_INT *bb,
const CCTK_INT *gsize,
const CCTK_REAL *delta,
CCTK_REAL *dvar);
- DECLARE_CCTK_PARAMETERS
- if (CCTK_GrouplshVN(cctkGH, 3, lsh, var_name) < 0)
- {
- CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Error getting local size of grid for variable '%s'",var_name);
- }
- if (CCTK_GroupbboxVN(cctkGH, 6, bbox, var_name) < 0)
- {
- CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Error getting bounding box for variable '%s'",var_name);
- }
- if (CCTK_GroupnghostzonesVN(cctkGH, 3, nghostzones, var_name) < 0)
- {
- CCTK_VWarn(1, __LINE__, __FILE__, CCTK_THORNSTRING,
- "Error getting number of ghostzones for variable '%s'",var_name);
- }
-
- ni = lsh[0]; nj = lsh[1]; nk = lsh[2];
-
- ierr = GetDomainSpecification ( 3, phys_min, phys_max, int_min, int_max,
- ext_min, ext_max, spacing );
+ ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
switch(dir) {
case 0: {
- delta = spacing[0];
- bb[0] = bbox[0]; bb[1] = bbox[1];
- gsize = nghostzones[0];
+ delta = CCTK_DELTA_SPACE(0);
+ bb[0] = cctk_bbox[0]; bb[1] = cctk_bbox[1];
+ gsize = cctk_nghostzones[0];
break;
}
case 1: {
- delta = spacing[1];
- bb[0] = bbox[2]; bb[1] = bbox[3];
- gsize = nghostzones[1];
+ delta = CCTK_DELTA_SPACE(1);
+ bb[0] = cctk_bbox[2]; bb[1] = cctk_bbox[3];
+ gsize = cctk_nghostzones[1];
break;
}
case 2: {
- delta = spacing[2];
- bb[0] = bbox[4]; bb[1] = bbox[5];
- gsize = nghostzones[2];
+ delta = CCTK_DELTA_SPACE(2);
+ bb[0] = cctk_bbox[4]; bb[1] = cctk_bbox[5];
+ gsize = cctk_nghostzones[2];
break;
}
default:
@@ -86,20 +91,35 @@ void DiffGf ( const CCTK_POINTER cctkGH, const CCTK_INT dir,
var = (CCTK_REAL *)(CCTK_VarDataPtr(cctkGH,0,var_name));
dvar = (CCTK_REAL *)(CCTK_VarDataPtr(cctkGH,0,dvar_name));
- switch(order) {
- case 2: {
- CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
- break;
- }
- case 4: {
- CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
- break;
- }
- case 6: {
- CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
- break;
- }
- default:
- assert(0);
+ if ( CCTK_Equals(norm_type,"Diagonal") ) {
+ switch(order) {
+ case 2: {
+ CCTK_FNAME(deriv_gf_2_1)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ case 4: {
+ CCTK_FNAME(deriv_gf_4_2)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ case 6: {
+ CCTK_FNAME(deriv_gf_6_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ case 8: {
+ CCTK_FNAME(deriv_gf_8_4)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ default:
+ assert(0);
+ }
+ } else {
+ switch(order) {
+ case 4: {
+ CCTK_FNAME(deriv_gf_4_3)(var,&ni,&nj,&nk,&dir,bb,&gsize,&delta,dvar);
+ break;
+ }
+ default:
+ assert(0);
+ }
}
}
diff --git a/src/call_norm_mask.c b/src/call_norm_mask.c
new file mode 100644
index 0000000..5e7996a
--- /dev/null
+++ b/src/call_norm_mask.c
@@ -0,0 +1,24 @@
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+#include <assert.h>
+
+void SetNMask ( const CCTK_POINTER cctkGH_, CCTK_REAL *mask )
+{
+ cGH const * restrict const cctkGH = cctkGH_;
+ DECLARE_CCTK_PARAMETERS
+ DECLARE_CCTK_ARGUMENTS
+
+ CCTK_INT ni, nj, nk;
+ void CCTK_FCALL CCTK_FNAME(sbp_set_norm_mask)(const CCTK_POINTER cctkGH,
+ const CCTK_INT *ni,
+ const CCTK_INT *nj,
+ const CCTK_INT *nk,
+ CCTK_REAL *mask);
+
+
+ ni = cctk_lsh[0]; nj = cctk_lsh[1]; nk = cctk_lsh[2];
+
+ CCTK_FNAME(sbp_set_norm_mask)(cctkGH_,&ni,&nj,&nk,mask);
+}
diff --git a/src/make.code.defn b/src/make.code.defn
index a5db2d8..893ee3f 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,7 +2,7 @@
# $Header$
# Source files in this directory
-SRCS = Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 Derivatives_8_4.F90 call_derivs.c call_derivs_name.c
+SRCS = Derivatives_2_1.F90 Derivatives_4_2.F90 Derivatives_6_3.F90 Derivatives_8_4.F90 call_derivs.c Derivatives_4_3.F90 call_derivs_name.c Get_Coeff.F90 call_norm_mask.c set_norm_mask.F90
# Subdirectories containing source files
SUBDIRS =
diff --git a/src/set_norm_mask.F90 b/src/set_norm_mask.F90
new file mode 100644
index 0000000..30b3c63
--- /dev/null
+++ b/src/set_norm_mask.F90
@@ -0,0 +1,124 @@
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Functions.h"
+#include "cctk_Parameters.h"
+
+subroutine sbp_set_norm_mask (cctkGH, ni, nj, nk, mask)
+
+ implicit none
+
+ DECLARE_CCTK_FUNCTIONS
+ DECLARE_CCTK_PARAMETERS
+
+ CCTK_POINTER, intent(IN) :: cctkGH
+ CCTK_INT, intent(IN) :: ni, nj, nk
+ CCTK_REAL, dimension(ni,nj,nk), intent(OUT) :: mask
+
+ CCTK_REAL, parameter :: zero = 0.0
+ integer, parameter :: wp = kind(zero)
+ CCTK_INT :: symtable, n_elements, nchar, pen_sym_handle, np
+ CCTK_INT, dimension(6) :: symbnd
+ CCTK_POINTER :: psym_name
+ character(len=256) :: symmetry_name
+ CCTK_REAL, dimension(:), allocatable :: mask_x, mask_y, mask_z
+
+! Note: The first number is twice the value from the paper, since Carpet
+! Multiplies by 1/2 at the boundary when doing the sum reduction.
+ CCTK_REAL, dimension(2), parameter :: bmask_2 = (/ 1.0_wp, 1.0_wp /)
+ CCTK_REAL, dimension(3), parameter :: bmask_3 = (/ 1.0_wp, 1.0_wp, 1.0_wp /)
+ CCTK_REAL, dimension(4), parameter :: bmask_4 = (/ 34.0_wp/48.0_wp, &
+ 59.0_wp/48.0_wp, &
+ 43.0_wp/48.0_wp, &
+ 49.0_wp/48.0_wp /)
+ CCTK_REAL, dimension(6), parameter :: bmask_6 = (/ 27298.0_wp/43200._wp, &
+ 12013.0_wp/8640._wp, &
+ 2711.0_wp/4320.0_wp, &
+ 5359.0_wp/4320.0_wp, &
+ 7877.0_wp/8640.0_wp, &
+ 43801.0_wp/43200.0_wp /)
+ CCTK_REAL, dimension(8), parameter :: bmask_8 = (/ 2996278.0_wp/5080320.0_wp,&
+ 1107307.0_wp/725760.0_wp, &
+ 20761.0_wp/80640.0_wp, &
+ 1304999.0_wp/725760.0_wp, &
+ 299527.0_wp/725760.0_wp, &
+ 103097.0_wp/80640.0_wp, &
+ 670091.0_wp/725760.0_wp, &
+ 5127739.0_wp/5080320.0_wp/)
+
+ CCTK_REAL, dimension(8) :: bmask
+
+ symtable = SymmetryTableHandleForGrid ( cctkGH )
+ call Util_TableGetIntArray ( n_elements, symtable, 6, &
+ symbnd, 'symmetry_handle' )
+
+ pen_sym_handle = SymmetryHandleOfName ( 'penalty inter-patch' )
+
+ if ( any ( symbnd == pen_sym_handle ) ) then
+ allocate ( mask_x(ni), mask_y(nj), mask_z(nk) )
+ mask_x = 1.0d0; mask_y = 1.0d0; mask_z = 1.0d0
+
+ if ( CCTK_EQUALS(norm_type,'Diagonal') ) then
+ select case (order)
+ case (2)
+ bmask(1:2) = bmask_2
+ case (4)
+ bmask(1:4) = bmask_4
+ case(6)
+ bmask(1:6) = bmask_6
+ case(8)
+ bmask(1:8) = bmask_8
+ end select
+ np = order
+ else
+ select case (order)
+ case (4)
+ bmask(1:3) = bmask_3
+ np = 3
+ end select
+ end if
+
+ if ( symbnd(1) == pen_sym_handle ) then
+ mask_x(1:np) = bmask(1:np)
+ end if
+ if ( symbnd(2) == pen_sym_handle ) then
+ mask_x(ni:ni-np+1:-1) = bmask(1:np)
+ end if
+ if ( symbnd(3) == pen_sym_handle ) then
+ mask_y(1:np) = bmask(1:np)
+ end if
+ if ( symbnd(4) == pen_sym_handle ) then
+ mask_y(nj:nj-np+1:-1) = bmask(1:np)
+ end if
+ if ( symbnd(5) == pen_sym_handle ) then
+ mask_z(1:np) = bmask(1:np)
+ end if
+ if ( symbnd(6) == pen_sym_handle ) then
+ mask_z(nk:nk-np+1:-1) = bmask(1:np)
+ end if
+
+ mask = outer_prod_3d ( mask_x, mask_y, mask_z )
+ else
+ mask = 1.0_wp
+ end if
+
+contains
+
+ function outer_prod_3d ( ax, ay, az )
+ CCTK_REAL, dimension(:), intent(IN) :: ax
+ CCTK_REAL, dimension(:), intent(IN) :: ay
+ CCTK_REAL, dimension(:), intent(IN) :: az
+ CCTK_REAL, dimension(size(ax),size(ay),size(az)) :: outer_prod_3d
+ CCTK_INT :: nx, ny, nz, i, j, k
+
+ nx = size(ax); ny = size(ay); nz = size(az)
+ do k = 1, nz
+ do j = 1, ny
+ do i = 1, nx
+ outer_prod_3d(i,j,k) = ax(i) * ay(j) * az(k)
+ end do
+ end do
+ end do
+ end function outer_prod_3d
+
+end subroutine sbp_set_norm_mask
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-12 12:20:08 +0000