diff options
24 files changed, 1129 insertions, 806 deletions
diff --git a/ML_ADM/param.ccl b/ML_ADM/param.ccl index 6ca561f..df67e04 100644 --- a/ML_ADM/param.ccl +++ b/ML_ADM/param.ccl @@ -32,11 +32,11 @@ KEYWORD my_initial_data "my_initial_data" } "ADMBase" private: -KEYWORD SpaceTime "SpaceTime" +KEYWORD my_boundary_condition "my_boundary_condition" { - "Space" :: "Space" - "Space+Matter" :: "Space+Matter" -} "Space" + "none" :: "none" + "Minkowski" :: "Minkowski" +} "none" restricted: CCTK_INT ML_ADM_MaxNumEvolvedVars "Number of evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars @@ -50,6 +50,12 @@ CCTK_INT ML_ADM_MaxNumConstrainedVars "Number of constrained variables used by t 38:38 :: "Number of constrained variables used by this thorn" } 38 +private: +CCTK_INT timelevels "Number of active timelevels" +{ + 0:3 :: "" +} 3 + restricted: CCTK_INT ML_ADM_Minkowski_calc_every "ML_ADM_Minkowski_calc_every" { @@ -69,6 +75,12 @@ CCTK_INT ML_ADM_RHS_calc_every "ML_ADM_RHS_calc_every" } 1 restricted: +CCTK_INT ML_ADM_boundary_calc_every "ML_ADM_boundary_calc_every" +{ + *:* :: "" +} 1 + +restricted: CCTK_INT ML_ADM_convertToADMBase_calc_every "ML_ADM_convertToADMBase_calc_every" { *:* :: "" @@ -81,6 +93,12 @@ CCTK_INT ML_ADM_constraints_calc_every "ML_ADM_constraints_calc_every" } 1 restricted: +CCTK_INT ML_ADM_constraints_boundary_calc_every "ML_ADM_constraints_boundary_calc_every" +{ + *:* :: "" +} 1 + +restricted: CCTK_INT ML_ADM_Minkowski_calc_offset "ML_ADM_Minkowski_calc_offset" { *:* :: "" @@ -99,6 +117,12 @@ CCTK_INT ML_ADM_RHS_calc_offset "ML_ADM_RHS_calc_offset" } 0 restricted: +CCTK_INT ML_ADM_boundary_calc_offset "ML_ADM_boundary_calc_offset" +{ + *:* :: "" +} 0 + +restricted: CCTK_INT ML_ADM_convertToADMBase_calc_offset "ML_ADM_convertToADMBase_calc_offset" { *:* :: "" @@ -110,6 +134,12 @@ CCTK_INT ML_ADM_constraints_calc_offset "ML_ADM_constraints_calc_offset" *:* :: "" } 0 +restricted: +CCTK_INT ML_ADM_constraints_boundary_calc_offset "ML_ADM_constraints_boundary_calc_offset" +{ + *:* :: "" +} 0 + private: KEYWORD K11_bound "Boundary condition to implement" { diff --git a/ML_ADM/schedule.ccl b/ML_ADM/schedule.ccl index 36a4cdd..5dbca42 100644 --- a/ML_ADM/schedule.ccl +++ b/ML_ADM/schedule.ccl @@ -16,13 +16,57 @@ STORAGE: ml_metricrhs[1] STORAGE: ml_shiftrhs[1] -STORAGE: ml_curv[3] +if (timelevels == 1) +{ + STORAGE: ml_curv[1] +} +if (timelevels == 2) +{ + STORAGE: ml_curv[2] +} +if (timelevels == 3) +{ + STORAGE: ml_curv[3] +} -STORAGE: ml_lapse[3] +if (timelevels == 1) +{ + STORAGE: ml_lapse[1] +} +if (timelevels == 2) +{ + STORAGE: ml_lapse[2] +} +if (timelevels == 3) +{ + STORAGE: ml_lapse[3] +} -STORAGE: ml_metric[3] +if (timelevels == 1) +{ + STORAGE: ml_metric[1] +} +if (timelevels == 2) +{ + STORAGE: ml_metric[2] +} +if (timelevels == 3) +{ + STORAGE: ml_metric[3] +} -STORAGE: ml_shift[3] +if (timelevels == 1) +{ + STORAGE: ml_shift[1] +} +if (timelevels == 2) +{ + STORAGE: ml_shift[2] +} +if (timelevels == 3) +{ + STORAGE: ml_shift[3] +} schedule ML_ADM_Startup at STARTUP { @@ -48,7 +92,6 @@ if (CCTK_EQUALS(my_initial_data, "Minkowski")) schedule ML_ADM_Minkowski IN ADMBase_InitialData { LANG: C - } "ML_ADM_Minkowski" } @@ -58,40 +101,51 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase")) schedule ML_ADM_convertFromADMBase AT initial AFTER ADMBase_PostInitial { LANG: C - } "ML_ADM_convertFromADMBase" } schedule ML_ADM_RHS IN MoL_CalcRHS { LANG: C - } "ML_ADM_RHS" schedule ML_ADM_RHS AT analysis { LANG: C - SYNC: ml_curvrhs SYNC: ml_lapserhs SYNC: ml_metricrhs SYNC: ml_shiftrhs } "ML_ADM_RHS" -schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER ML_ADM_ApplyBCs + +if (CCTK_EQUALS(my_boundary_condition, "Minkowski")) +{ + schedule ML_ADM_boundary IN MoL_PostStep + { + LANG: C + } "ML_ADM_boundary" +} + +schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER (ML_ADM_ApplyBCs ML_ADM_boundary) { LANG: C - } "ML_ADM_convertToADMBase" schedule ML_ADM_constraints AT analysis { LANG: C - SYNC: Ham SYNC: mom + TRIGGERS: Ham + TRIGGERS: mom } "ML_ADM_constraints" +schedule ML_ADM_constraints_boundary AT analysis AFTER ML_ADM_constraints +{ + LANG: C +} "ML_ADM_constraints_boundary" + schedule ML_ADM_ApplyBoundConds in MoL_PostStep { LANG: C @@ -111,5 +165,4 @@ schedule ML_ADM_CheckBoundaries at BASEGRID schedule group ApplyBCs as ML_ADM_ApplyBCs in MoL_PostStep after ML_ADM_ApplyBoundConds { # no language specified - } "Apply boundary conditions controlled by thorn Boundary" diff --git a/ML_ADM/src/ML_ADM_boundary.c b/ML_ADM/src/ML_ADM_boundary.c new file mode 100644 index 0000000..02c77aa --- /dev/null +++ b/ML_ADM/src/ML_ADM_boundary.c @@ -0,0 +1,181 @@ +/* File produced by user eschnett */ +/* Produced with Mathematica Version 6.0 for Mac OS X x86 (32-bit) (April 20, 2007) */ + +/* Mathematica script written by Ian Hinder and Sascha Husa */ + +#define KRANC_C + +#include <math.h> +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "GenericFD.h" +#include "Differencing.h" +#include "loopcontrol.h" + +/* Define macros used in calculations */ +#define INITVALUE (42) +#define INV(x) ((1.0) / (x)) +#define SQR(x) ((x) * (x)) +#define CUB(x) ((x) * (x) * (x)) +#define QAD(x) ((x) * (x) * (x) * (x)) + +void ML_ADM_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[]) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + + /* Declare finite differencing variables */ + CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; + CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; + CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; + CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; + + + /* Declare predefined quantities */ + CCTK_REAL p1o12dx = INITVALUE; + CCTK_REAL p1o12dy = INITVALUE; + CCTK_REAL p1o12dz = INITVALUE; + CCTK_REAL p1o144dxdy = INITVALUE; + CCTK_REAL p1o144dxdz = INITVALUE; + CCTK_REAL p1o144dydz = INITVALUE; + CCTK_REAL pm1o12dx2 = INITVALUE; + CCTK_REAL pm1o12dy2 = INITVALUE; + CCTK_REAL pm1o12dz2 = INITVALUE; + + if (verbose > 1) + { + CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_boundary_Body"); + } + + if (cctk_iteration % ML_ADM_boundary_calc_every != ML_ADM_boundary_calc_offset) + { + return; + } + + /* Include user-supplied include files */ + + /* Initialise finite differencing variables */ + dx = CCTK_DELTA_SPACE(0); + dy = CCTK_DELTA_SPACE(1); + dz = CCTK_DELTA_SPACE(2); + dxi = 1.0 / dx; + dyi = 1.0 / dy; + dzi = 1.0 / dz; + khalf = 0.5; + kthird = 1/3.0; + ktwothird = 2.0/3.0; + kfourthird = 4.0/3.0; + keightthird = 8.0/3.0; + hdxi = 0.5 * dxi; + hdyi = 0.5 * dyi; + hdzi = 0.5 * dzi; + + /* Initialize predefined quantities */ + p1o12dx = INV(dx)/12.; + p1o12dy = INV(dy)/12.; + p1o12dz = INV(dz)/12.; + p1o144dxdy = (INV(dx)*INV(dy))/144.; + p1o144dxdz = (INV(dx)*INV(dz))/144.; + p1o144dydz = (INV(dy)*INV(dz))/144.; + pm1o12dx2 = -pow(dx,-2)/12.; + pm1o12dy2 = -pow(dy,-2)/12.; + pm1o12dz2 = -pow(dz,-2)/12.; + + /* Loop over the grid points */ + _Pragma ("omp parallel") + LC_LOOP3 (ML_ADM_boundary, + i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int index = INITVALUE; + int subblock_index = INITVALUE; + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + + /* Declare shorthands */ + + /* Declare local copies of grid functions */ + CCTK_REAL alphaL = INITVALUE; + CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; + CCTK_REAL g11L = INITVALUE, g12L = INITVALUE, g13L = INITVALUE, g22L = INITVALUE, g23L = INITVALUE, g33L = INITVALUE; + CCTK_REAL K11L = INITVALUE, K12L = INITVALUE, K13L = INITVALUE, K22L = INITVALUE, K23L = INITVALUE, K33L = INITVALUE; + /* Declare precomputed derivatives*/ + + /* Declare derivatives */ + + /* Assign local copies of grid functions */ + + /* Assign local copies of subblock grid functions */ + + /* Include user supplied include files */ + + /* Precompute derivatives (new style) */ + + /* Precompute derivatives (old style) */ + + /* Calculate temporaries and grid functions */ + g11L = 1; + + g12L = 0; + + g13L = 0; + + g22L = 1; + + g23L = 0; + + g33L = 1; + + K11L = 0; + + K12L = 0; + + K13L = 0; + + K22L = 0; + + K23L = 0; + + K33L = 0; + + alphaL = 1; + + beta1L = 0; + + beta2L = 0; + + beta3L = 0; + + + /* Copy local copies back to grid functions */ + alpha[index] = alphaL; + beta1[index] = beta1L; + beta2[index] = beta2L; + beta3[index] = beta3L; + g11[index] = g11L; + g12[index] = g12L; + g13[index] = g13L; + g22[index] = g22L; + g23[index] = g23L; + g33[index] = g33L; + K11[index] = K11L; + K12[index] = K12L; + K13[index] = K13L; + K22[index] = K22L; + K23[index] = K23L; + K33[index] = K33L; + + /* Copy local copies back to subblock grid functions */ + } + LC_ENDLOOP3 (ML_ADM_boundary); +} + +void ML_ADM_boundary(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + GenericFD_LoopOverBoundary(cctkGH, &ML_ADM_boundary_Body); +} diff --git a/ML_ADM/src/ML_ADM_constraints_boundary.c b/ML_ADM/src/ML_ADM_constraints_boundary.c new file mode 100644 index 0000000..7a7053d --- /dev/null +++ b/ML_ADM/src/ML_ADM_constraints_boundary.c @@ -0,0 +1,143 @@ +/* File produced by user eschnett */ +/* Produced with Mathematica Version 6.0 for Mac OS X x86 (32-bit) (April 20, 2007) */ + +/* Mathematica script written by Ian Hinder and Sascha Husa */ + +#define KRANC_C + +#include <math.h> +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "GenericFD.h" +#include "Differencing.h" +#include "loopcontrol.h" + +/* Define macros used in calculations */ +#define INITVALUE (42) +#define INV(x) ((1.0) / (x)) +#define SQR(x) ((x) * (x)) +#define CUB(x) ((x) * (x) * (x)) +#define QAD(x) ((x) * (x) * (x) * (x)) + +void ML_ADM_constraints_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[]) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + + /* Declare finite differencing variables */ + CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; + CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; + CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; + CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; + + + /* Declare predefined quantities */ + CCTK_REAL p1o12dx = INITVALUE; + CCTK_REAL p1o12dy = INITVALUE; + CCTK_REAL p1o12dz = INITVALUE; + CCTK_REAL p1o144dxdy = INITVALUE; + CCTK_REAL p1o144dxdz = INITVALUE; + CCTK_REAL p1o144dydz = INITVALUE; + CCTK_REAL pm1o12dx2 = INITVALUE; + CCTK_REAL pm1o12dy2 = INITVALUE; + CCTK_REAL pm1o12dz2 = INITVALUE; + + if (verbose > 1) + { + CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_ADM_constraints_boundary_Body"); + } + + if (cctk_iteration % ML_ADM_constraints_boundary_calc_every != ML_ADM_constraints_boundary_calc_offset) + { + return; + } + + /* Include user-supplied include files */ + + /* Initialise finite differencing variables */ + dx = CCTK_DELTA_SPACE(0); + dy = CCTK_DELTA_SPACE(1); + dz = CCTK_DELTA_SPACE(2); + dxi = 1.0 / dx; + dyi = 1.0 / dy; + dzi = 1.0 / dz; + khalf = 0.5; + kthird = 1/3.0; + ktwothird = 2.0/3.0; + kfourthird = 4.0/3.0; + keightthird = 8.0/3.0; + hdxi = 0.5 * dxi; + hdyi = 0.5 * dyi; + hdzi = 0.5 * dzi; + + /* Initialize predefined quantities */ + p1o12dx = INV(dx)/12.; + p1o12dy = INV(dy)/12.; + p1o12dz = INV(dz)/12.; + p1o144dxdy = (INV(dx)*INV(dy))/144.; + p1o144dxdz = (INV(dx)*INV(dz))/144.; + p1o144dydz = (INV(dy)*INV(dz))/144.; + pm1o12dx2 = -pow(dx,-2)/12.; + pm1o12dy2 = -pow(dy,-2)/12.; + pm1o12dz2 = -pow(dz,-2)/12.; + + /* Loop over the grid points */ + _Pragma ("omp parallel") + LC_LOOP3 (ML_ADM_constraints_boundary, + i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int index = INITVALUE; + int subblock_index = INITVALUE; + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + + /* Declare shorthands */ + + /* Declare local copies of grid functions */ + CCTK_REAL HL = INITVALUE; + CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; + /* Declare precomputed derivatives*/ + + /* Declare derivatives */ + + /* Assign local copies of grid functions */ + + /* Assign local copies of subblock grid functions */ + + /* Include user supplied include files */ + + /* Precompute derivatives (new style) */ + + /* Precompute derivatives (old style) */ + + /* Calculate temporaries and grid functions */ + HL = 0; + + M1L = 0; + + M2L = 0; + + M3L = 0; + + + /* Copy local copies back to grid functions */ + H[index] = HL; + M1[index] = M1L; + M2[index] = M2L; + M3[index] = M3L; + + /* Copy local copies back to subblock grid functions */ + } + LC_ENDLOOP3 (ML_ADM_constraints_boundary); +} + +void ML_ADM_constraints_boundary(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + GenericFD_LoopOverBoundary(cctkGH, &ML_ADM_constraints_boundary_Body); +} diff --git a/ML_ADM/src/make.code.defn b/ML_ADM/src/make.code.defn index 6e2a3f8..9762fc1 100644 --- a/ML_ADM/src/make.code.defn +++ b/ML_ADM/src/make.code.defn @@ -3,4 +3,4 @@ # Mathematica script written by Ian Hinder and Sascha Husa -SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_ADM_Minkowski.c ML_ADM_convertFromADMBase.c ML_ADM_RHS.c ML_ADM_convertToADMBase.c ML_ADM_constraints.c Boundaries.c +SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_ADM_Minkowski.c ML_ADM_convertFromADMBase.c ML_ADM_RHS.c ML_ADM_boundary.c ML_ADM_convertToADMBase.c ML_ADM_constraints.c ML_ADM_constraints_boundary.c Boundaries.c diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl index bd69853..af97786 100644 --- a/ML_BSSN/param.ccl +++ b/ML_BSSN/param.ccl @@ -7,36 +7,6 @@ shares: ADMBase -EXTENDS CCTK_KEYWORD initial_data "initial_data" -{ - ML_BSSN__Minkowski :: "" -} - - -EXTENDS CCTK_KEYWORD initial_lapse "initial_lapse" -{ - ML_BSSN__Minkowski :: "" -} - - -EXTENDS CCTK_KEYWORD initial_shift "initial_shift" -{ - ML_BSSN__Minkowski :: "" -} - - -EXTENDS CCTK_KEYWORD initial_dtlapse "initial_dtlapse" -{ - ML_BSSN__Minkowski :: "" -} - - -EXTENDS CCTK_KEYWORD initial_dtshift "initial_dtshift" -{ - ML_BSSN__Minkowski :: "" -} - - EXTENDS CCTK_KEYWORD evolution_method "evolution_method" { ML_BSSN :: "" @@ -104,13 +74,13 @@ restricted: CCTK_REAL LapseAdvectionCoeff "Factor in front of the shift advection terms in 1+log" { "*:*" :: "" -} 1. +} 1 restricted: CCTK_REAL ShiftAdvectionCoeff "Factor in front of the shift advection terms in gamma driver" { "*:*" :: "" -} 1. +} 1 restricted: CCTK_INT harmonicN "d/dt alpha = - f alpha^n K (harmonic=2, 1+log=1)" @@ -132,11 +102,11 @@ KEYWORD my_initial_data "my_initial_data" } "ADMBase" private: -KEYWORD SpaceTime "SpaceTime" +KEYWORD my_boundary_condition "my_boundary_condition" { - "Space" :: "Space" - "Space+Matter" :: "Space+Matter" -} "Space" + "none" :: "none" + "Minkowski" :: "Minkowski" +} "none" restricted: CCTK_INT ML_BSSN_MaxNumEvolvedVars "Number of evolved variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Evolved_Vars @@ -150,6 +120,12 @@ CCTK_INT ML_BSSN_MaxNumConstrainedVars "Number of constrained variables used by 43:43 :: "Number of constrained variables used by this thorn" } 43 +private: +CCTK_INT timelevels "Number of active timelevels" +{ + 0:3 :: "" +} 3 + restricted: CCTK_INT ML_BSSN_Minkowski_calc_every "ML_BSSN_Minkowski_calc_every" { @@ -175,13 +151,13 @@ CCTK_INT ML_BSSN_RHS_calc_every "ML_BSSN_RHS_calc_every" } 1 restricted: -CCTK_INT ML_BSSN_matter_calc_every "ML_BSSN_matter_calc_every" +CCTK_INT ML_BSSN_enforce_calc_every "ML_BSSN_enforce_calc_every" { *:* :: "" } 1 restricted: -CCTK_INT ML_BSSN_enforce_calc_every "ML_BSSN_enforce_calc_every" +CCTK_INT ML_BSSN_boundary_calc_every "ML_BSSN_boundary_calc_every" { *:* :: "" } 1 @@ -199,7 +175,7 @@ CCTK_INT ML_BSSN_constraints_calc_every "ML_BSSN_constraints_calc_every" } 1 restricted: -CCTK_INT ML_BSSN_matter_constraints_calc_every "ML_BSSN_matter_constraints_calc_every" +CCTK_INT ML_BSSN_constraints_boundary_calc_every "ML_BSSN_constraints_boundary_calc_every" { *:* :: "" } 1 @@ -229,13 +205,13 @@ CCTK_INT ML_BSSN_RHS_calc_offset "ML_BSSN_RHS_calc_offset" } 0 restricted: -CCTK_INT ML_BSSN_matter_calc_offset "ML_BSSN_matter_calc_offset" +CCTK_INT ML_BSSN_enforce_calc_offset "ML_BSSN_enforce_calc_offset" { *:* :: "" } 0 restricted: -CCTK_INT ML_BSSN_enforce_calc_offset "ML_BSSN_enforce_calc_offset" +CCTK_INT ML_BSSN_boundary_calc_offset "ML_BSSN_boundary_calc_offset" { *:* :: "" } 0 @@ -253,7 +229,7 @@ CCTK_INT ML_BSSN_constraints_calc_offset "ML_BSSN_constraints_calc_offset" } 0 restricted: -CCTK_INT ML_BSSN_matter_constraints_calc_offset "ML_BSSN_matter_constraints_calc_offset" +CCTK_INT ML_BSSN_constraints_boundary_calc_offset "ML_BSSN_constraints_boundary_calc_offset" { *:* :: "" } 0 diff --git a/ML_BSSN/schedule.ccl b/ML_BSSN/schedule.ccl index a000a82..1b62db3 100644 --- a/ML_BSSN/schedule.ccl +++ b/ML_BSSN/schedule.ccl @@ -32,23 +32,122 @@ STORAGE: ML_shiftrhs[1] STORAGE: ML_trace_curvrhs[1] -STORAGE: ML_curv[3] +if (timelevels == 1) +{ + STORAGE: ML_curv[1] +} +if (timelevels == 2) +{ + STORAGE: ML_curv[2] +} +if (timelevels == 3) +{ + STORAGE: ML_curv[3] +} -STORAGE: ML_dtlapse[3] +if (timelevels == 1) +{ + STORAGE: ML_dtlapse[1] +} +if (timelevels == 2) +{ + STORAGE: ML_dtlapse[2] +} +if (timelevels == 3) +{ + STORAGE: ML_dtlapse[3] +} -STORAGE: ML_dtshift[3] +if (timelevels == 1) +{ + STORAGE: ML_dtshift[1] +} +if (timelevels == 2) +{ + STORAGE: ML_dtshift[2] +} +if (timelevels == 3) +{ + STORAGE: ML_dtshift[3] +} -STORAGE: ML_Gamma[3] +if (timelevels == 1) +{ + STORAGE: ML_Gamma[1] +} +if (timelevels == 2) +{ + STORAGE: ML_Gamma[2] +} +if (timelevels == 3) +{ + STORAGE: ML_Gamma[3] +} -STORAGE: ML_lapse[3] +if (timelevels == 1) +{ + STORAGE: ML_lapse[1] +} +if (timelevels == 2) +{ + STORAGE: ML_lapse[2] +} +if (timelevels == 3) +{ + STORAGE: ML_lapse[3] +} -STORAGE: ML_log_confac[3] +if (timelevels == 1) +{ + STORAGE: ML_log_confac[1] +} +if (timelevels == 2) +{ + STORAGE: ML_log_confac[2] +} +if (timelevels == 3) +{ + STORAGE: ML_log_confac[3] +} -STORAGE: ML_metric[3] +if (timelevels == 1) +{ + STORAGE: ML_metric[1] +} +if (timelevels == 2) +{ + STORAGE: ML_metric[2] +} +if (timelevels == 3) +{ + STORAGE: ML_metric[3] +} -STORAGE: ML_shift[3] +if (timelevels == 1) +{ + STORAGE: ML_shift[1] +} +if (timelevels == 2) +{ + STORAGE: ML_shift[2] +} +if (timelevels == 3) +{ + STORAGE: ML_shift[3] +} -STORAGE: ML_trace_curv[3] +if (timelevels == 1) +{ + STORAGE: ML_trace_curv[1] +} +if (timelevels == 2) +{ + STORAGE: ML_trace_curv[2] +} +if (timelevels == 3) +{ + STORAGE: ML_trace_curv[3] +} schedule ML_BSSN_Startup at STARTUP { @@ -69,12 +168,11 @@ schedule ML_BSSN_RegisterSymmetries at BASEGRID } "register symmetries" -if (CCTK_EQUALS(initial_data, "ML_BSSN__Minkowski")) +if (CCTK_EQUALS(my_initial_data, "Minkowski")) { schedule ML_BSSN_Minkowski IN ADMBase_InitialData { LANG: C - } "ML_BSSN_Minkowski" } @@ -84,7 +182,6 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase")) schedule ML_BSSN_convertFromADMBase AT initial AFTER ADMBase_PostInitial { LANG: C - } "ML_BSSN_convertFromADMBase" } @@ -94,87 +191,50 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase")) schedule ML_BSSN_convertFromADMBaseGamma AT initial AFTER ML_BSSN_convertFromADMBase { LANG: C - SYNC: ML_dtlapse SYNC: ML_dtshift SYNC: ML_Gamma } "ML_BSSN_convertFromADMBaseGamma" } - -if (CCTK_EQUALS(evolution_method, "ML_BSSN")) +schedule ML_BSSN_RHS IN ML_BSSN_evolCalcGroup { - schedule ML_BSSN_RHS IN MoL_CalcRHS - { - LANG: C - - TRIGGERS: ML_log_confacrhs - TRIGGERS: ML_metricrhs - TRIGGERS: ML_Gammarhs - TRIGGERS: ML_trace_curvrhs - TRIGGERS: ML_curvrhs - TRIGGERS: ML_lapserhs - TRIGGERS: ML_dtlapserhs - TRIGGERS: ML_shiftrhs - TRIGGERS: ML_dtshiftrhs - } "ML_BSSN_RHS" -} + LANG: C + SYNC: ML_curvrhs + SYNC: ML_dtlapserhs + SYNC: ML_dtshiftrhs + SYNC: ML_Gammarhs + SYNC: ML_lapserhs + SYNC: ML_log_confacrhs + SYNC: ML_metricrhs + SYNC: ML_shiftrhs + SYNC: ML_trace_curvrhs +} "ML_BSSN_RHS" if (CCTK_EQUALS(evolution_method, "ML_BSSN")) { - schedule ML_BSSN_RHS AT analysis - { - LANG: C - - SYNC: ML_curvrhs - SYNC: ML_dtlapserhs - SYNC: ML_dtshiftrhs - SYNC: ML_Gammarhs - SYNC: ML_lapserhs - SYNC: ML_log_confacrhs - SYNC: ML_metricrhs - SYNC: ML_shiftrhs - SYNC: ML_trace_curvrhs - TRIGGERS: ML_log_confacrhs - TRIGGERS: ML_metricrhs - TRIGGERS: ML_Gammarhs - TRIGGERS: ML_trace_curvrhs - TRIGGERS: ML_curvrhs - TRIGGERS: ML_lapserhs - TRIGGERS: ML_dtlapserhs - TRIGGERS: ML_shiftrhs - TRIGGERS: ML_dtshiftrhs - } "ML_BSSN_RHS" -} - - -if (CCTK_EQUALS(SpaceTime, "Space+Matter")) -{ - schedule ML_BSSN_matter IN MoL_CalcRHS AFTER ML_BSSN_RHS + schedule ML_BSSN_enforce IN MoL_PostStep BEFORE ML_BSSN_BoundConds { LANG: C - - } "ML_BSSN_matter" + } "ML_BSSN_enforce" } -if (CCTK_EQUALS(evolution_method, "ML_BSSN")) +if (CCTK_EQUALS(my_boundary_condition, "Minkowski")) { - schedule ML_BSSN_enforce IN MoL_PostStep BEFORE ML_BSSN_BoundConds + schedule ML_BSSN_boundary IN MoL_PostStep { LANG: C - - } "ML_BSSN_enforce" + } "ML_BSSN_boundary" } if (CCTK_EQUALS(evolution_method, "ML_BSSN")) { - schedule ML_BSSN_convertToADMBase IN MoL_PostStep AFTER ML_BSSN_ApplyBCs AFTER ML_BSSN_enforce + schedule ML_BSSN_convertToADMBase IN MoL_PostStep AFTER (ML_BSSN_ApplyBCs ML_BSSN_boundary ML_BSSN_enforce) { LANG: C - SYNC: ADMBase::curv SYNC: ADMBase::dtlapse SYNC: ADMBase::dtshift @@ -184,36 +244,20 @@ if (CCTK_EQUALS(evolution_method, "ML_BSSN")) } "ML_BSSN_convertToADMBase" } - -if (CCTK_EQUALS(evolution_method, "ML_BSSN")) +schedule ML_BSSN_constraints IN ML_BSSN_constraintsCalcGroup { - schedule ML_BSSN_constraints AT analysis - { - LANG: C - - SYNC: cons_detg - SYNC: cons_Gamma - SYNC: cons_traceA - SYNC: Ham - SYNC: mom - TRIGGERS: Ham - TRIGGERS: mom - } "ML_BSSN_constraints" -} - - -if (CCTK_EQUALS(SpaceTime, "Space+Matter")) + LANG: C + SYNC: cons_detg + SYNC: cons_Gamma + SYNC: cons_traceA + SYNC: Ham + SYNC: mom +} "ML_BSSN_constraints" + +schedule ML_BSSN_constraints_boundary IN ML_BSSN_constraintsCalcGroup AFTER ML_BSSN_constraints { - schedule ML_BSSN_matter_constraints AT analysis AFTER ML_BSSN_constraints - { - LANG: C - - SYNC: Ham - SYNC: mom - TRIGGERS: Ham - TRIGGERS: mom - } "ML_BSSN_matter_constraints" -} + LANG: C +} "ML_BSSN_constraints_boundary" schedule ML_BSSN_ApplyBoundConds in MoL_PostStep { @@ -239,5 +283,4 @@ schedule ML_BSSN_CheckBoundaries at BASEGRID schedule group ApplyBCs as ML_BSSN_ApplyBCs in MoL_PostStep after ML_BSSN_ApplyBoundConds { # no language specified - } "Apply boundary conditions controlled by thorn Boundary" diff --git a/ML_BSSN/src/ML_BSSN_matter_constraints.c b/ML_BSSN/src/ML_BSSN_boundary.c index b467210..41f2fda 100644 --- a/ML_BSSN/src/ML_BSSN_matter_constraints.c +++ b/ML_BSSN/src/ML_BSSN_boundary.c @@ -20,7 +20,7 @@ #define CUB(x) ((x) * (x) * (x)) #define QAD(x) ((x) * (x) * (x) * (x)) -void ML_BSSN_matter_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[]) +void ML_BSSN_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[]) { DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS @@ -46,10 +46,10 @@ void ML_BSSN_matter_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C if (verbose > 1) { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_matter_constraints_Body"); + CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_boundary_Body"); } - if (cctk_iteration % ML_BSSN_matter_constraints_calc_every != ML_BSSN_matter_constraints_calc_offset) + if (cctk_iteration % ML_BSSN_boundary_calc_every != ML_BSSN_boundary_calc_offset) { return; } @@ -85,7 +85,7 @@ void ML_BSSN_matter_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C /* Loop over the grid points */ _Pragma ("omp parallel") - LC_LOOP3 (ML_BSSN_matter_constraints, + LC_LOOP3 (ML_BSSN_boundary, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) { @@ -95,49 +95,22 @@ void ML_BSSN_matter_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); /* Declare shorthands */ - CCTK_REAL rho = INITVALUE; - CCTK_REAL S1 = INITVALUE, S2 = INITVALUE, S3 = INITVALUE; - CCTK_REAL T00 = INITVALUE, T01 = INITVALUE, T02 = INITVALUE, T03 = INITVALUE, T11 = INITVALUE, T12 = INITVALUE; - CCTK_REAL T13 = INITVALUE, T22 = INITVALUE, T23 = INITVALUE, T33 = INITVALUE; /* Declare local copies of grid functions */ + CCTK_REAL AL = INITVALUE; CCTK_REAL alphaL = INITVALUE; + CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE; + CCTK_REAL B1L = INITVALUE, B2L = INITVALUE, B3L = INITVALUE; CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL eTttL = INITVALUE; - CCTK_REAL eTtxL = INITVALUE; - CCTK_REAL eTtyL = INITVALUE; - CCTK_REAL eTtzL = INITVALUE; - CCTK_REAL eTxxL = INITVALUE; - CCTK_REAL eTxyL = INITVALUE; - CCTK_REAL eTxzL = INITVALUE; - CCTK_REAL eTyyL = INITVALUE; - CCTK_REAL eTyzL = INITVALUE; - CCTK_REAL eTzzL = INITVALUE; - CCTK_REAL HL = INITVALUE; - CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; + CCTK_REAL gt11L = INITVALUE, gt12L = INITVALUE, gt13L = INITVALUE, gt22L = INITVALUE, gt23L = INITVALUE, gt33L = INITVALUE; + CCTK_REAL phiL = INITVALUE; + CCTK_REAL trKL = INITVALUE; + CCTK_REAL Xt1L = INITVALUE, Xt2L = INITVALUE, Xt3L = INITVALUE; /* Declare precomputed derivatives*/ /* Declare derivatives */ /* Assign local copies of grid functions */ - alphaL = alpha[index]; - beta1L = beta1[index]; - beta2L = beta2[index]; - beta3L = beta3[index]; - eTttL = eTtt[index]; - eTtxL = eTtx[index]; - eTtyL = eTty[index]; - eTtzL = eTtz[index]; - eTxxL = eTxx[index]; - eTxyL = eTxy[index]; - eTxzL = eTxz[index]; - eTyyL = eTyy[index]; - eTyzL = eTyz[index]; - eTzzL = eTzz[index]; - HL = H[index]; - M1L = M1[index]; - M2L = M2[index]; - M3L = M3[index]; /* Assign local copies of subblock grid functions */ @@ -148,60 +121,93 @@ void ML_BSSN_matter_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, C /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - T00 = eTttL; + phiL = 0; - T01 = eTtxL; + gt11L = 1; - T02 = eTtyL; + gt12L = 0; - T03 = eTtzL; + gt13L = 0; - T11 = eTxxL; + gt22L = 1; - T12 = eTxyL; + gt23L = 0; - T13 = eTxzL; + gt33L = 1; - T22 = eTyyL; + trKL = 0; - T23 = eTyzL; + At11L = 0; - T33 = eTzzL; + At12L = 0; - rho = pow(alphaL,-2)*(T00 - 2*(beta2L*T02 + beta3L*T03) + - 2*(beta1L*(-T01 + beta2L*T12 + beta3L*T13) + beta2L*beta3L*T23) + T11*SQR(beta1L) + T22*SQR(beta2L) + - T33*SQR(beta3L)); + At13L = 0; - S1 = (-T01 + beta1L*T11 + beta2L*T12 + beta3L*T13)*INV(alphaL); + At22L = 0; - S2 = (-T02 + beta1L*T12 + beta2L*T22 + beta3L*T23)*INV(alphaL); + At23L = 0; - S3 = (-T03 + beta1L*T13 + beta2L*T23 + beta3L*T33)*INV(alphaL); + At33L = 0; - HL = HL - 50.26548245743669181540229413247204614715*rho; + Xt1L = 0; - M1L = M1L - 25.13274122871834590770114706623602307358*S1; + Xt2L = 0; - M2L = M2L - 25.13274122871834590770114706623602307358*S2; + Xt3L = 0; - M3L = M3L - 25.13274122871834590770114706623602307358*S3; + alphaL = 1; + + AL = 0; + + beta1L = 0; + + beta2L = 0; + + beta3L = 0; + + B1L = 0; + + B2L = 0; + + B3L = 0; /* Copy local copies back to grid functions */ - H[index] = HL; - M1[index] = M1L; - M2[index] = M2L; - M3[index] = M3L; + A[index] = AL; + alpha[index] = alphaL; + At11[index] = At11L; + At12[index] = At12L; + At13[index] = At13L; + At22[index] = At22L; + At23[index] = At23L; + At33[index] = At33L; + B1[index] = B1L; + B2[index] = B2L; + B3[index] = B3L; + beta1[index] = beta1L; + beta2[index] = beta2L; + beta3[index] = beta3L; + gt11[index] = gt11L; + gt12[index] = gt12L; + gt13[index] = gt13L; + gt22[index] = gt22L; + gt23[index] = gt23L; + gt33[index] = gt33L; + phi[index] = phiL; + trK[index] = trKL; + Xt1[index] = Xt1L; + Xt2[index] = Xt2L; + Xt3[index] = Xt3L; /* Copy local copies back to subblock grid functions */ } - LC_ENDLOOP3 (ML_BSSN_matter_constraints); + LC_ENDLOOP3 (ML_BSSN_boundary); } -void ML_BSSN_matter_constraints(CCTK_ARGUMENTS) +void ML_BSSN_boundary(CCTK_ARGUMENTS) { DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS - GenericFD_LoopOverInterior(cctkGH, &ML_BSSN_matter_constraints_Body); + GenericFD_LoopOverBoundary(cctkGH, &ML_BSSN_boundary_Body); } diff --git a/ML_BSSN/src/ML_BSSN_constraints_boundary.c b/ML_BSSN/src/ML_BSSN_constraints_boundary.c new file mode 100644 index 0000000..19b6084 --- /dev/null +++ b/ML_BSSN/src/ML_BSSN_constraints_boundary.c @@ -0,0 +1,143 @@ +/* File produced by user eschnett */ +/* Produced with Mathematica Version 6.0 for Mac OS X x86 (32-bit) (April 20, 2007) */ + +/* Mathematica script written by Ian Hinder and Sascha Husa */ + +#define KRANC_C + +#include <math.h> +#include "cctk.h" +#include "cctk_Arguments.h" +#include "cctk_Parameters.h" +#include "GenericFD.h" +#include "Differencing.h" +#include "loopcontrol.h" + +/* Define macros used in calculations */ +#define INITVALUE (42) +#define INV(x) ((1.0) / (x)) +#define SQR(x) ((x) * (x)) +#define CUB(x) ((x) * (x) * (x)) +#define QAD(x) ((x) * (x) * (x) * (x)) + +void ML_BSSN_constraints_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[]) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + + /* Declare finite differencing variables */ + CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; + CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; + CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; + CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; + + + /* Declare predefined quantities */ + CCTK_REAL p1o12dx = INITVALUE; + CCTK_REAL p1o12dy = INITVALUE; + CCTK_REAL p1o12dz = INITVALUE; + CCTK_REAL p1o144dxdy = INITVALUE; + CCTK_REAL p1o144dxdz = INITVALUE; + CCTK_REAL p1o144dydz = INITVALUE; + CCTK_REAL pm1o12dx2 = INITVALUE; + CCTK_REAL pm1o12dy2 = INITVALUE; + CCTK_REAL pm1o12dz2 = INITVALUE; + + if (verbose > 1) + { + CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_constraints_boundary_Body"); + } + + if (cctk_iteration % ML_BSSN_constraints_boundary_calc_every != ML_BSSN_constraints_boundary_calc_offset) + { + return; + } + + /* Include user-supplied include files */ + + /* Initialise finite differencing variables */ + dx = CCTK_DELTA_SPACE(0); + dy = CCTK_DELTA_SPACE(1); + dz = CCTK_DELTA_SPACE(2); + dxi = 1.0 / dx; + dyi = 1.0 / dy; + dzi = 1.0 / dz; + khalf = 0.5; + kthird = 1/3.0; + ktwothird = 2.0/3.0; + kfourthird = 4.0/3.0; + keightthird = 8.0/3.0; + hdxi = 0.5 * dxi; + hdyi = 0.5 * dyi; + hdzi = 0.5 * dzi; + + /* Initialize predefined quantities */ + p1o12dx = INV(dx)/12.; + p1o12dy = INV(dy)/12.; + p1o12dz = INV(dz)/12.; + p1o144dxdy = (INV(dx)*INV(dy))/144.; + p1o144dxdz = (INV(dx)*INV(dz))/144.; + p1o144dydz = (INV(dy)*INV(dz))/144.; + pm1o12dx2 = -pow(dx,-2)/12.; + pm1o12dy2 = -pow(dy,-2)/12.; + pm1o12dz2 = -pow(dz,-2)/12.; + + /* Loop over the grid points */ + _Pragma ("omp parallel") + LC_LOOP3 (ML_BSSN_constraints_boundary, + i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], + cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) + { + int index = INITVALUE; + int subblock_index = INITVALUE; + index = CCTK_GFINDEX3D(cctkGH,i,j,k); + subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); + + /* Declare shorthands */ + + /* Declare local copies of grid functions */ + CCTK_REAL HL = INITVALUE; + CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; + /* Declare precomputed derivatives*/ + + /* Declare derivatives */ + + /* Assign local copies of grid functions */ + + /* Assign local copies of subblock grid functions */ + + /* Include user supplied include files */ + + /* Precompute derivatives (new style) */ + + /* Precompute derivatives (old style) */ + + /* Calculate temporaries and grid functions */ + HL = 0; + + M1L = 0; + + M2L = 0; + + M3L = 0; + + + /* Copy local copies back to grid functions */ + H[index] = HL; + M1[index] = M1L; + M2[index] = M2L; + M3[index] = M3L; + + /* Copy local copies back to subblock grid functions */ + } + LC_ENDLOOP3 (ML_BSSN_constraints_boundary); +} + +void ML_BSSN_constraints_boundary(CCTK_ARGUMENTS) +{ + DECLARE_CCTK_ARGUMENTS + DECLARE_CCTK_PARAMETERS + + GenericFD_LoopOverBoundary(cctkGH, &ML_BSSN_constraints_boundary_Body); +} diff --git a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c index 861ee58..576deb1 100644 --- a/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c +++ b/ML_BSSN/src/ML_BSSN_convertFromADMBaseGamma.c @@ -115,15 +115,15 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa /* Declare precomputed derivatives*/ /* Declare derivatives */ - CCTK_REAL PDstandardNth1B1 = INITVALUE; - CCTK_REAL PDstandardNth2B1 = INITVALUE; - CCTK_REAL PDstandardNth3B1 = INITVALUE; - CCTK_REAL PDstandardNth1B2 = INITVALUE; - CCTK_REAL PDstandardNth2B2 = INITVALUE; - CCTK_REAL PDstandardNth3B2 = INITVALUE; - CCTK_REAL PDstandardNth1B3 = INITVALUE; - CCTK_REAL PDstandardNth2B3 = INITVALUE; - CCTK_REAL PDstandardNth3B3 = INITVALUE; + CCTK_REAL PDstandardNth1beta1 = INITVALUE; + CCTK_REAL PDstandardNth2beta1 = INITVALUE; + CCTK_REAL PDstandardNth3beta1 = INITVALUE; + CCTK_REAL PDstandardNth1beta2 = INITVALUE; + CCTK_REAL PDstandardNth2beta2 = INITVALUE; + CCTK_REAL PDstandardNth3beta2 = INITVALUE; + CCTK_REAL PDstandardNth1beta3 = INITVALUE; + CCTK_REAL PDstandardNth2beta3 = INITVALUE; + CCTK_REAL PDstandardNth3beta3 = INITVALUE; CCTK_REAL PDstandardNth1gt11 = INITVALUE; CCTK_REAL PDstandardNth2gt11 = INITVALUE; CCTK_REAL PDstandardNth3gt11 = INITVALUE; @@ -145,9 +145,6 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa /* Assign local copies of grid functions */ alphaL = alpha[index]; - B1L = B1[index]; - B2L = B2[index]; - B3L = B3[index]; beta1L = beta1[index]; beta2L = beta2[index]; beta3L = beta3[index]; @@ -167,15 +164,15 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa /* Include user supplied include files */ /* Precompute derivatives (new style) */ - PDstandardNth1B1 = PDstandardNth1(B1, i, j, k); - PDstandardNth2B1 = PDstandardNth2(B1, i, j, k); - PDstandardNth3B1 = PDstandardNth3(B1, i, j, k); - PDstandardNth1B2 = PDstandardNth1(B2, i, j, k); - PDstandardNth2B2 = PDstandardNth2(B2, i, j, k); - PDstandardNth3B2 = PDstandardNth3(B2, i, j, k); - PDstandardNth1B3 = PDstandardNth1(B3, i, j, k); - PDstandardNth2B3 = PDstandardNth2(B3, i, j, k); - PDstandardNth3B3 = PDstandardNth3(B3, i, j, k); + PDstandardNth1beta1 = PDstandardNth1(beta1, i, j, k); + PDstandardNth2beta1 = PDstandardNth2(beta1, i, j, k); + PDstandardNth3beta1 = PDstandardNth3(beta1, i, j, k); + PDstandardNth1beta2 = PDstandardNth1(beta2, i, j, k); + PDstandardNth2beta2 = PDstandardNth2(beta2, i, j, k); + PDstandardNth3beta2 = PDstandardNth3(beta2, i, j, k); + PDstandardNth1beta3 = PDstandardNth1(beta3, i, j, k); + PDstandardNth2beta3 = PDstandardNth2(beta3, i, j, k); + PDstandardNth3beta3 = PDstandardNth3(beta3, i, j, k); PDstandardNth1gt11 = PDstandardNth1(gt11, i, j, k); PDstandardNth2gt11 = PDstandardNth2(gt11, i, j, k); PDstandardNth3gt11 = PDstandardNth3(gt11, i, j, k); @@ -272,16 +269,16 @@ void ML_BSSN_convertFromADMBaseGamma_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT fa Xt3L = Gt311*gtu11 + Gt322*gtu22 + 2*(Gt312*gtu21 + Gt313*gtu31 + Gt323*gtu32) + Gt333*gtu33; - AL = -1.*dtalpL*(-1. + LapseAdvectionCoeff)*INV(harmonicF)*pow(alphaL,-harmonicN); + AL = -(dtalpL*(-1 + LapseAdvectionCoeff)*INV(harmonicF)*pow(alphaL,-harmonicN)); - B1L = (dtbetaxL - (beta1L*PDstandardNth1B1 + beta2L*PDstandardNth2B1 + beta3L*PDstandardNth3B1)*ShiftAdvectionCoeff)* - INV(ShiftGammaCoeff); + B1L = (dtbetaxL - (beta1L*PDstandardNth1beta1 + beta2L*PDstandardNth2beta1 + beta3L*PDstandardNth3beta1)* + ShiftAdvectionCoeff)*INV(ShiftGammaCoeff); - B2L = (dtbetayL - (beta1L*PDstandardNth1B2 + beta2L*PDstandardNth2B2 + beta3L*PDstandardNth3B2)*ShiftAdvectionCoeff)* - INV(ShiftGammaCoeff); + B2L = (dtbetayL - (beta1L*PDstandardNth1beta2 + beta2L*PDstandardNth2beta2 + beta3L*PDstandardNth3beta2)* + ShiftAdvectionCoeff)*INV(ShiftGammaCoeff); - B3L = (dtbetazL - (beta1L*PDstandardNth1B3 + beta2L*PDstandardNth2B3 + beta3L*PDstandardNth3B3)*ShiftAdvectionCoeff)* - INV(ShiftGammaCoeff); + B3L = (dtbetazL - (beta1L*PDstandardNth1beta3 + beta2L*PDstandardNth2beta3 + beta3L*PDstandardNth3beta3)* + ShiftAdvectionCoeff)*INV(ShiftGammaCoeff); /* Copy local copies back to grid functions */ diff --git a/ML_BSSN/src/ML_BSSN_matter.c b/ML_BSSN/src/ML_BSSN_matter.c deleted file mode 100644 index d08c254..0000000 --- a/ML_BSSN/src/ML_BSSN_matter.c +++ /dev/null @@ -1,298 +0,0 @@ -/* File produced by user eschnett */ -/* Produced with Mathematica Version 6.0 for Mac OS X x86 (32-bit) (April 20, 2007) */ - -/* Mathematica script written by Ian Hinder and Sascha Husa */ - -#define KRANC_C - -#include <math.h> -#include "cctk.h" -#include "cctk_Arguments.h" -#include "cctk_Parameters.h" -#include "GenericFD.h" -#include "Differencing.h" -#include "loopcontrol.h" - -/* Define macros used in calculations */ -#define INITVALUE (42) -#define INV(x) ((1.0) / (x)) -#define SQR(x) ((x) * (x)) -#define CUB(x) ((x) * (x) * (x)) -#define QAD(x) ((x) * (x) * (x) * (x)) - -void ML_BSSN_matter_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[3], CCTK_REAL tangentA[3], CCTK_REAL tangentB[3], CCTK_INT min[3], CCTK_INT max[3], CCTK_INT n_subblock_gfs, CCTK_REAL *subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - - /* Declare finite differencing variables */ - CCTK_REAL dx = INITVALUE, dy = INITVALUE, dz = INITVALUE; - CCTK_REAL dxi = INITVALUE, dyi = INITVALUE, dzi = INITVALUE; - CCTK_REAL khalf = INITVALUE, kthird = INITVALUE, ktwothird = INITVALUE, kfourthird = INITVALUE, keightthird = INITVALUE; - CCTK_REAL hdxi = INITVALUE, hdyi = INITVALUE, hdzi = INITVALUE; - - - /* Declare predefined quantities */ - CCTK_REAL p1o12dx = INITVALUE; - CCTK_REAL p1o12dy = INITVALUE; - CCTK_REAL p1o12dz = INITVALUE; - CCTK_REAL p1o144dxdy = INITVALUE; - CCTK_REAL p1o144dxdz = INITVALUE; - CCTK_REAL p1o144dydz = INITVALUE; - CCTK_REAL pm1o12dx2 = INITVALUE; - CCTK_REAL pm1o12dy2 = INITVALUE; - CCTK_REAL pm1o12dz2 = INITVALUE; - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_matter_Body"); - } - - if (cctk_iteration % ML_BSSN_matter_calc_every != ML_BSSN_matter_calc_offset) - { - return; - } - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - dx = CCTK_DELTA_SPACE(0); - dy = CCTK_DELTA_SPACE(1); - dz = CCTK_DELTA_SPACE(2); - dxi = 1.0 / dx; - dyi = 1.0 / dy; - dzi = 1.0 / dz; - khalf = 0.5; - kthird = 1/3.0; - ktwothird = 2.0/3.0; - kfourthird = 4.0/3.0; - keightthird = 8.0/3.0; - hdxi = 0.5 * dxi; - hdyi = 0.5 * dyi; - hdzi = 0.5 * dzi; - - /* Initialize predefined quantities */ - p1o12dx = INV(dx)/12.; - p1o12dy = INV(dy)/12.; - p1o12dz = INV(dz)/12.; - p1o144dxdy = (INV(dx)*INV(dy))/144.; - p1o144dxdz = (INV(dx)*INV(dz))/144.; - p1o144dydz = (INV(dy)*INV(dz))/144.; - pm1o12dx2 = -pow(dx,-2)/12.; - pm1o12dy2 = -pow(dy,-2)/12.; - pm1o12dz2 = -pow(dz,-2)/12.; - - /* Loop over the grid points */ - _Pragma ("omp parallel") - LC_LOOP3 (ML_BSSN_matter, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) - { - int index = INITVALUE; - int subblock_index = INITVALUE; - index = CCTK_GFINDEX3D(cctkGH,i,j,k); - subblock_index = i - min[0] + (max[0] - min[0]) * (j - min[1] + (max[1]-min[1]) * (k - min[2])); - - /* Declare shorthands */ - CCTK_REAL detgt = INITVALUE; - CCTK_REAL e4phi = INITVALUE; - CCTK_REAL em4phi = INITVALUE; - CCTK_REAL g11 = INITVALUE, g12 = INITVALUE, g13 = INITVALUE, g22 = INITVALUE, g23 = INITVALUE, g33 = INITVALUE; - CCTK_REAL gtu11 = INITVALUE, gtu21 = INITVALUE, gtu22 = INITVALUE, gtu31 = INITVALUE, gtu32 = INITVALUE, gtu33 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE; - CCTK_REAL rho = INITVALUE; - CCTK_REAL S1 = INITVALUE, S2 = INITVALUE, S3 = INITVALUE; - CCTK_REAL T00 = INITVALUE, T01 = INITVALUE, T02 = INITVALUE, T03 = INITVALUE, T11 = INITVALUE, T12 = INITVALUE; - CCTK_REAL T13 = INITVALUE, T22 = INITVALUE, T23 = INITVALUE, T33 = INITVALUE; - CCTK_REAL trS = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL At11rhsL = INITVALUE, At12rhsL = INITVALUE, At13rhsL = INITVALUE, At22rhsL = INITVALUE, At23rhsL = INITVALUE, At33rhsL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL eTttL = INITVALUE; - CCTK_REAL eTtxL = INITVALUE; - CCTK_REAL eTtyL = INITVALUE; - CCTK_REAL eTtzL = INITVALUE; - CCTK_REAL eTxxL = INITVALUE; - CCTK_REAL eTxyL = INITVALUE; - CCTK_REAL eTxzL = INITVALUE; - CCTK_REAL eTyyL = INITVALUE; - CCTK_REAL eTyzL = INITVALUE; - CCTK_REAL eTzzL = INITVALUE; - CCTK_REAL gt11L = INITVALUE, gt12L = INITVALUE, gt13L = INITVALUE, gt22L = INITVALUE, gt23L = INITVALUE, gt33L = INITVALUE; - CCTK_REAL phiL = INITVALUE; - CCTK_REAL trKrhsL = INITVALUE; - CCTK_REAL Xt1rhsL = INITVALUE, Xt2rhsL = INITVALUE, Xt3rhsL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - alphaL = alpha[index]; - At11rhsL = At11rhs[index]; - At12rhsL = At12rhs[index]; - At13rhsL = At13rhs[index]; - At22rhsL = At22rhs[index]; - At23rhsL = At23rhs[index]; - At33rhsL = At33rhs[index]; - beta1L = beta1[index]; - beta2L = beta2[index]; - beta3L = beta3[index]; - eTttL = eTtt[index]; - eTtxL = eTtx[index]; - eTtyL = eTty[index]; - eTtzL = eTtz[index]; - eTxxL = eTxx[index]; - eTxyL = eTxy[index]; - eTxzL = eTxz[index]; - eTyyL = eTyy[index]; - eTyzL = eTyz[index]; - eTzzL = eTzz[index]; - gt11L = gt11[index]; - gt12L = gt12[index]; - gt13L = gt13[index]; - gt22L = gt22[index]; - gt23L = gt23[index]; - gt33L = gt33[index]; - phiL = phi[index]; - trKrhsL = trKrhs[index]; - Xt1rhsL = Xt1rhs[index]; - Xt2rhsL = Xt2rhs[index]; - Xt3rhsL = Xt3rhs[index]; - - /* Assign local copies of subblock grid functions */ - - /* Include user supplied include files */ - - /* Precompute derivatives (new style) */ - - /* Precompute derivatives (old style) */ - - /* Calculate temporaries and grid functions */ - T00 = eTttL; - - T01 = eTtxL; - - T02 = eTtyL; - - T03 = eTtzL; - - T11 = eTxxL; - - T12 = eTxyL; - - T13 = eTxzL; - - T22 = eTyyL; - - T23 = eTyzL; - - T33 = eTzzL; - - rho = pow(alphaL,-2)*(T00 - 2*(beta2L*T02 + beta3L*T03) + - 2*(beta1L*(-T01 + beta2L*T12 + beta3L*T13) + beta2L*beta3L*T23) + T11*SQR(beta1L) + T22*SQR(beta2L) + - T33*SQR(beta3L)); - - S1 = (-T01 + beta1L*T11 + beta2L*T12 + beta3L*T13)*INV(alphaL); - - S2 = (-T02 + beta1L*T12 + beta2L*T22 + beta3L*T23)*INV(alphaL); - - S3 = (-T03 + beta1L*T13 + beta2L*T23 + beta3L*T33)*INV(alphaL); - - detgt = 1; - - gtu11 = INV(detgt)*(gt22L*gt33L - SQR(gt23L)); - - gtu21 = (gt13L*gt23L - gt12L*gt33L)*INV(detgt); - - gtu31 = (-(gt13L*gt22L) + gt12L*gt23L)*INV(detgt); - - gtu22 = INV(detgt)*(gt11L*gt33L - SQR(gt13L)); - - gtu32 = (gt12L*gt13L - gt11L*gt23L)*INV(detgt); - - gtu33 = INV(detgt)*(gt11L*gt22L - SQR(gt12L)); - - e4phi = exp(4*phiL); - - em4phi = INV(e4phi); - - g11 = e4phi*gt11L; - - g12 = e4phi*gt12L; - - g13 = e4phi*gt13L; - - g22 = e4phi*gt22L; - - g23 = e4phi*gt23L; - - g33 = e4phi*gt33L; - - gu11 = em4phi*gtu11; - - gu21 = em4phi*gtu21; - - gu31 = em4phi*gtu31; - - gu22 = em4phi*gtu22; - - gu32 = em4phi*gtu32; - - gu33 = em4phi*gtu33; - - trS = gu11*T11 + gu22*T22 + 2*(gu21*T12 + gu31*T13 + gu32*T23) + gu33*T33; - - trKrhsL = trKrhsL + 12.56637061435917295385057353311801153679*alphaL*(rho + trS); - - At11rhsL = At11rhsL + alphaL*em4phi*(-25.13274122871834590770114706623602307358*T11 + - 8.377580409572781969233715688745341024526*g11*trS); - - At12rhsL = At12rhsL + alphaL*em4phi*(-25.13274122871834590770114706623602307358*T12 + - 8.377580409572781969233715688745341024526*g12*trS); - - At13rhsL = At13rhsL + alphaL*em4phi*(-25.13274122871834590770114706623602307358*T13 + - 8.377580409572781969233715688745341024526*g13*trS); - - At22rhsL = At22rhsL + alphaL*em4phi*(-25.13274122871834590770114706623602307358*T22 + - 8.377580409572781969233715688745341024526*g22*trS); - - At23rhsL = At23rhsL + alphaL*em4phi*(-25.13274122871834590770114706623602307358*T23 + - 8.377580409572781969233715688745341024526*g23*trS); - - At33rhsL = At33rhsL + alphaL*em4phi*(-25.13274122871834590770114706623602307358*T33 + - 8.377580409572781969233715688745341024526*g33*trS); - - Xt1rhsL = -50.26548245743669181540229413247204614715*alphaL*(gtu11*S1 + gtu21*S2 + gtu31*S3) + Xt1rhsL; - - Xt2rhsL = -50.26548245743669181540229413247204614715*alphaL*(gtu21*S1 + gtu22*S2 + gtu32*S3) + Xt2rhsL; - - Xt3rhsL = -50.26548245743669181540229413247204614715*alphaL*(gtu31*S1 + gtu32*S2 + gtu33*S3) + Xt3rhsL; - - - /* Copy local copies back to grid functions */ - At11rhs[index] = At11rhsL; - At12rhs[index] = At12rhsL; - At13rhs[index] = At13rhsL; - At22rhs[index] = At22rhsL; - At23rhs[index] = At23rhsL; - At33rhs[index] = At33rhsL; - trKrhs[index] = trKrhsL; - Xt1rhs[index] = Xt1rhsL; - Xt2rhs[index] = Xt2rhsL; - Xt3rhs[index] = Xt3rhsL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (ML_BSSN_matter); -} - -void ML_BSSN_matter(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverInterior(cctkGH, &ML_BSSN_matter_Body); -} diff --git a/ML_BSSN/src/make.code.defn b/ML_BSSN/src/make.code.defn index 2d1b2a0..869e3ab 100644 --- a/ML_BSSN/src/make.code.defn +++ b/ML_BSSN/src/make.code.defn @@ -3,4 +3,4 @@ # Mathematica script written by Ian Hinder and Sascha Husa -SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_Minkowski.c ML_BSSN_convertFromADMBase.c ML_BSSN_convertFromADMBaseGamma.c ML_BSSN_RHS.c ML_BSSN_matter.c ML_BSSN_enforce.c ML_BSSN_convertToADMBase.c ML_BSSN_constraints.c ML_BSSN_matter_constraints.c Boundaries.c +SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_Minkowski.c ML_BSSN_convertFromADMBase.c ML_BSSN_convertFromADMBaseGamma.c ML_BSSN_RHS.c ML_BSSN_enforce.c ML_BSSN_boundary.c ML_BSSN_convertToADMBase.c ML_BSSN_constraints.c ML_BSSN_constraints_boundary.c Boundaries.c diff --git a/ML_BSSN_Helper/configuration.ccl b/ML_BSSN_Helper/configuration.ccl new file mode 100644 index 0000000..aecb1c4 --- /dev/null +++ b/ML_BSSN_Helper/configuration.ccl @@ -0,0 +1 @@ +REQUIRES THORNS: CoordGauge diff --git a/ML_BSSN_Helper/interface.ccl b/ML_BSSN_Helper/interface.ccl index 7fda368..0e82975 100644 --- a/ML_BSSN_Helper/interface.ccl +++ b/ML_BSSN_Helper/interface.ccl @@ -1 +1,3 @@ IMPLEMENTS: ML_BSSN_Helper + +INHERITS: ADMBase CoordGauge ML_BSSN diff --git a/ML_BSSN_Helper/schedule.ccl b/ML_BSSN_Helper/schedule.ccl index 49edf5e..ec388d1 100644 --- a/ML_BSSN_Helper/schedule.ccl +++ b/ML_BSSN_Helper/schedule.ccl @@ -1,7 +1,48 @@ if (CCTK_EQUALS (evolution_method, "ML_BSSN")) { + + STORAGE: ADMBase::metric[3] + STORAGE: ADMBase::curv[3] + STORAGE: ADMBase::lapse[3] + STORAGE: ADMBase::shift[3] + STORAGE: ADMBase::dtlapse[3] + STORAGE: ADMBase::dtshift[3] + + SCHEDULE ML_BSSN_RegisterSlicing AT startup + { + LANG: C + OPTIONS: meta + } "Register slicing" + SCHEDULE ML_BSSN_UnsetCheckpointTags AT basegrid { LANG: C - OPTIONS: global + OPTIONS: meta } "Don't checkpoint ADMBase variables" + + + + SCHEDULE GROUP ML_BSSN_evolCalcGroup IN MoL_CalcRHS + { + } "Calculate BSSN RHS" + + SCHEDULE GROUP ML_BSSN_evolCalcGroup AT analysis + { + TRIGGERS: ML_BSSN::ML_log_confacrhs + TRIGGERS: ML_BSSN::ML_metricrhs + TRIGGERS: ML_BSSN::ML_Gammarhs + TRIGGERS: ML_BSSN::ML_trace_curvrhs + TRIGGERS: ML_BSSN::ML_curvrhs + TRIGGERS: ML_BSSN::ML_lapserhs + TRIGGERS: ML_BSSN::ML_dtlapserhs + TRIGGERS: ML_BSSN::ML_shiftrhs + TRIGGERS: ML_BSSN::ML_dtshiftrhs + } "Calculate BSSN RHS" + + + + SCHEDULE GROUP ML_BSSN_constraintsCalcGroup AT analysis + { + TRIGGERS: ML_BSSN::Ham + TRIGGERS: ML_BSSN::mom + } "Calculate BSSN constraints" } diff --git a/ML_BSSN_Helper/src/RegisterSlicing.c b/ML_BSSN_Helper/src/RegisterSlicing.c new file mode 100644 index 0000000..313bf32 --- /dev/null +++ b/ML_BSSN_Helper/src/RegisterSlicing.c @@ -0,0 +1,10 @@ +#include <cctk.h> + +#include "CactusEinstein/CoordGauge/src/Slicing.h" + +int +ML_BSSN_RegisterSlicing (void) +{ + Einstein_RegisterSlicing ("ML_BSSN"); + return 0; +} diff --git a/ML_BSSN_Helper/src/make.code.defn b/ML_BSSN_Helper/src/make.code.defn index 7df7001..ec09a08 100644 --- a/ML_BSSN_Helper/src/make.code.defn +++ b/ML_BSSN_Helper/src/make.code.defn @@ -1,2 +1,2 @@ # -*-Makefile-*- -SRCS = UnsetCheckpointTags.c +SRCS = RegisterSlicing.c UnsetCheckpointTags.c diff --git a/ML_FOWaveToy/param.ccl b/ML_FOWaveToy/param.ccl index 8a18614..7256e37 100644 --- a/ML_FOWaveToy/param.ccl +++ b/ML_FOWaveToy/param.ccl @@ -36,6 +36,12 @@ CCTK_INT ML_FOWaveToy_MaxNumConstrainedVars "Number of constrained variables use 61:61 :: "Number of constrained variables used by this thorn" } 61 +private: +CCTK_INT timelevels "Number of active timelevels" +{ + 0:2 :: "" +} 2 + restricted: CCTK_INT WTFO_Gaussian_calc_every "WTFO_Gaussian_calc_every" { diff --git a/ML_FOWaveToy/schedule.ccl b/ML_FOWaveToy/schedule.ccl index e4b8085..7703b78 100644 --- a/ML_FOWaveToy/schedule.ccl +++ b/ML_FOWaveToy/schedule.ccl @@ -12,11 +12,32 @@ STORAGE: WT_urhs[1] STORAGE: WT_vrhs[1] -STORAGE: WT_rho[2] +if (timelevels == 1) +{ + STORAGE: WT_rho[1] +} +if (timelevels == 2) +{ + STORAGE: WT_rho[2] +} -STORAGE: WT_u[2] +if (timelevels == 1) +{ + STORAGE: WT_u[1] +} +if (timelevels == 2) +{ + STORAGE: WT_u[2] +} -STORAGE: WT_v[2] +if (timelevels == 1) +{ + STORAGE: WT_v[1] +} +if (timelevels == 2) +{ + STORAGE: WT_v[2] +} schedule ML_FOWaveToy_Startup at STARTUP { @@ -39,19 +60,16 @@ schedule ML_FOWaveToy_RegisterSymmetries at BASEGRID schedule WTFO_Gaussian AT initial { LANG: C - } "WTFO_Gaussian" schedule WTFO_RHS IN MoL_CalcRHS { LANG: C - } "WTFO_RHS" schedule WTFO_RHS AT analysis { LANG: C - SYNC: WT_rhorhs SYNC: WT_urhs SYNC: WT_vrhs @@ -60,7 +78,6 @@ schedule WTFO_RHS AT analysis schedule WTFO_constraints AT analysis { LANG: C - SYNC: WT_w } "WTFO_constraints" @@ -82,5 +99,4 @@ schedule ML_FOWaveToy_CheckBoundaries at BASEGRID schedule group ApplyBCs as ML_FOWaveToy_ApplyBCs in MoL_PostStep after ML_FOWaveToy_ApplyBoundConds { # no language specified - } "Apply boundary conditions controlled by thorn Boundary" diff --git a/ML_FOWavetoy/param.ccl b/ML_FOWavetoy/param.ccl index 8a18614..7256e37 100644 --- a/ML_FOWavetoy/param.ccl +++ b/ML_FOWavetoy/param.ccl @@ -36,6 +36,12 @@ CCTK_INT ML_FOWaveToy_MaxNumConstrainedVars "Number of constrained variables use 61:61 :: "Number of constrained variables used by this thorn" } 61 +private: +CCTK_INT timelevels "Number of active timelevels" +{ + 0:2 :: "" +} 2 + restricted: CCTK_INT WTFO_Gaussian_calc_every "WTFO_Gaussian_calc_every" { diff --git a/ML_FOWavetoy/schedule.ccl b/ML_FOWavetoy/schedule.ccl index e4b8085..7703b78 100644 --- a/ML_FOWavetoy/schedule.ccl +++ b/ML_FOWavetoy/schedule.ccl @@ -12,11 +12,32 @@ STORAGE: WT_urhs[1] STORAGE: WT_vrhs[1] -STORAGE: WT_rho[2] +if (timelevels == 1) +{ + STORAGE: WT_rho[1] +} +if (timelevels == 2) +{ + STORAGE: WT_rho[2] +} -STORAGE: WT_u[2] +if (timelevels == 1) +{ + STORAGE: WT_u[1] +} +if (timelevels == 2) +{ + STORAGE: WT_u[2] +} -STORAGE: WT_v[2] +if (timelevels == 1) +{ + STORAGE: WT_v[1] +} +if (timelevels == 2) +{ + STORAGE: WT_v[2] +} schedule ML_FOWaveToy_Startup at STARTUP { @@ -39,19 +60,16 @@ schedule ML_FOWaveToy_RegisterSymmetries at BASEGRID schedule WTFO_Gaussian AT initial { LANG: C - } "WTFO_Gaussian" schedule WTFO_RHS IN MoL_CalcRHS { LANG: C - } "WTFO_RHS" schedule WTFO_RHS AT analysis { LANG: C - SYNC: WT_rhorhs SYNC: WT_urhs SYNC: WT_vrhs @@ -60,7 +78,6 @@ schedule WTFO_RHS AT analysis schedule WTFO_constraints AT analysis { LANG: C - SYNC: WT_w } "WTFO_constraints" @@ -82,5 +99,4 @@ schedule ML_FOWaveToy_CheckBoundaries at BASEGRID schedule group ApplyBCs as ML_FOWaveToy_ApplyBCs in MoL_PostStep after ML_FOWaveToy_ApplyBoundConds { # no language specified - } "Apply boundary conditions controlled by thorn Boundary" diff --git a/ML_WaveToy/param.ccl b/ML_WaveToy/param.ccl index 5964f1e..36cb15e 100644 --- a/ML_WaveToy/param.ccl +++ b/ML_WaveToy/param.ccl @@ -36,6 +36,12 @@ CCTK_INT ML_WaveToy_MaxNumConstrainedVars "Number of constrained variables used 58:58 :: "Number of constrained variables used by this thorn" } 58 +private: +CCTK_INT timelevels "Number of active timelevels" +{ + 0:2 :: "" +} 2 + restricted: CCTK_INT WT_Gaussian_calc_every "WT_Gaussian_calc_every" { diff --git a/ML_WaveToy/schedule.ccl b/ML_WaveToy/schedule.ccl index 5308b4c..b52fd02 100644 --- a/ML_WaveToy/schedule.ccl +++ b/ML_WaveToy/schedule.ccl @@ -8,9 +8,23 @@ STORAGE: WT_rhorhs[1] STORAGE: WT_urhs[1] -STORAGE: WT_rho[2] +if (timelevels == 1) +{ + STORAGE: WT_rho[1] +} +if (timelevels == 2) +{ + STORAGE: WT_rho[2] +} -STORAGE: WT_u[2] +if (timelevels == 1) +{ + STORAGE: WT_u[1] +} +if (timelevels == 2) +{ + STORAGE: WT_u[2] +} schedule ML_WaveToy_Startup at STARTUP { @@ -33,19 +47,16 @@ schedule ML_WaveToy_RegisterSymmetries at BASEGRID schedule WT_Gaussian AT initial { LANG: C - } "WT_Gaussian" schedule WT_RHS IN MoL_CalcRHS { LANG: C - } "WT_RHS" schedule WT_RHS AT analysis { LANG: C - SYNC: WT_rhorhs SYNC: WT_urhs } "WT_RHS" @@ -67,5 +78,4 @@ schedule ML_WaveToy_CheckBoundaries at BASEGRID schedule group ApplyBCs as ML_WaveToy_ApplyBCs in MoL_PostStep after ML_WaveToy_ApplyBoundConds { # no language specified - } "Apply boundary conditions controlled by thorn Boundary" diff --git a/m/McLachlan.m b/m/McLachlan.m index e1c4e60..0767025 100644 --- a/m/McLachlan.m +++ b/m/McLachlan.m @@ -10,24 +10,13 @@ SetSourceLanguage["C"]; (* Derivatives *) (******************************************************************************) +KD = KroneckerDelta; + +(* derivative order: 2, 4, 6, 8, ... *) derivOrder = 4; derivatives = { - (* - PDstandard2nd[i_] -> StandardCenteredDifferenceOperator[1,1,i], - PDstandard2nd[i_, i_] -> StandardCenteredDifferenceOperator[2,1,i], - PDstandard2nd[i_, j_] -> StandardCenteredDifferenceOperator[1,1,i] - StandardCenteredDifferenceOperator[1,1,j] - *) - - (* - PDstandard4th[i_] -> StandardCenteredDifferenceOperator[1,2,i], - PDstandard4th[i_, i_] -> StandardCenteredDifferenceOperator[2,2,i], - PDstandard4th[i_, j_] -> StandardCenteredDifferenceOperator[1,2,i] - StandardCenteredDifferenceOperator[1,2,j] - *) - PDstandardNth[i_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i], PDstandardNth[i_, i_] -> StandardCenteredDifferenceOperator[2,derivOrder/2,i], PDstandardNth[i_, j_] -> StandardCenteredDifferenceOperator[1,derivOrder/2,i] @@ -45,10 +34,11 @@ PDglob[var_,lx_,ly_] := UseGlobalDerivs = False; PD := If [UseGlobalDerivs, PDglob, PDloc]; -(* timelevels *) +(* timelevels: 2 or 3 *) evolutionTimelevels = 3; -KD = KroneckerDelta; +(* matter: 0 or 1 *) +addMatter = 0; (******************************************************************************) (* Tensors *) @@ -176,8 +166,6 @@ initialCalc = Name -> "ML_ADM_Minkowski", Schedule -> {"IN ADMBase_InitialData"}, ConditionalOnKeyword -> {"my_initial_data", "Minkowski"}, - (* Where -> Boundary, *) - (* Where -> Interior, *) Equations -> { g[la,lb] -> KD[la,lb], @@ -191,9 +179,7 @@ initialCalcBSSN = { Name -> "ML_BSSN_Minkowski", Schedule -> {"IN ADMBase_InitialData"}, - ConditionalOnKeyword -> {"initial_data", "ML_BSSN__Minkowski"}, - (* Where -> Boundary, *) - (* Where -> Interior, *) + ConditionalOnKeyword -> {"my_initial_data", "Minkowski"}, Equations -> { phi -> 0, @@ -243,7 +229,6 @@ convertFromADMBaseCalcBSSN = { Name -> "ML_BSSN_convertFromADMBase", Schedule -> {"AT initial AFTER ADMBase_PostInitial"}, - (* Should only happen if ADMBase::initial_data != ML_BSSN *) ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi, K[la,lb]}, Equations -> @@ -284,7 +269,6 @@ convertFromADMBaseCalcBSSNGamma = { Name -> "ML_BSSN_convertFromADMBaseGamma", Schedule -> {"AT initial AFTER ML_BSSN_convertFromADMBase"}, - (* Should only happen if ADMBase::initial_data != ML_BSSN *) ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, Where -> Interior, Shorthands -> {detgt, gtu[ua,ub], Gt[ua,lb,lc]}, @@ -296,23 +280,14 @@ convertFromADMBaseCalcBSSNGamma = (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]), Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc], - (* TODO: check this *) - (* A -> - (dtalp - Lie[alpha, beta]) / (harmonicF alpha^harmonicN), *) - A -> - dtalp / (harmonicF alpha^harmonicN) ( LapseAdvectionCoeff - 1.0 ), - - (* TODO: check this *) - (* B1 -> dtbetax / (ShiftGammaCoeff alpha^ShiftAlphaPower) *) - (* B2 -> dtbetay / (ShiftGammaCoeff alpha^ShiftAlphaPower) *) - (* B3 -> dtbetaz / (ShiftGammaCoeff alpha^ShiftAlphaPower) *) - (* B1 -> ShiftGammaCoeff dtbetax / ((ShiftGammaCoeff^2 + 1.0*^-100) alpha^ShiftAlphaPower), - B2 -> ShiftGammaCoeff dtbetay / ((ShiftGammaCoeff^2 + 1.0*^-100) alpha^ShiftAlphaPower), - B3 -> ShiftGammaCoeff dtbetaz / ((ShiftGammaCoeff^2 + 1.0*^-100) alpha^ShiftAlphaPower) *) - B1 -> 1/ShiftGammaCoeff ( + dtbetax - - ShiftAdvectionCoeff beta[ua] PD[B1,la] ), - B2 -> 1/ShiftGammaCoeff ( + dtbetay - - ShiftAdvectionCoeff beta[ua] PD[B2,la] ), - B3 -> 1/ShiftGammaCoeff ( + dtbetaz - - ShiftAdvectionCoeff beta[ua] PD[B3,la] ) + A -> - dtalp / (harmonicF alpha^harmonicN) (LapseAdvectionCoeff - 1), + + B1 -> 1/ShiftGammaCoeff + (dtbetax - ShiftAdvectionCoeff beta[ua] PD[beta1,la]), + B2 -> 1/ShiftGammaCoeff + (dtbetay - ShiftAdvectionCoeff beta[ua] PD[beta2,la]), + B3 -> 1/ShiftGammaCoeff + (dtbetaz - ShiftAdvectionCoeff beta[ua] PD[beta3,la]) } } @@ -323,7 +298,7 @@ convertFromADMBaseCalcBSSNGamma = convertToADMBaseCalc = { Name -> "ML_ADM_convertToADMBase", - Schedule -> {"IN MoL_PostStep AFTER ML_ADM_ApplyBCs"}, + Schedule -> {"IN MoL_PostStep AFTER (ML_ADM_ApplyBCs ML_ADM_boundary)"}, Equations -> { gxx -> g11, @@ -353,7 +328,7 @@ convertToADMBaseCalc = convertToADMBaseCalcBSSN = { Name -> "ML_BSSN_convertToADMBase", - Schedule -> {"IN MoL_PostStep AFTER ML_BSSN_ApplyBCs AFTER ML_BSSN_enforce"}, + Schedule -> {"IN MoL_PostStep AFTER (ML_BSSN_ApplyBCs ML_BSSN_boundary ML_BSSN_enforce)"}, ConditionalOnKeyword -> {"evolution_method", "ML_BSSN"}, Where -> Interior, Shorthands -> {e4phi, g[la,lb], K[la,lb]}, @@ -379,13 +354,8 @@ convertToADMBaseCalcBSSN = betay -> beta2, betaz -> beta3, (* see RHS *) -(* dtalp -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], - dtbetax -> ShiftGammaCoeff alpha^ShiftAlphaPower B1, - dtbetay -> ShiftGammaCoeff alpha^ShiftAlphaPower B2, - dtbetaz -> ShiftGammaCoeff alpha^ShiftAlphaPower B3 *) - - dtalp -> - harmonicF alpha^harmonicN ( - ( 1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK ) + dtalp -> - harmonicF alpha^harmonicN + ((1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK) + LapseAdvectionCoeff beta[ua] PD[alpha,la], dtbetax -> + ShiftGammaCoeff B1 + ShiftAdvectionCoeff beta[ub] PD[beta[ua],lb], @@ -431,17 +401,7 @@ evolCalc = evolCalcBSSN = { Name -> "ML_BSSN_RHS", - Schedule -> {"IN MoL_CalcRHS", "AT analysis"}, - ConditionalOnKeyword -> {"evolution_method", "ML_BSSN"}, - TriggerGroups -> {"ML_log_confacrhs", - "ML_metricrhs" , - "ML_Gammarhs" , - "ML_trace_curvrhs", - "ML_curvrhs" , - "ML_lapserhs" , - "ML_dtlapserhs" , - "ML_shiftrhs" , - "ML_dtshiftrhs" }, + Schedule -> {"IN ML_BSSN_evolCalcGroup"}, Where -> Interior, Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], dgtu[ua,ub,lc], ddgtu[ua,ub,lc,ld], Gt[ua,lb,lc], @@ -449,20 +409,9 @@ evolCalcBSSN = Atm[ua,lb], Atu[ua,ub], e4phi, em4phi, g[la,lb], detg, ddetg[la], gu[ua,ub], G[ua,lb,lc], Ats[la,lb], trAts, - T00, T0[la], T[la,lb]}, + T00, T0[la], T[la,lb], rho, S[la], trS}, Equations -> { - T00 -> eTtt, - T01 -> eTtx, - T02 -> eTty, - T03 -> eTtz, - T11 -> eTxx, - T12 -> eTxy, - T13 -> eTxz, - T22 -> eTyy, - T23 -> eTyz, - T33 -> eTzz, - detgt -> 1 (* detgtExpr *), ddetgt[la] -> 0 (* ddetgtExpr[la] *), @@ -488,11 +437,11 @@ evolCalcBSSN = + gtu[ul,um] (+ Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,lm] + Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,lm] + Gt[uk,li,lm] gt[lk,ln] Gt[un,ll,lj]), -(* Rt[li,lj] -> (1/2) ( - gtu[ul,um] PD[gt[li,lj],ll,lm] - + gt[lk,li] PD[Xt[uk],lj] + - + gt[lk,lj] PD[Xt[uk],li] - + Xtn[uk] gt[li,ln] Gt[un,lj,lk] - + Xtn[uk] gt[lj,ln] Gt[un,li,lk] ) +(* Rt[li,lj] -> (1/2) (- gtu[ul,um] PD[gt[li,lj],ll,lm] + + gt[lk,li] PD[Xt[uk],lj] + + + gt[lk,lj] PD[Xt[uk],li] + + Xtn[uk] gt[li,ln] Gt[un,lj,lk] + + Xtn[uk] gt[lj,ln] Gt[un,li,lk]) + gtu[ul,um] (+ Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,lm] + Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,lm] + Gt[uk,li,lm] gt[lk,ln] Gt[un,ll,lj]), *) @@ -514,14 +463,39 @@ evolCalcBSSN = gu[ua,ub] -> em4phi gtu[ua,ub], (* ddetg[la] -> 12 detg PD[phi,la], G[ua,lb,lc] -> Gt[ua,lb,lc] - + 1/(6 detg) ( KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb] - - gtu[ua,ud] gt[lb,lc] ddetg[ld] ), *) + + 1/(6 detg) (KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb] + - gtu[ua,ud] gt[lb,lc] ddetg[ld]), *) G[ua,lb,lc] -> Gt[ua,lb,lc] - + 2 ( KD[ua,lb] PD[phi,lc] + KD[ua,lc] PD[phi,lb] - - gtu[ua,ud] gt[lb,lc] PD[phi,ld] ), + + 2 (KD[ua,lb] PD[phi,lc] + KD[ua,lc] PD[phi,lb] + - gtu[ua,ud] gt[lb,lc] PD[phi,ld]), R[la,lb] -> Rt[la,lb] + Rphi[la,lb], + (* Matter terms *) + + T00 -> addMatter eTtt, + T01 -> addMatter eTtx, + T02 -> addMatter eTty, + T03 -> addMatter eTtz, + T11 -> addMatter eTxx, + T12 -> addMatter eTxy, + T13 -> addMatter eTxz, + T22 -> addMatter eTyy, + T23 -> addMatter eTyz, + T33 -> addMatter eTzz, + + (* rho = n^a n^b T_ab *) + rho -> addMatter + (1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj])), + + (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) + S[li] -> addMatter (-1/alpha (T0[li] - beta[uj] T[li,lj])), + + (* trS = gamma^ij T_ij *) + trS -> addMatter (gu[ui,uj] T[li,lj]), + + (* RHS terms *) + (* PRD 62, 044034 (2000), eqn. (10) *) (* PRD 67 084023 (2003), eqn. (16) and (23) *) dot[phi] -> - (1/6) alpha trK @@ -551,34 +525,36 @@ evolCalcBSSN = + (1/3) gtu[ui,uj] PD[beta[ul],lj,ll] + beta[uj] PD[Xt[ui],lj] - Xtn[uj] PD[beta[ui],lj] - + (2/3) Xtn[ui] PD[beta[uj],lj], + + (2/3) Xtn[ui] PD[beta[uj],lj] + (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) + + addMatter (- 16 pi alpha gtu[ui,uj] S[lj]), (* PRD 62, 044034 (2000), eqn. (11) *) dot[trK] -> - gu[ua,ub] CD[alpha,la,lb] + alpha (Atm[ua,lb] Atm[ub,la] + (1/3) trK^2) - + Lie[trK, beta], - + + Lie[trK, beta] + (* Equation (4.21) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) + + addMatter (4 pi alpha (rho + trS)), + (* PRD 62, 044034 (2000), eqn. (12) *) (* TODO: use Hamiltonian constraint to make tracefree *) Ats[la,lb] -> - CD[alpha,la,lb] + alpha R[la,lb], trAts -> gu[ua,ub] Ats[la,lb], dot[At[la,lb]] -> + em4phi (+ Ats[la,lb] - (1/3) g[la,lb] trAts ) + alpha (trK At[la,lb] - 2 At[la,lc] Atm[uc,lb]) - + Lie[At[la,lb], beta] - (2/3) At[la,lb] PD[beta[uc],lc], -(* dot[At[la,lb]] -> + em4phi (+ (- CD[alpha,la,lb] + alpha R[la,lb]) - - (1/3) g[la,lb] gu[uc,ud] - (- CD[alpha,lc,ld] + alpha R[lc,ld])) - + alpha (trK At[la,lb] - 2 At[la,lc] Atm[uc,lb]) - + Lie[At[la,lb], beta] - (2/3) At[la,lb] PD[beta[uc],lc], *) + + Lie[At[la,lb], beta] - (2/3) At[la,lb] PD[beta[uc],lc] + (* Equation (4.23) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) + + addMatter (- em4phi alpha 8 pi + (T[la,lb] - (1/3) g[la,lb] trS)), (* dot[alpha] -> - harmonicF alpha^harmonicN trK, *) (* dot[alpha] -> - harmonicF alpha^harmonicN A + Lie[alpha, beta], *) - dot[alpha] -> - harmonicF alpha^harmonicN ( - ( 1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK ) + dot[alpha] -> - harmonicF alpha^harmonicN ( + (1 - LapseAdvectionCoeff) A + LapseAdvectionCoeff trK) + LapseAdvectionCoeff beta[ua] PD[alpha,la], (* TODO: is the above Lie derivative correct? *) - dot[A] -> ( 1 - LapseAdvectionCoeff ) ( dot[trK] - AlphaDriver A ), + dot[A] -> (1 - LapseAdvectionCoeff) (dot[trK] - AlphaDriver A), (* dot[beta[ua]] -> eta Xt[ua], *) (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *) @@ -586,85 +562,12 @@ evolCalcBSSN = + ShiftAdvectionCoeff beta[ub] PD[beta[ua],lb], dot[B[ua]] -> + dot[Xt[ua]] - BetaDriver B[ua] - + ShiftAdvectionCoeff beta[ub] ( + PD[B[ua],lb] - - PD[Xt[ua],lb] ) + + ShiftAdvectionCoeff beta[ub] (+ PD[B[ua],lb] + - PD[Xt[ua],lb]) (* TODO: is there a Lie derivative of the shift missing? *) } } -addMatterBSSN = -{ - Name -> "ML_BSSN_matter", - Schedule -> {"IN MoL_CalcRHS AFTER ML_BSSN_RHS"}, - (* ConditionalOnKeyword -> {"evolution_method", "ML_BSSN"}, *) - ConditionalOnKeyword -> {"SpaceTime", "Space+Matter"}, - (* Conditional -> {Textual -> "stress_energy_storage && stress_energy_at_RHS"}, *) - (*TriggerGroups -> {"ML_log_confacrhs", - "ML_metricrhs" , - "ML_Gammarhs" , - "ML_trace_curvrhs", - "ML_curvrhs" , - "ML_lapserhs" , - "ML_dtlapserhs" , - "ML_shiftrhs" , - "ML_dtshiftrhs" },*) - Where -> Interior, - Shorthands -> {T00, T0[la], T[la,lb], rho, S[la], trS, - detgt, gtu[ua,ub], e4phi, em4phi, g[la,lb], gu[ua,ub]}, - Equations -> - { - - T00 -> eTtt, - T01 -> eTtx, - T02 -> eTty, - T03 -> eTtz, - T11 -> eTxx, - T12 -> eTxy, - T13 -> eTxz, - T22 -> eTyy, - T23 -> eTyz, - T33 -> eTzz, - - (* rho = n^a n^b T_ab *) - rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj]), - - (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) - S[li] -> -1/alpha ( T0[li] - beta[uj] T[li,lj] ), - - (* trS = gamma^ij T_ij *) - detgt -> 1, - gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], - e4phi -> Exp [4 phi], - em4phi -> 1 / e4phi, - g[la,lb] -> e4phi gt[la,lb], - gu[ua,ub] -> em4phi gtu[ua,ub], - trS -> gu[ui,uj] T[li,lj], - - (* Equation (4.21) in Baumgarte & Shapiro (Phys.Rept. 376 (2003) 41-131) *) - dot[trK] -> dot[trK] + 4 pi alpha ( rho + trS ), - - (* Equation (4.23) in Baumgarte & Shapiro (Phys.Rept. 376 (2003) 41-131) *) - dot[At[la,lb]] -> dot[At[la,lb]] - - em4phi alpha 8 pi ( T[la,lb] - (1/3) g[la,lb] trS ), - - (* Equation (4.28) in Baumgarte & Shapiro (Phys.Rept. 376 (2003) 41-131) *) - dot[Xt[ui]] -> dot[Xt[ui]] - 16 pi alpha gtu[ui,uj] S[lj] - } -} - -(* -evolCalcAnalysisBSSN = - evolCalcBSSN - // DeleteCases [#, Schedule -> _] & - // Join [#, { Schedule -> {"AT analysis"}, - TriggerGroups -> {"ML_metricrhs"}}] &; -*) - -(* -addMatterAnalysisBSSN = - "AT analysis AFTER ML_BSSN_RHS"}, - *) - enforceCalcBSSN = { Name -> "ML_BSSN_enforce", @@ -685,6 +588,45 @@ enforceCalcBSSN = } (******************************************************************************) +(* Boundary conditions *) +(******************************************************************************) + +boundaryCalc = +{ + Name -> "ML_ADM_boundary", + Schedule -> {"IN MoL_PostStep"}, + ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"}, + Where -> Boundary, + Equations -> + { + g[la,lb] -> KD[la,lb], + K[la,lb] -> 0, + alpha -> 1, + beta[ua] -> 0 + } +} + +boundaryCalcBSSN = +{ + Name -> "ML_BSSN_boundary", + Schedule -> {"IN MoL_PostStep"}, + ConditionalOnKeyword -> {"my_boundary_condition", "Minkowski"}, + Where -> Boundary, + Equations -> + { + phi -> 0, + gt[la,lb] -> KD[la,lb], + trK -> 0, + At[la,lb] -> 0, + Xt[ua] -> 0, + alpha -> 1, + A -> 0, + beta[ua] -> 0, + B[ua] -> 0 + } +} + +(******************************************************************************) (* Constraint equations *) (******************************************************************************) @@ -692,6 +634,7 @@ constraintsCalc = { Name -> "ML_ADM_constraints", Schedule -> {"AT analysis"}, + TriggerGroups -> {"Ham", "mom"}, Where -> Interior, Shorthands -> {detg, gu[ua,ub], G[ua,lb,lc], R[la,lb], trR, Km[ua,lb], trK}, Equations -> @@ -712,17 +655,29 @@ constraintsCalc = } } +constraintsBoundaryCalc = +{ + Name -> "ML_ADM_constraints_boundary", + Schedule -> {"AT analysis AFTER ML_ADM_constraints"}, + (* TriggerGroups -> {"Ham", "mom"}, *) + Where -> Boundary, + Equations -> + { + H -> 0, + M[la] -> 0 + } +} + constraintsCalcBSSN = { Name -> "ML_BSSN_constraints", - Schedule -> {"AT analysis"}, - ConditionalOnKeyword -> {"evolution_method", "ML_BSSN"}, - TriggerGroups -> {"Ham", "mom"}, + Schedule -> {"IN ML_BSSN_constraintsCalcGroup"}, Where -> Interior, Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], Gt[ua,lb,lc], e4phi, em4phi, g[la,lb], detg, gu[ua,ub], ddetg[la], G[ua,lb,lc], Rt[la,lb], Rphi[la,lb], R[la,lb], trR, Atm[la,lb], - gK[la,lb,lc]}, + gK[la,lb,lc], + T00, T0[la], T[la,lb], rho, S[la]}, Equations -> { detgt -> 1 (* detgtExpr *), @@ -793,18 +748,40 @@ constraintsCalcBSSN = (* Km[ua,lb] -> gu[ua,uc] K[lc,lb], *) Atm[ua,lb] -> gtu[ua,uc] At[lc,lb], + (* Matter terms *) + + T00 -> eTtt, + T01 -> eTtx, + T02 -> eTty, + T03 -> eTtz, + T11 -> eTxx, + T12 -> eTxy, + T13 -> eTxz, + T22 -> eTyy, + T23 -> eTyz, + T33 -> eTzz, + + (* rho = n^a n^b T_ab *) + rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj]), + + (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) + S[li] -> -1/alpha (T0[li] - beta[uj] T[li,lj]), + + (* Constraints *) + (* H -> trR - Km[ua,lb] Km[ub,la] + trK^2, *) (* PRD 67, 084023 (2003), eqn. (19) *) - H -> trR - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2, + H -> trR - Atm[ua,lb] Atm[ub,la] + (2/3) trK^2 - addMatter 16 pi rho, -(* (* gK[la,lb,lc] -> CD[K[la,lb],lc], *) - gK[la,lb,lc] -> + 4 e4phi PD[phi,lc] At[la,lb] + e4phi CD[At[la,lb],lc] + (* gK[la,lb,lc] -> CD[K[la,lb],lc], *) +(* gK[la,lb,lc] -> + 4 e4phi PD[phi,lc] At[la,lb] + e4phi CD[At[la,lb],lc] + (1/3) g[la,lb] PD[trK,lc], M[la] -> gu[ub,uc] (gK[lc,la,lb] - gK[lc,lb,la]), *) - M[li] -> + gtu[uj,uk] ( CDt[At[li,lj],lk] + 6 At[li,lj] PD[phi,lk] ) - - (2/3) PD[trK,li], + M[li] -> + gtu[uj,uk] (CDt[At[li,lj],lk] + 6 At[li,lj] PD[phi,lk]) + - (2/3) PD[trK,li] + - addMatter 8 pi S[li], (* TODO: use PRD 67, 084023 (2003), eqn. (20) *) (* det gamma-tilde *) @@ -818,38 +795,15 @@ constraintsCalcBSSN = } } -constraintsMatterBSSN = +constraintsBoundaryCalcBSSN = { - Name -> "ML_BSSN_matter_constraints", - Schedule -> {"AT analysis AFTER ML_BSSN_constraints"}, - (* ConditionalOnKeyword -> {"evolution_method", "ML_BSSN"}, *) - ConditionalOnKeyword -> {"SpaceTime", "Space+Matter"}, - TriggerGroups -> {"Ham", "mom"}, - Where -> Interior, - Shorthands -> {T00, T0[la], T[la,lb], rho, S[la]}, + Name -> "ML_BSSN_constraints_boundary", + Schedule -> {"IN ML_BSSN_constraintsCalcGroup AFTER ML_BSSN_constraints"}, + Where -> Boundary, Equations -> { - - T00 -> eTtt, - T01 -> eTtx, - T02 -> eTty, - T03 -> eTtz, - T11 -> eTxx, - T12 -> eTxy, - T13 -> eTxz, - T22 -> eTyy, - T23 -> eTyz, - T33 -> eTzz, - - (* rho = n^a n^b T_ab *) - rho -> 1/alpha^2 (T00 - 2 beta[ui] T0[li] + beta[ui] beta[uj] T[li,lj]), - - (* S_i = -p^a_i n^b T_ab, where p^a_i = delta^a_i + n^a n_i *) - S[li] -> -1/alpha ( T0[li] - beta[uj] T[li,lj] ), - - H -> H - 16 pi rho, - - M[li] -> M[li] - 8 pi S[li] + H -> 0, + M[la] -> 0 } } @@ -863,32 +817,11 @@ inheritedImplementations = {"ADMBase", "TmunuBase"}; (* Parameters *) (******************************************************************************) -(* inheritedKeywordParameters = { "ADMBase::initial_data" }; *) inheritedKeywordParameters = {}; extendedKeywordParameters = { { - Name -> "ADMBase::initial_data", - AllowedValues -> {"ML_BSSN__Minkowski"} - }, - { - Name -> "ADMBase::initial_lapse", - AllowedValues -> {"ML_BSSN__Minkowski"} - }, - { - Name -> "ADMBase::initial_shift", - AllowedValues -> {"ML_BSSN__Minkowski"} - }, - { - Name -> "ADMBase::initial_dtlapse", - AllowedValues -> {"ML_BSSN__Minkowski"} - }, - { - Name -> "ADMBase::initial_dtshift", - AllowedValues -> {"ML_BSSN__Minkowski"} - }, - { Name -> "ADMBase::evolution_method", AllowedValues -> {"ML_BSSN"} }, @@ -912,11 +845,11 @@ keywordParameters = Default -> "ADMBase" }, { - Name -> "SpaceTime", + Name -> "my_boundary_condition", (* Visibility -> "restricted", *) (* Description -> "ddd", *) - AllowedValues -> {"Space", "Space+Matter"}, - Default -> "Space" + AllowedValues -> {"none", "Minkowski"}, + Default -> "none" } }; @@ -955,12 +888,12 @@ realParameters = { Name -> LapseAdvectionCoeff, Description -> "Factor in front of the shift advection terms in 1+log", - Default -> 1. + Default -> 1 }, { Name -> ShiftAdvectionCoeff, Description -> "Factor in front of the shift advection terms in gamma driver", - Default -> 1. + Default -> 1 } }; @@ -974,8 +907,10 @@ calculations = initialCalc, convertFromADMBaseCalc, evolCalc, + boundaryCalc, convertToADMBaseCalc, - constraintsCalc + constraintsCalc, + constraintsBoundaryCalc }; CreateKrancThornTT [groups, ".", "ML_ADM", @@ -994,11 +929,11 @@ calculationsBSSN = convertFromADMBaseCalcBSSN, convertFromADMBaseCalcBSSNGamma, evolCalcBSSN, - addMatterBSSN, enforceCalcBSSN, + boundaryCalcBSSN, convertToADMBaseCalcBSSN, constraintsCalcBSSN, - constraintsMatterBSSN + constraintsBoundaryCalcBSSN }; CreateKrancThornTT [groupsBSSN, ".", "ML_BSSN", |