diff options
32 files changed, 62 insertions, 1077 deletions
diff --git a/ML_BSSN/interface.ccl b/ML_BSSN/interface.ccl index a0d38b5..e257aee 100644 --- a/ML_BSSN/interface.ccl +++ b/ML_BSSN/interface.ccl @@ -25,12 +25,6 @@ CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT USES FUNCTION Boundary_SelectVarForBC public: -CCTK_REAL ML_BetaDriver type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' -{ - eta -} "ML_BetaDriver" - -public: CCTK_REAL ML_cons_detg type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=2.0000000000000000000' { cS diff --git a/ML_BSSN/param.ccl b/ML_BSSN/param.ccl index 2a90c6a..d13eea4 100644 --- a/ML_BSSN/param.ccl +++ b/ML_BSSN/param.ccl @@ -207,18 +207,6 @@ CCTK_INT ML_BSSN_convertFromADMBaseGamma_calc_every "ML_BSSN_convertFromADMBaseG } 1 restricted: -CCTK_INT ML_BSSN_setBetaDriverConstant_calc_every "ML_BSSN_setBetaDriverConstant_calc_every" -{ - *:* :: "" -} 1 - -restricted: -CCTK_INT ML_BSSN_setBetaDriverSpatial_calc_every "ML_BSSN_setBetaDriverSpatial_calc_every" -{ - *:* :: "" -} 1 - -restricted: CCTK_INT ML_BSSN_RHS_calc_every "ML_BSSN_RHS_calc_every" { *:* :: "" @@ -321,18 +309,6 @@ CCTK_INT ML_BSSN_convertFromADMBaseGamma_calc_offset "ML_BSSN_convertFromADMBase } 0 restricted: -CCTK_INT ML_BSSN_setBetaDriverConstant_calc_offset "ML_BSSN_setBetaDriverConstant_calc_offset" -{ - *:* :: "" -} 0 - -restricted: -CCTK_INT ML_BSSN_setBetaDriverSpatial_calc_offset "ML_BSSN_setBetaDriverSpatial_calc_offset" -{ - *:* :: "" -} 0 - -restricted: CCTK_INT ML_BSSN_RHS_calc_offset "ML_BSSN_RHS_calc_offset" { *:* :: "" diff --git a/ML_BSSN/schedule.ccl b/ML_BSSN/schedule.ccl index bd0f54d..f7a8179 100644 --- a/ML_BSSN/schedule.ccl +++ b/ML_BSSN/schedule.ccl @@ -1,8 +1,6 @@ # File produced by Kranc -STORAGE: ML_BetaDriver[1] - STORAGE: ML_cons_detg[1] STORAGE: ML_cons_Gamma[1] @@ -367,24 +365,6 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase")) } "ML_BSSN_convertFromADMBaseGamma" } - -if (CCTK_EQUALS(UseSpatialBetaDriver, "no")) -{ - schedule ML_BSSN_setBetaDriverConstant IN ML_BSSN_InitEta - { - LANG: C - } "ML_BSSN_setBetaDriverConstant" -} - - -if (CCTK_EQUALS(UseSpatialBetaDriver, "yes")) -{ - schedule ML_BSSN_setBetaDriverSpatial IN ML_BSSN_InitEta - { - LANG: C - } "ML_BSSN_setBetaDriverSpatial" -} - schedule ML_BSSN_RHS IN NoSuchGroup { LANG: C diff --git a/ML_BSSN/src/ML_BSSN_RHS.c b/ML_BSSN/src/ML_BSSN_RHS.c index 241c214..70d492b 100644 --- a/ML_BSSN/src/ML_BSSN_RHS.c +++ b/ML_BSSN/src/ML_BSSN_RHS.c @@ -109,6 +109,7 @@ void ML_BSSN_RHS_Body(cGH const * restrict const cctkGH, int const dir, int cons // CCTK_REAL dir1 = INITVALUE, dir2 = INITVALUE, dir3 = INITVALUE; // CCTK_REAL e4phi = INITVALUE; // CCTK_REAL em4phi = INITVALUE; + // CCTK_REAL eta = INITVALUE; // CCTK_REAL fac1 = INITVALUE, fac2 = INITVALUE; // CCTK_REAL g11 = INITVALUE; // CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE; @@ -144,7 +145,6 @@ void ML_BSSN_RHS_Body(cGH const * restrict const cctkGH, int const dir, int cons // CCTK_REAL At22L = INITVALUE, At22rhsL = INITVALUE, At23L = INITVALUE, At23rhsL = INITVALUE, At33L = INITVALUE, At33rhsL = INITVALUE; // CCTK_REAL B1L = INITVALUE, B1rhsL = INITVALUE, B2L = INITVALUE, B2rhsL = INITVALUE, B3L = INITVALUE, B3rhsL = INITVALUE; // CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE; - // CCTK_REAL etaL = INITVALUE; // CCTK_REAL eTttL = INITVALUE; // CCTK_REAL eTtxL = INITVALUE; // CCTK_REAL eTtyL = INITVALUE; @@ -158,6 +158,7 @@ void ML_BSSN_RHS_Body(cGH const * restrict const cctkGH, int const dir, int cons // CCTK_REAL gt11L = INITVALUE, gt11rhsL = INITVALUE, gt12L = INITVALUE, gt12rhsL = INITVALUE, gt13L = INITVALUE, gt13rhsL = INITVALUE; // CCTK_REAL gt22L = INITVALUE, gt22rhsL = INITVALUE, gt23L = INITVALUE, gt23rhsL = INITVALUE, gt33L = INITVALUE, gt33rhsL = INITVALUE; // CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE; + // CCTK_REAL rL = INITVALUE; // CCTK_REAL trKL = INITVALUE, trKrhsL = INITVALUE; // CCTK_REAL Xt1L = INITVALUE, Xt1rhsL = INITVALUE, Xt2L = INITVALUE, Xt2rhsL = INITVALUE, Xt3L = INITVALUE, Xt3rhsL = INITVALUE; /* Declare precomputed derivatives*/ @@ -290,7 +291,6 @@ void ML_BSSN_RHS_Body(cGH const * restrict const cctkGH, int const dir, int cons CCTK_REAL const beta1L = beta1[index]; CCTK_REAL const beta2L = beta2[index]; CCTK_REAL const beta3L = beta3[index]; - CCTK_REAL const etaL = eta[index]; CCTK_REAL const eTttL = (*stress_energy_state) ? (eTtt[index]) : 0.0; CCTK_REAL const eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0; CCTK_REAL const eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0; @@ -308,6 +308,7 @@ void ML_BSSN_RHS_Body(cGH const * restrict const cctkGH, int const dir, int cons CCTK_REAL const gt23L = gt23[index]; CCTK_REAL const gt33L = gt33[index]; CCTK_REAL const phiL = phi[index]; + CCTK_REAL const rL = r[index]; CCTK_REAL const trKL = trK[index]; CCTK_REAL const Xt1L = Xt1[index]; CCTK_REAL const Xt2L = Xt2[index]; @@ -1088,6 +1089,8 @@ void ML_BSSN_RHS_Body(cGH const * restrict const cctkGH, int const dir, int cons CCTK_REAL const ArhsL = (-1 + LapseAdvectionCoeff)*(AL*AlphaDriver - trKrhsL); + CCTK_REAL const eta = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); + CCTK_REAL const beta1rhsL = (PDupwindNth1(beta1, i, j, k)*beta1L + PDupwindNth2(beta1, i, j, k)*beta2L + PDupwindNth3(beta1, i, j, k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff; @@ -1097,15 +1100,15 @@ void ML_BSSN_RHS_Body(cGH const * restrict const cctkGH, int const dir, int cons CCTK_REAL const beta3rhsL = (PDupwindNth1(beta3, i, j, k)*beta1L + PDupwindNth2(beta3, i, j, k)*beta2L + PDupwindNth3(beta3, i, j, k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff; - CCTK_REAL const B1rhsL = -(B1L*etaL) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + + CCTK_REAL const B1rhsL = -(B1L*eta) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*beta2L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt1rhsL; - CCTK_REAL const B2rhsL = -(B2L*etaL) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + + CCTK_REAL const B2rhsL = -(B2L*eta) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*beta2L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt2rhsL; - CCTK_REAL const B3rhsL = -(B3L*etaL) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + + CCTK_REAL const B3rhsL = -(B3L*eta) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + (PDupwindNth2(B3, i, j, k) - PDupwindNth2(Xt3, i, j, k))*beta2L + (PDupwindNth3(B3, i, j, k) - PDupwindNth3(Xt3, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt3rhsL; diff --git a/ML_BSSN/src/ML_BSSN_RHS1.c b/ML_BSSN/src/ML_BSSN_RHS1.c index d0ded82..23f9856 100644 --- a/ML_BSSN/src/ML_BSSN_RHS1.c +++ b/ML_BSSN/src/ML_BSSN_RHS1.c @@ -105,6 +105,7 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con // CCTK_REAL cdphi1 = INITVALUE, cdphi2 = INITVALUE, cdphi3 = INITVALUE; // CCTK_REAL detgt = INITVALUE; // CCTK_REAL dir1 = INITVALUE, dir2 = INITVALUE, dir3 = INITVALUE; + // CCTK_REAL eta = INITVALUE; // CCTK_REAL fac1 = INITVALUE; // CCTK_REAL Gt111 = INITVALUE, Gt112 = INITVALUE, Gt113 = INITVALUE, Gt122 = INITVALUE, Gt123 = INITVALUE, Gt133 = INITVALUE; // CCTK_REAL Gt211 = INITVALUE, Gt212 = INITVALUE, Gt213 = INITVALUE, Gt222 = INITVALUE, Gt223 = INITVALUE, Gt233 = INITVALUE; @@ -118,7 +119,6 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con // CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE; // CCTK_REAL B1L = INITVALUE, B1rhsL = INITVALUE, B2L = INITVALUE, B2rhsL = INITVALUE, B3L = INITVALUE, B3rhsL = INITVALUE; // CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE; - // CCTK_REAL etaL = INITVALUE; // CCTK_REAL eTtxL = INITVALUE; // CCTK_REAL eTtyL = INITVALUE; // CCTK_REAL eTtzL = INITVALUE; @@ -131,6 +131,7 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con // CCTK_REAL gt11L = INITVALUE, gt11rhsL = INITVALUE, gt12L = INITVALUE, gt12rhsL = INITVALUE, gt13L = INITVALUE, gt13rhsL = INITVALUE; // CCTK_REAL gt22L = INITVALUE, gt22rhsL = INITVALUE, gt23L = INITVALUE, gt23rhsL = INITVALUE, gt33L = INITVALUE, gt33rhsL = INITVALUE; // CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE; + // CCTK_REAL rL = INITVALUE; // CCTK_REAL trKL = INITVALUE; // CCTK_REAL Xt1L = INITVALUE, Xt1rhsL = INITVALUE, Xt2L = INITVALUE, Xt2rhsL = INITVALUE, Xt3L = INITVALUE, Xt3rhsL = INITVALUE; /* Declare precomputed derivatives*/ @@ -205,7 +206,6 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con CCTK_REAL const beta1L = beta1[index]; CCTK_REAL const beta2L = beta2[index]; CCTK_REAL const beta3L = beta3[index]; - CCTK_REAL const etaL = eta[index]; CCTK_REAL const eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0; CCTK_REAL const eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0; CCTK_REAL const eTtzL = (*stress_energy_state) ? (eTtz[index]) : 0.0; @@ -222,6 +222,7 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con CCTK_REAL const gt23L = gt23[index]; CCTK_REAL const gt33L = gt33[index]; CCTK_REAL const phiL = phi[index]; + CCTK_REAL const rL = r[index]; CCTK_REAL const trKL = trK[index]; CCTK_REAL const Xt1L = Xt1[index]; CCTK_REAL const Xt2L = Xt2[index]; @@ -486,6 +487,8 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con 2*(PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3)*Xtn3 - 3*(PDstandardNth1beta3*Xtn1 + PDstandardNth2beta3*Xtn2 + PDstandardNth3beta3*Xtn3)); + CCTK_REAL const eta = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); + CCTK_REAL const beta1rhsL = (PDupwindNth1(beta1, i, j, k)*beta1L + PDupwindNth2(beta1, i, j, k)*beta2L + PDupwindNth3(beta1, i, j, k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff; @@ -495,15 +498,15 @@ void ML_BSSN_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int con CCTK_REAL const beta3rhsL = (PDupwindNth1(beta3, i, j, k)*beta1L + PDupwindNth2(beta3, i, j, k)*beta2L + PDupwindNth3(beta3, i, j, k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff; - CCTK_REAL const B1rhsL = -(B1L*etaL) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + + CCTK_REAL const B1rhsL = -(B1L*eta) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*beta2L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt1rhsL; - CCTK_REAL const B2rhsL = -(B2L*etaL) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + + CCTK_REAL const B2rhsL = -(B2L*eta) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*beta2L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt2rhsL; - CCTK_REAL const B3rhsL = -(B3L*etaL) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + + CCTK_REAL const B3rhsL = -(B3L*eta) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + (PDupwindNth2(B3, i, j, k) - PDupwindNth2(Xt3, i, j, k))*beta2L + (PDupwindNth3(B3, i, j, k) - PDupwindNth3(Xt3, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt3rhsL; diff --git a/ML_BSSN/src/ML_BSSN_setBetaDriverConstant.c b/ML_BSSN/src/ML_BSSN_setBetaDriverConstant.c deleted file mode 100644 index 18721ad..0000000 --- a/ML_BSSN/src/ML_BSSN_setBetaDriverConstant.c +++ /dev/null @@ -1,137 +0,0 @@ -/* File produced by Kranc */ - -#define KRANC_C - -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.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_setBetaDriverConstant_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - /* Declare finite differencing variables */ - - /* 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 p1odx = INITVALUE; - // CCTK_REAL p1ody = INITVALUE; - // CCTK_REAL p1odz = INITVALUE; - // CCTK_REAL pm1o12dx2 = INITVALUE; - // CCTK_REAL pm1o12dy2 = INITVALUE; - // CCTK_REAL pm1o12dz2 = INITVALUE; - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_setBetaDriverConstant_Body"); - } - - if (cctk_iteration % ML_BSSN_setBetaDriverConstant_calc_every != ML_BSSN_setBetaDriverConstant_calc_offset) - { - return; - } - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - CCTK_REAL const dx = CCTK_DELTA_SPACE(0); - CCTK_REAL const dy = CCTK_DELTA_SPACE(1); - CCTK_REAL const dz = CCTK_DELTA_SPACE(2); - int const di = 1; - int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); - int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); - CCTK_REAL const dxi = 1.0 / dx; - CCTK_REAL const dyi = 1.0 / dy; - CCTK_REAL const dzi = 1.0 / dz; - CCTK_REAL const khalf = 0.5; - CCTK_REAL const kthird = 1/3.0; - CCTK_REAL const ktwothird = 2.0/3.0; - CCTK_REAL const kfourthird = 4.0/3.0; - CCTK_REAL const keightthird = 8.0/3.0; - CCTK_REAL const hdxi = 0.5 * dxi; - CCTK_REAL const hdyi = 0.5 * dyi; - CCTK_REAL const hdzi = 0.5 * dzi; - - /* Initialize predefined quantities */ - CCTK_REAL const p1o12dx = INV(dx)/12.; - CCTK_REAL const p1o12dy = INV(dy)/12.; - CCTK_REAL const p1o12dz = INV(dz)/12.; - CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; - CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; - CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; - CCTK_REAL const p1odx = INV(dx); - CCTK_REAL const p1ody = INV(dy); - CCTK_REAL const p1odz = INV(dz); - CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; - CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; - CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; - - /* Loop over the grid points */ - #pragma omp parallel - LC_LOOP3 (ML_BSSN_setBetaDriverConstant, - 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; - int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); - int const 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 etaL = 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 */ - CCTK_REAL const etaL = BetaDriver; - - - /* Copy local copies back to grid functions */ - eta[index] = etaL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (ML_BSSN_setBetaDriverConstant); -} - -void ML_BSSN_setBetaDriverConstant(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - GenericFD_LoopOverEverything(cctkGH, &ML_BSSN_setBetaDriverConstant_Body); -} diff --git a/ML_BSSN/src/ML_BSSN_setBetaDriverSpatial.c b/ML_BSSN/src/ML_BSSN_setBetaDriverSpatial.c deleted file mode 100644 index 1e726e2..0000000 --- a/ML_BSSN/src/ML_BSSN_setBetaDriverSpatial.c +++ /dev/null @@ -1,139 +0,0 @@ -/* File produced by Kranc */ - -#define KRANC_C - -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.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_setBetaDriverSpatial_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - /* Declare finite differencing variables */ - - /* 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 p1odx = INITVALUE; - // CCTK_REAL p1ody = INITVALUE; - // CCTK_REAL p1odz = INITVALUE; - // CCTK_REAL pm1o12dx2 = INITVALUE; - // CCTK_REAL pm1o12dy2 = INITVALUE; - // CCTK_REAL pm1o12dz2 = INITVALUE; - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_setBetaDriverSpatial_Body"); - } - - if (cctk_iteration % ML_BSSN_setBetaDriverSpatial_calc_every != ML_BSSN_setBetaDriverSpatial_calc_offset) - { - return; - } - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - CCTK_REAL const dx = CCTK_DELTA_SPACE(0); - CCTK_REAL const dy = CCTK_DELTA_SPACE(1); - CCTK_REAL const dz = CCTK_DELTA_SPACE(2); - int const di = 1; - int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); - int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); - CCTK_REAL const dxi = 1.0 / dx; - CCTK_REAL const dyi = 1.0 / dy; - CCTK_REAL const dzi = 1.0 / dz; - CCTK_REAL const khalf = 0.5; - CCTK_REAL const kthird = 1/3.0; - CCTK_REAL const ktwothird = 2.0/3.0; - CCTK_REAL const kfourthird = 4.0/3.0; - CCTK_REAL const keightthird = 8.0/3.0; - CCTK_REAL const hdxi = 0.5 * dxi; - CCTK_REAL const hdyi = 0.5 * dyi; - CCTK_REAL const hdzi = 0.5 * dzi; - - /* Initialize predefined quantities */ - CCTK_REAL const p1o12dx = INV(dx)/12.; - CCTK_REAL const p1o12dy = INV(dy)/12.; - CCTK_REAL const p1o12dz = INV(dz)/12.; - CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; - CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; - CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; - CCTK_REAL const p1odx = INV(dx); - CCTK_REAL const p1ody = INV(dy); - CCTK_REAL const p1odz = INV(dz); - CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; - CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; - CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; - - /* Loop over the grid points */ - #pragma omp parallel - LC_LOOP3 (ML_BSSN_setBetaDriverSpatial, - 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; - int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); - int const 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 etaL = INITVALUE; - // CCTK_REAL rL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - CCTK_REAL const rL = r[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 */ - CCTK_REAL const etaL = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); - - - /* Copy local copies back to grid functions */ - eta[index] = etaL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (ML_BSSN_setBetaDriverSpatial); -} - -void ML_BSSN_setBetaDriverSpatial(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - GenericFD_LoopOverEverything(cctkGH, &ML_BSSN_setBetaDriverSpatial_Body); -} diff --git a/ML_BSSN/src/RegisterSymmetries.c b/ML_BSSN/src/RegisterSymmetries.c index 16bd0d9..20db895 100644 --- a/ML_BSSN/src/RegisterSymmetries.c +++ b/ML_BSSN/src/RegisterSymmetries.c @@ -144,11 +144,6 @@ void ML_BSSN_RegisterSymmetries(CCTK_ARGUMENTS) sym[0] = 1; sym[1] = 1; sym[2] = 1; - SetCartSymVN(cctkGH, sym, "ML_BSSN::eta"); - - sym[0] = 1; - sym[1] = 1; - sym[2] = 1; SetCartSymVN(cctkGH, sym, "ML_BSSN::cS"); sym[0] = -1; diff --git a/ML_BSSN/src/make.code.defn b/ML_BSSN/src/make.code.defn index eed6fb1..e78720b 100644 --- a/ML_BSSN/src/make.code.defn +++ b/ML_BSSN/src/make.code.defn @@ -1,3 +1,3 @@ # File produced by Kranc -SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_Minkowski.c ML_BSSN_convertFromADMBase.c ML_BSSN_convertFromADMBaseGamma.c ML_BSSN_setBetaDriverConstant.c ML_BSSN_setBetaDriverSpatial.c ML_BSSN_RHS.c ML_BSSN_RHS1.c ML_BSSN_RHS2.c ML_BSSN_RHSStaticBoundary.c ML_BSSN_RHSRadiativeBoundary.c ML_BSSN_enforce.c ML_BSSN_enforce2.c ML_BSSN_boundary.c ML_BSSN_convertToADMBase.c ML_BSSN_convertToADMBaseDtLapseShift.c ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_convertToADMBaseFakeDtLapseShift.c ML_BSSN_constraints.c ML_BSSN_constraints_boundary.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_RHS1.c ML_BSSN_RHS2.c ML_BSSN_RHSStaticBoundary.c ML_BSSN_RHSRadiativeBoundary.c ML_BSSN_enforce.c ML_BSSN_enforce2.c ML_BSSN_boundary.c ML_BSSN_convertToADMBase.c ML_BSSN_convertToADMBaseDtLapseShift.c ML_BSSN_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_convertToADMBaseFakeDtLapseShift.c ML_BSSN_constraints.c ML_BSSN_constraints_boundary.c Boundaries.c diff --git a/ML_BSSN_Helper/src/SetGroupTags.c b/ML_BSSN_Helper/src/SetGroupTags.c index 3113637..0666c1b 100644 --- a/ML_BSSN_Helper/src/SetGroupTags.c +++ b/ML_BSSN_Helper/src/SetGroupTags.c @@ -28,8 +28,6 @@ ML_BSSN_SetGroupTags (void) set_group_tags (0, 0, 0, "ML_BSSN::ML_Ham"); set_group_tags (0, 0, 0, "ML_BSSN::ML_mom"); - set_group_tags (0, 1, 0, "ML_BSSN::ML_BetaDriver"); - int const checkpoint = rhs_timelevels > 1; set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_dtlapserhs"); set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_dtshiftrhs"); diff --git a/ML_BSSN_MP/interface.ccl b/ML_BSSN_MP/interface.ccl index d284d84..2291b0c 100644 --- a/ML_BSSN_MP/interface.ccl +++ b/ML_BSSN_MP/interface.ccl @@ -25,12 +25,6 @@ CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT USES FUNCTION Boundary_SelectVarForBC public: -CCTK_REAL ML_BetaDriver type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' -{ - eta -} "ML_BetaDriver" - -public: CCTK_REAL ML_cons_detg type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=2.0000000000000000000' { cS diff --git a/ML_BSSN_MP/param.ccl b/ML_BSSN_MP/param.ccl index c62baf4..99680bb 100644 --- a/ML_BSSN_MP/param.ccl +++ b/ML_BSSN_MP/param.ccl @@ -207,18 +207,6 @@ CCTK_INT ML_BSSN_MP_convertFromADMBaseGamma_calc_every "ML_BSSN_MP_convertFromAD } 1 restricted: -CCTK_INT ML_BSSN_MP_setBetaDriverConstant_calc_every "ML_BSSN_MP_setBetaDriverConstant_calc_every" -{ - *:* :: "" -} 1 - -restricted: -CCTK_INT ML_BSSN_MP_setBetaDriverSpatial_calc_every "ML_BSSN_MP_setBetaDriverSpatial_calc_every" -{ - *:* :: "" -} 1 - -restricted: CCTK_INT ML_BSSN_MP_RHS_calc_every "ML_BSSN_MP_RHS_calc_every" { *:* :: "" @@ -321,18 +309,6 @@ CCTK_INT ML_BSSN_MP_convertFromADMBaseGamma_calc_offset "ML_BSSN_MP_convertFromA } 0 restricted: -CCTK_INT ML_BSSN_MP_setBetaDriverConstant_calc_offset "ML_BSSN_MP_setBetaDriverConstant_calc_offset" -{ - *:* :: "" -} 0 - -restricted: -CCTK_INT ML_BSSN_MP_setBetaDriverSpatial_calc_offset "ML_BSSN_MP_setBetaDriverSpatial_calc_offset" -{ - *:* :: "" -} 0 - -restricted: CCTK_INT ML_BSSN_MP_RHS_calc_offset "ML_BSSN_MP_RHS_calc_offset" { *:* :: "" diff --git a/ML_BSSN_MP/schedule.ccl b/ML_BSSN_MP/schedule.ccl index e5cf2ae..17fb41c 100644 --- a/ML_BSSN_MP/schedule.ccl +++ b/ML_BSSN_MP/schedule.ccl @@ -1,8 +1,6 @@ # File produced by Kranc -STORAGE: ML_BetaDriver[1] - STORAGE: ML_cons_detg[1] STORAGE: ML_cons_Gamma[1] @@ -367,24 +365,6 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase")) } "ML_BSSN_MP_convertFromADMBaseGamma" } - -if (CCTK_EQUALS(UseSpatialBetaDriver, "no")) -{ - schedule ML_BSSN_MP_setBetaDriverConstant IN ML_BSSN_MP_InitEta - { - LANG: C - } "ML_BSSN_MP_setBetaDriverConstant" -} - - -if (CCTK_EQUALS(UseSpatialBetaDriver, "yes")) -{ - schedule ML_BSSN_MP_setBetaDriverSpatial IN ML_BSSN_MP_InitEta - { - LANG: C - } "ML_BSSN_MP_setBetaDriverSpatial" -} - schedule ML_BSSN_MP_RHS IN NoSuchGroup { LANG: C diff --git a/ML_BSSN_MP/src/ML_BSSN_MP_RHS.c b/ML_BSSN_MP/src/ML_BSSN_MP_RHS.c index 272a573..b45d1c3 100644 --- a/ML_BSSN_MP/src/ML_BSSN_MP_RHS.c +++ b/ML_BSSN_MP/src/ML_BSSN_MP_RHS.c @@ -109,6 +109,7 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c // CCTK_REAL dir1 = INITVALUE, dir2 = INITVALUE, dir3 = INITVALUE; // CCTK_REAL e4phi = INITVALUE; // CCTK_REAL em4phi = INITVALUE; + // CCTK_REAL eta = INITVALUE; // CCTK_REAL fac1 = INITVALUE, fac2 = INITVALUE; // CCTK_REAL g11 = INITVALUE; // CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE; @@ -147,7 +148,6 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c // CCTK_REAL dJ111L = INITVALUE, dJ112L = INITVALUE, dJ113L = INITVALUE, dJ122L = INITVALUE, dJ123L = INITVALUE, dJ133L = INITVALUE; // CCTK_REAL dJ211L = INITVALUE, dJ212L = INITVALUE, dJ213L = INITVALUE, dJ222L = INITVALUE, dJ223L = INITVALUE, dJ233L = INITVALUE; // CCTK_REAL dJ311L = INITVALUE, dJ312L = INITVALUE, dJ313L = INITVALUE, dJ322L = INITVALUE, dJ323L = INITVALUE, dJ333L = INITVALUE; - // CCTK_REAL etaL = INITVALUE; // CCTK_REAL eTttL = INITVALUE; // CCTK_REAL eTtxL = INITVALUE; // CCTK_REAL eTtyL = INITVALUE; @@ -163,6 +163,7 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c // CCTK_REAL J11L = INITVALUE, J12L = INITVALUE, J13L = INITVALUE, J21L = INITVALUE, J22L = INITVALUE, J23L = INITVALUE; // CCTK_REAL J31L = INITVALUE, J32L = INITVALUE, J33L = INITVALUE; // CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE; + // CCTK_REAL rL = INITVALUE; // CCTK_REAL trKL = INITVALUE, trKrhsL = INITVALUE; // CCTK_REAL Xt1L = INITVALUE, Xt1rhsL = INITVALUE, Xt2L = INITVALUE, Xt2rhsL = INITVALUE, Xt3L = INITVALUE, Xt3rhsL = INITVALUE; /* Declare precomputed derivatives*/ @@ -313,7 +314,6 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c CCTK_REAL const dJ322L = dJ322[index]; CCTK_REAL const dJ323L = dJ323[index]; CCTK_REAL const dJ333L = dJ333[index]; - CCTK_REAL const etaL = eta[index]; CCTK_REAL const eTttL = (*stress_energy_state) ? (eTtt[index]) : 0.0; CCTK_REAL const eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0; CCTK_REAL const eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0; @@ -340,6 +340,7 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c CCTK_REAL const J32L = J32[index]; CCTK_REAL const J33L = J33[index]; CCTK_REAL const phiL = phi[index]; + CCTK_REAL const rL = r[index]; CCTK_REAL const trKL = trK[index]; CCTK_REAL const Xt1L = Xt1[index]; CCTK_REAL const Xt2L = Xt2[index]; @@ -1640,6 +1641,8 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c CCTK_REAL const ArhsL = (-1 + LapseAdvectionCoeff)*(AL*AlphaDriver - trKrhsL); + CCTK_REAL const eta = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); + CCTK_REAL const beta1rhsL = (PDupwindNth1(beta1, i, j, k)*(beta1L*J11L + beta2L*J12L + beta3L*J13L) + PDupwindNth2(beta1, i, j, k)*(beta1L*J21L + beta2L*J22L + beta3L*J23L) + PDupwindNth3(beta1, i, j, k)*(beta1L*J31L + beta2L*J32L + beta3L*J33L))*ShiftAdvectionCoeff + @@ -1655,7 +1658,7 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c PDupwindNth3(beta3, i, j, k)*(beta1L*J31L + beta2L*J32L + beta3L*J33L))*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff; - CCTK_REAL const B1rhsL = -(B1L*etaL) + (beta1L*((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*J11L + + CCTK_REAL const B1rhsL = -(B1L*eta) + (beta1L*((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*J11L + (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*J21L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*J31L) + beta2L*((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*J12L + @@ -1665,7 +1668,7 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*J23L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*J33L))*ShiftAdvectionCoeff + Xt1rhsL; - CCTK_REAL const B2rhsL = -(B2L*etaL) + (beta1L*((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*J11L + + CCTK_REAL const B2rhsL = -(B2L*eta) + (beta1L*((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*J11L + (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*J21L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*J31L) + beta2L*((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*J12L + @@ -1675,7 +1678,7 @@ void ML_BSSN_MP_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*J23L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*J33L))*ShiftAdvectionCoeff + Xt2rhsL; - CCTK_REAL const B3rhsL = -(B3L*etaL) + (beta1L*((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*J11L + + CCTK_REAL const B3rhsL = -(B3L*eta) + (beta1L*((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*J11L + (PDupwindNth2(B3, i, j, k) - PDupwindNth2(Xt3, i, j, k))*J21L + (PDupwindNth3(B3, i, j, k) - PDupwindNth3(Xt3, i, j, k))*J31L) + beta2L*((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*J12L + diff --git a/ML_BSSN_MP/src/ML_BSSN_MP_RHS1.c b/ML_BSSN_MP/src/ML_BSSN_MP_RHS1.c index 3207d05..0e48762 100644 --- a/ML_BSSN_MP/src/ML_BSSN_MP_RHS1.c +++ b/ML_BSSN_MP/src/ML_BSSN_MP_RHS1.c @@ -105,6 +105,7 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int // CCTK_REAL cdphi1 = INITVALUE, cdphi2 = INITVALUE, cdphi3 = INITVALUE; // CCTK_REAL detgt = INITVALUE; // CCTK_REAL dir1 = INITVALUE, dir2 = INITVALUE, dir3 = INITVALUE; + // CCTK_REAL eta = INITVALUE; // CCTK_REAL fac1 = INITVALUE; // CCTK_REAL Gt111 = INITVALUE, Gt112 = INITVALUE, Gt113 = INITVALUE, Gt122 = INITVALUE, Gt123 = INITVALUE, Gt133 = INITVALUE; // CCTK_REAL Gt211 = INITVALUE, Gt212 = INITVALUE, Gt213 = INITVALUE, Gt222 = INITVALUE, Gt223 = INITVALUE, Gt233 = INITVALUE; @@ -121,7 +122,6 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int // CCTK_REAL dJ111L = INITVALUE, dJ112L = INITVALUE, dJ113L = INITVALUE, dJ122L = INITVALUE, dJ123L = INITVALUE, dJ133L = INITVALUE; // CCTK_REAL dJ211L = INITVALUE, dJ212L = INITVALUE, dJ213L = INITVALUE, dJ222L = INITVALUE, dJ223L = INITVALUE, dJ233L = INITVALUE; // CCTK_REAL dJ311L = INITVALUE, dJ312L = INITVALUE, dJ313L = INITVALUE, dJ322L = INITVALUE, dJ323L = INITVALUE, dJ333L = INITVALUE; - // CCTK_REAL etaL = INITVALUE; // CCTK_REAL eTtxL = INITVALUE; // CCTK_REAL eTtyL = INITVALUE; // CCTK_REAL eTtzL = INITVALUE; @@ -136,6 +136,7 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int // CCTK_REAL J11L = INITVALUE, J12L = INITVALUE, J13L = INITVALUE, J21L = INITVALUE, J22L = INITVALUE, J23L = INITVALUE; // CCTK_REAL J31L = INITVALUE, J32L = INITVALUE, J33L = INITVALUE; // CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE; + // CCTK_REAL rL = INITVALUE; // CCTK_REAL trKL = INITVALUE; // CCTK_REAL Xt1L = INITVALUE, Xt1rhsL = INITVALUE, Xt2L = INITVALUE, Xt2rhsL = INITVALUE, Xt3L = INITVALUE, Xt3rhsL = INITVALUE; /* Declare precomputed derivatives*/ @@ -228,7 +229,6 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int CCTK_REAL const dJ322L = dJ322[index]; CCTK_REAL const dJ323L = dJ323[index]; CCTK_REAL const dJ333L = dJ333[index]; - CCTK_REAL const etaL = eta[index]; CCTK_REAL const eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0; CCTK_REAL const eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0; CCTK_REAL const eTtzL = (*stress_energy_state) ? (eTtz[index]) : 0.0; @@ -254,6 +254,7 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int CCTK_REAL const J32L = J32[index]; CCTK_REAL const J33L = J33[index]; CCTK_REAL const phiL = phi[index]; + CCTK_REAL const rL = r[index]; CCTK_REAL const trKL = trK[index]; CCTK_REAL const Xt1L = Xt1[index]; CCTK_REAL const Xt2L = Xt2[index]; @@ -785,6 +786,8 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int dJ313L*PDstandardNth3beta1 + dJ323L*PDstandardNth3beta2 + dJ333L*PDstandardNth3beta3 + PDstandardNth11beta3*SQR(J13L) + PDstandardNth22beta3*SQR(J23L) + PDstandardNth33beta3*SQR(J33L))); + CCTK_REAL const eta = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); + CCTK_REAL const beta1rhsL = (PDupwindNth1(beta1, i, j, k)*(beta1L*J11L + beta2L*J12L + beta3L*J13L) + PDupwindNth2(beta1, i, j, k)*(beta1L*J21L + beta2L*J22L + beta3L*J23L) + PDupwindNth3(beta1, i, j, k)*(beta1L*J31L + beta2L*J32L + beta3L*J33L))*ShiftAdvectionCoeff + @@ -800,7 +803,7 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int PDupwindNth3(beta3, i, j, k)*(beta1L*J31L + beta2L*J32L + beta3L*J33L))*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff; - CCTK_REAL const B1rhsL = -(B1L*etaL) + (beta1L*((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*J11L + + CCTK_REAL const B1rhsL = -(B1L*eta) + (beta1L*((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*J11L + (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*J21L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*J31L) + beta2L*((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*J12L + @@ -810,7 +813,7 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*J23L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*J33L))*ShiftAdvectionCoeff + Xt1rhsL; - CCTK_REAL const B2rhsL = -(B2L*etaL) + (beta1L*((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*J11L + + CCTK_REAL const B2rhsL = -(B2L*eta) + (beta1L*((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*J11L + (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*J21L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*J31L) + beta2L*((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*J12L + @@ -820,7 +823,7 @@ void ML_BSSN_MP_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*J23L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*J33L))*ShiftAdvectionCoeff + Xt2rhsL; - CCTK_REAL const B3rhsL = -(B3L*etaL) + (beta1L*((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*J11L + + CCTK_REAL const B3rhsL = -(B3L*eta) + (beta1L*((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*J11L + (PDupwindNth2(B3, i, j, k) - PDupwindNth2(Xt3, i, j, k))*J21L + (PDupwindNth3(B3, i, j, k) - PDupwindNth3(Xt3, i, j, k))*J31L) + beta2L*((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*J12L + diff --git a/ML_BSSN_MP/src/ML_BSSN_MP_setBetaDriverConstant.c b/ML_BSSN_MP/src/ML_BSSN_MP_setBetaDriverConstant.c deleted file mode 100644 index 99a3e11..0000000 --- a/ML_BSSN_MP/src/ML_BSSN_MP_setBetaDriverConstant.c +++ /dev/null @@ -1,137 +0,0 @@ -/* File produced by Kranc */ - -#define KRANC_C - -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.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_MP_setBetaDriverConstant_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - /* Declare finite differencing variables */ - - /* 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 p1odx = INITVALUE; - // CCTK_REAL p1ody = INITVALUE; - // CCTK_REAL p1odz = INITVALUE; - // CCTK_REAL pm1o12dx2 = INITVALUE; - // CCTK_REAL pm1o12dy2 = INITVALUE; - // CCTK_REAL pm1o12dz2 = INITVALUE; - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_MP_setBetaDriverConstant_Body"); - } - - if (cctk_iteration % ML_BSSN_MP_setBetaDriverConstant_calc_every != ML_BSSN_MP_setBetaDriverConstant_calc_offset) - { - return; - } - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - CCTK_REAL const dx = CCTK_DELTA_SPACE(0); - CCTK_REAL const dy = CCTK_DELTA_SPACE(1); - CCTK_REAL const dz = CCTK_DELTA_SPACE(2); - int const di = 1; - int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); - int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); - CCTK_REAL const dxi = 1.0 / dx; - CCTK_REAL const dyi = 1.0 / dy; - CCTK_REAL const dzi = 1.0 / dz; - CCTK_REAL const khalf = 0.5; - CCTK_REAL const kthird = 1/3.0; - CCTK_REAL const ktwothird = 2.0/3.0; - CCTK_REAL const kfourthird = 4.0/3.0; - CCTK_REAL const keightthird = 8.0/3.0; - CCTK_REAL const hdxi = 0.5 * dxi; - CCTK_REAL const hdyi = 0.5 * dyi; - CCTK_REAL const hdzi = 0.5 * dzi; - - /* Initialize predefined quantities */ - CCTK_REAL const p1o12dx = INV(dx)/12.; - CCTK_REAL const p1o12dy = INV(dy)/12.; - CCTK_REAL const p1o12dz = INV(dz)/12.; - CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; - CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; - CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; - CCTK_REAL const p1odx = INV(dx); - CCTK_REAL const p1ody = INV(dy); - CCTK_REAL const p1odz = INV(dz); - CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; - CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; - CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; - - /* Loop over the grid points */ - #pragma omp parallel - LC_LOOP3 (ML_BSSN_MP_setBetaDriverConstant, - 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; - int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); - int const 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 etaL = 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 */ - CCTK_REAL const etaL = BetaDriver; - - - /* Copy local copies back to grid functions */ - eta[index] = etaL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (ML_BSSN_MP_setBetaDriverConstant); -} - -void ML_BSSN_MP_setBetaDriverConstant(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - GenericFD_LoopOverEverything(cctkGH, &ML_BSSN_MP_setBetaDriverConstant_Body); -} diff --git a/ML_BSSN_MP/src/ML_BSSN_MP_setBetaDriverSpatial.c b/ML_BSSN_MP/src/ML_BSSN_MP_setBetaDriverSpatial.c deleted file mode 100644 index d517184..0000000 --- a/ML_BSSN_MP/src/ML_BSSN_MP_setBetaDriverSpatial.c +++ /dev/null @@ -1,139 +0,0 @@ -/* File produced by Kranc */ - -#define KRANC_C - -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.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_MP_setBetaDriverSpatial_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - /* Declare finite differencing variables */ - - /* 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 p1odx = INITVALUE; - // CCTK_REAL p1ody = INITVALUE; - // CCTK_REAL p1odz = INITVALUE; - // CCTK_REAL pm1o12dx2 = INITVALUE; - // CCTK_REAL pm1o12dy2 = INITVALUE; - // CCTK_REAL pm1o12dz2 = INITVALUE; - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_MP_setBetaDriverSpatial_Body"); - } - - if (cctk_iteration % ML_BSSN_MP_setBetaDriverSpatial_calc_every != ML_BSSN_MP_setBetaDriverSpatial_calc_offset) - { - return; - } - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - CCTK_REAL const dx = CCTK_DELTA_SPACE(0); - CCTK_REAL const dy = CCTK_DELTA_SPACE(1); - CCTK_REAL const dz = CCTK_DELTA_SPACE(2); - int const di = 1; - int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); - int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); - CCTK_REAL const dxi = 1.0 / dx; - CCTK_REAL const dyi = 1.0 / dy; - CCTK_REAL const dzi = 1.0 / dz; - CCTK_REAL const khalf = 0.5; - CCTK_REAL const kthird = 1/3.0; - CCTK_REAL const ktwothird = 2.0/3.0; - CCTK_REAL const kfourthird = 4.0/3.0; - CCTK_REAL const keightthird = 8.0/3.0; - CCTK_REAL const hdxi = 0.5 * dxi; - CCTK_REAL const hdyi = 0.5 * dyi; - CCTK_REAL const hdzi = 0.5 * dzi; - - /* Initialize predefined quantities */ - CCTK_REAL const p1o12dx = INV(dx)/12.; - CCTK_REAL const p1o12dy = INV(dy)/12.; - CCTK_REAL const p1o12dz = INV(dz)/12.; - CCTK_REAL const p1o144dxdy = (INV(dx)*INV(dy))/144.; - CCTK_REAL const p1o144dxdz = (INV(dx)*INV(dz))/144.; - CCTK_REAL const p1o144dydz = (INV(dy)*INV(dz))/144.; - CCTK_REAL const p1odx = INV(dx); - CCTK_REAL const p1ody = INV(dy); - CCTK_REAL const p1odz = INV(dz); - CCTK_REAL const pm1o12dx2 = -pow(dx,-2)/12.; - CCTK_REAL const pm1o12dy2 = -pow(dy,-2)/12.; - CCTK_REAL const pm1o12dz2 = -pow(dz,-2)/12.; - - /* Loop over the grid points */ - #pragma omp parallel - LC_LOOP3 (ML_BSSN_MP_setBetaDriverSpatial, - 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; - int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); - int const 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 etaL = INITVALUE; - // CCTK_REAL rL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - CCTK_REAL const rL = r[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 */ - CCTK_REAL const etaL = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); - - - /* Copy local copies back to grid functions */ - eta[index] = etaL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (ML_BSSN_MP_setBetaDriverSpatial); -} - -void ML_BSSN_MP_setBetaDriverSpatial(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - GenericFD_LoopOverEverything(cctkGH, &ML_BSSN_MP_setBetaDriverSpatial_Body); -} diff --git a/ML_BSSN_MP/src/RegisterSymmetries.c b/ML_BSSN_MP/src/RegisterSymmetries.c index d6d178b..e2081a6 100644 --- a/ML_BSSN_MP/src/RegisterSymmetries.c +++ b/ML_BSSN_MP/src/RegisterSymmetries.c @@ -144,11 +144,6 @@ void ML_BSSN_MP_RegisterSymmetries(CCTK_ARGUMENTS) sym[0] = 1; sym[1] = 1; sym[2] = 1; - SetCartSymVN(cctkGH, sym, "ML_BSSN_MP::eta"); - - sym[0] = 1; - sym[1] = 1; - sym[2] = 1; SetCartSymVN(cctkGH, sym, "ML_BSSN_MP::cS"); sym[0] = -1; diff --git a/ML_BSSN_MP/src/make.code.defn b/ML_BSSN_MP/src/make.code.defn index 792fefe..070fee3 100644 --- a/ML_BSSN_MP/src/make.code.defn +++ b/ML_BSSN_MP/src/make.code.defn @@ -1,3 +1,3 @@ # File produced by Kranc -SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_MP_Minkowski.c ML_BSSN_MP_convertFromADMBase.c ML_BSSN_MP_convertFromADMBaseGamma.c ML_BSSN_MP_setBetaDriverConstant.c ML_BSSN_MP_setBetaDriverSpatial.c ML_BSSN_MP_RHS.c ML_BSSN_MP_RHS1.c ML_BSSN_MP_RHS2.c ML_BSSN_MP_RHSStaticBoundary.c ML_BSSN_MP_RHSRadiativeBoundary.c ML_BSSN_MP_enforce.c ML_BSSN_MP_enforce2.c ML_BSSN_MP_boundary.c ML_BSSN_MP_convertToADMBase.c ML_BSSN_MP_convertToADMBaseDtLapseShift.c ML_BSSN_MP_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_MP_convertToADMBaseFakeDtLapseShift.c ML_BSSN_MP_constraints.c ML_BSSN_MP_constraints_boundary.c Boundaries.c +SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_MP_Minkowski.c ML_BSSN_MP_convertFromADMBase.c ML_BSSN_MP_convertFromADMBaseGamma.c ML_BSSN_MP_RHS.c ML_BSSN_MP_RHS1.c ML_BSSN_MP_RHS2.c ML_BSSN_MP_RHSStaticBoundary.c ML_BSSN_MP_RHSRadiativeBoundary.c ML_BSSN_MP_enforce.c ML_BSSN_MP_enforce2.c ML_BSSN_MP_boundary.c ML_BSSN_MP_convertToADMBase.c ML_BSSN_MP_convertToADMBaseDtLapseShift.c ML_BSSN_MP_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_MP_convertToADMBaseFakeDtLapseShift.c ML_BSSN_MP_constraints.c ML_BSSN_MP_constraints_boundary.c Boundaries.c diff --git a/ML_BSSN_MP_Helper/src/SetGroupTags.c b/ML_BSSN_MP_Helper/src/SetGroupTags.c index 4cfc394..91845dd 100644 --- a/ML_BSSN_MP_Helper/src/SetGroupTags.c +++ b/ML_BSSN_MP_Helper/src/SetGroupTags.c @@ -28,8 +28,6 @@ ML_BSSN_MP_SetGroupTags (void) set_group_tags (0, 0, 0, "ML_BSSN_MP::ML_Ham"); set_group_tags (0, 0, 0, "ML_BSSN_MP::ML_mom"); - set_group_tags (0, 1, 0, "ML_BSSN_MP::ML_BetaDriver"); - int const checkpoint = rhs_timelevels > 1; set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN_MP::ML_dtlapserhs"); set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN_MP::ML_dtshiftrhs"); diff --git a/ML_BSSN_O2/interface.ccl b/ML_BSSN_O2/interface.ccl index 7a35162..e0300fc 100644 --- a/ML_BSSN_O2/interface.ccl +++ b/ML_BSSN_O2/interface.ccl @@ -25,12 +25,6 @@ CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT USES FUNCTION Boundary_SelectVarForBC public: -CCTK_REAL ML_BetaDriver type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' -{ - eta -} "ML_BetaDriver" - -public: CCTK_REAL ML_cons_detg type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=2.0000000000000000000' { cS diff --git a/ML_BSSN_O2/param.ccl b/ML_BSSN_O2/param.ccl index 95b8e34..79956f9 100644 --- a/ML_BSSN_O2/param.ccl +++ b/ML_BSSN_O2/param.ccl @@ -207,18 +207,6 @@ CCTK_INT ML_BSSN_O2_convertFromADMBaseGamma_calc_every "ML_BSSN_O2_convertFromAD } 1 restricted: -CCTK_INT ML_BSSN_O2_setBetaDriverConstant_calc_every "ML_BSSN_O2_setBetaDriverConstant_calc_every" -{ - *:* :: "" -} 1 - -restricted: -CCTK_INT ML_BSSN_O2_setBetaDriverSpatial_calc_every "ML_BSSN_O2_setBetaDriverSpatial_calc_every" -{ - *:* :: "" -} 1 - -restricted: CCTK_INT ML_BSSN_O2_RHS_calc_every "ML_BSSN_O2_RHS_calc_every" { *:* :: "" @@ -321,18 +309,6 @@ CCTK_INT ML_BSSN_O2_convertFromADMBaseGamma_calc_offset "ML_BSSN_O2_convertFromA } 0 restricted: -CCTK_INT ML_BSSN_O2_setBetaDriverConstant_calc_offset "ML_BSSN_O2_setBetaDriverConstant_calc_offset" -{ - *:* :: "" -} 0 - -restricted: -CCTK_INT ML_BSSN_O2_setBetaDriverSpatial_calc_offset "ML_BSSN_O2_setBetaDriverSpatial_calc_offset" -{ - *:* :: "" -} 0 - -restricted: CCTK_INT ML_BSSN_O2_RHS_calc_offset "ML_BSSN_O2_RHS_calc_offset" { *:* :: "" diff --git a/ML_BSSN_O2/schedule.ccl b/ML_BSSN_O2/schedule.ccl index 80e2597..3e90062 100644 --- a/ML_BSSN_O2/schedule.ccl +++ b/ML_BSSN_O2/schedule.ccl @@ -1,8 +1,6 @@ # File produced by Kranc -STORAGE: ML_BetaDriver[1] - STORAGE: ML_cons_detg[1] STORAGE: ML_cons_Gamma[1] @@ -367,24 +365,6 @@ if (CCTK_EQUALS(my_initial_data, "ADMBase")) } "ML_BSSN_O2_convertFromADMBaseGamma" } - -if (CCTK_EQUALS(UseSpatialBetaDriver, "no")) -{ - schedule ML_BSSN_O2_setBetaDriverConstant IN ML_BSSN_O2_InitEta - { - LANG: C - } "ML_BSSN_O2_setBetaDriverConstant" -} - - -if (CCTK_EQUALS(UseSpatialBetaDriver, "yes")) -{ - schedule ML_BSSN_O2_setBetaDriverSpatial IN ML_BSSN_O2_InitEta - { - LANG: C - } "ML_BSSN_O2_setBetaDriverSpatial" -} - schedule ML_BSSN_O2_RHS IN NoSuchGroup { LANG: C diff --git a/ML_BSSN_O2/src/ML_BSSN_O2_RHS.c b/ML_BSSN_O2/src/ML_BSSN_O2_RHS.c index 1ee4dac..a080f35 100644 --- a/ML_BSSN_O2/src/ML_BSSN_O2_RHS.c +++ b/ML_BSSN_O2/src/ML_BSSN_O2_RHS.c @@ -115,6 +115,7 @@ void ML_BSSN_O2_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c // CCTK_REAL dir1 = INITVALUE, dir2 = INITVALUE, dir3 = INITVALUE; // CCTK_REAL e4phi = INITVALUE; // CCTK_REAL em4phi = INITVALUE; + // CCTK_REAL eta = INITVALUE; // CCTK_REAL fac1 = INITVALUE, fac2 = INITVALUE; // CCTK_REAL g11 = INITVALUE; // CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE; @@ -150,7 +151,6 @@ void ML_BSSN_O2_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c // CCTK_REAL At22L = INITVALUE, At22rhsL = INITVALUE, At23L = INITVALUE, At23rhsL = INITVALUE, At33L = INITVALUE, At33rhsL = INITVALUE; // CCTK_REAL B1L = INITVALUE, B1rhsL = INITVALUE, B2L = INITVALUE, B2rhsL = INITVALUE, B3L = INITVALUE, B3rhsL = INITVALUE; // CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE; - // CCTK_REAL etaL = INITVALUE; // CCTK_REAL eTttL = INITVALUE; // CCTK_REAL eTtxL = INITVALUE; // CCTK_REAL eTtyL = INITVALUE; @@ -164,6 +164,7 @@ void ML_BSSN_O2_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c // CCTK_REAL gt11L = INITVALUE, gt11rhsL = INITVALUE, gt12L = INITVALUE, gt12rhsL = INITVALUE, gt13L = INITVALUE, gt13rhsL = INITVALUE; // CCTK_REAL gt22L = INITVALUE, gt22rhsL = INITVALUE, gt23L = INITVALUE, gt23rhsL = INITVALUE, gt33L = INITVALUE, gt33rhsL = INITVALUE; // CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE; + // CCTK_REAL rL = INITVALUE; // CCTK_REAL trKL = INITVALUE, trKrhsL = INITVALUE; // CCTK_REAL Xt1L = INITVALUE, Xt1rhsL = INITVALUE, Xt2L = INITVALUE, Xt2rhsL = INITVALUE, Xt3L = INITVALUE, Xt3rhsL = INITVALUE; /* Declare precomputed derivatives*/ @@ -296,7 +297,6 @@ void ML_BSSN_O2_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c CCTK_REAL const beta1L = beta1[index]; CCTK_REAL const beta2L = beta2[index]; CCTK_REAL const beta3L = beta3[index]; - CCTK_REAL const etaL = eta[index]; CCTK_REAL const eTttL = (*stress_energy_state) ? (eTtt[index]) : 0.0; CCTK_REAL const eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0; CCTK_REAL const eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0; @@ -314,6 +314,7 @@ void ML_BSSN_O2_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c CCTK_REAL const gt23L = gt23[index]; CCTK_REAL const gt33L = gt33[index]; CCTK_REAL const phiL = phi[index]; + CCTK_REAL const rL = r[index]; CCTK_REAL const trKL = trK[index]; CCTK_REAL const Xt1L = Xt1[index]; CCTK_REAL const Xt2L = Xt2[index]; @@ -1094,6 +1095,8 @@ void ML_BSSN_O2_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c CCTK_REAL const ArhsL = (-1 + LapseAdvectionCoeff)*(AL*AlphaDriver - trKrhsL); + CCTK_REAL const eta = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); + CCTK_REAL const beta1rhsL = (PDupwindNth1(beta1, i, j, k)*beta1L + PDupwindNth2(beta1, i, j, k)*beta2L + PDupwindNth3(beta1, i, j, k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff; @@ -1103,15 +1106,15 @@ void ML_BSSN_O2_RHS_Body(cGH const * restrict const cctkGH, int const dir, int c CCTK_REAL const beta3rhsL = (PDupwindNth1(beta3, i, j, k)*beta1L + PDupwindNth2(beta3, i, j, k)*beta2L + PDupwindNth3(beta3, i, j, k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff; - CCTK_REAL const B1rhsL = -(B1L*etaL) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + + CCTK_REAL const B1rhsL = -(B1L*eta) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*beta2L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt1rhsL; - CCTK_REAL const B2rhsL = -(B2L*etaL) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + + CCTK_REAL const B2rhsL = -(B2L*eta) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*beta2L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt2rhsL; - CCTK_REAL const B3rhsL = -(B3L*etaL) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + + CCTK_REAL const B3rhsL = -(B3L*eta) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + (PDupwindNth2(B3, i, j, k) - PDupwindNth2(Xt3, i, j, k))*beta2L + (PDupwindNth3(B3, i, j, k) - PDupwindNth3(Xt3, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt3rhsL; diff --git a/ML_BSSN_O2/src/ML_BSSN_O2_RHS1.c b/ML_BSSN_O2/src/ML_BSSN_O2_RHS1.c index 29cef52..0ea2889 100644 --- a/ML_BSSN_O2/src/ML_BSSN_O2_RHS1.c +++ b/ML_BSSN_O2/src/ML_BSSN_O2_RHS1.c @@ -111,6 +111,7 @@ void ML_BSSN_O2_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int // CCTK_REAL cdphi1 = INITVALUE, cdphi2 = INITVALUE, cdphi3 = INITVALUE; // CCTK_REAL detgt = INITVALUE; // CCTK_REAL dir1 = INITVALUE, dir2 = INITVALUE, dir3 = INITVALUE; + // CCTK_REAL eta = INITVALUE; // CCTK_REAL fac1 = INITVALUE; // CCTK_REAL Gt111 = INITVALUE, Gt112 = INITVALUE, Gt113 = INITVALUE, Gt122 = INITVALUE, Gt123 = INITVALUE, Gt133 = INITVALUE; // CCTK_REAL Gt211 = INITVALUE, Gt212 = INITVALUE, Gt213 = INITVALUE, Gt222 = INITVALUE, Gt223 = INITVALUE, Gt233 = INITVALUE; @@ -124,7 +125,6 @@ void ML_BSSN_O2_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int // CCTK_REAL At11L = INITVALUE, At12L = INITVALUE, At13L = INITVALUE, At22L = INITVALUE, At23L = INITVALUE, At33L = INITVALUE; // CCTK_REAL B1L = INITVALUE, B1rhsL = INITVALUE, B2L = INITVALUE, B2rhsL = INITVALUE, B3L = INITVALUE, B3rhsL = INITVALUE; // CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE; - // CCTK_REAL etaL = INITVALUE; // CCTK_REAL eTtxL = INITVALUE; // CCTK_REAL eTtyL = INITVALUE; // CCTK_REAL eTtzL = INITVALUE; @@ -137,6 +137,7 @@ void ML_BSSN_O2_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int // CCTK_REAL gt11L = INITVALUE, gt11rhsL = INITVALUE, gt12L = INITVALUE, gt12rhsL = INITVALUE, gt13L = INITVALUE, gt13rhsL = INITVALUE; // CCTK_REAL gt22L = INITVALUE, gt22rhsL = INITVALUE, gt23L = INITVALUE, gt23rhsL = INITVALUE, gt33L = INITVALUE, gt33rhsL = INITVALUE; // CCTK_REAL phiL = INITVALUE, phirhsL = INITVALUE; + // CCTK_REAL rL = INITVALUE; // CCTK_REAL trKL = INITVALUE; // CCTK_REAL Xt1L = INITVALUE, Xt1rhsL = INITVALUE, Xt2L = INITVALUE, Xt2rhsL = INITVALUE, Xt3L = INITVALUE, Xt3rhsL = INITVALUE; /* Declare precomputed derivatives*/ @@ -211,7 +212,6 @@ void ML_BSSN_O2_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int CCTK_REAL const beta1L = beta1[index]; CCTK_REAL const beta2L = beta2[index]; CCTK_REAL const beta3L = beta3[index]; - CCTK_REAL const etaL = eta[index]; CCTK_REAL const eTtxL = (*stress_energy_state) ? (eTtx[index]) : 0.0; CCTK_REAL const eTtyL = (*stress_energy_state) ? (eTty[index]) : 0.0; CCTK_REAL const eTtzL = (*stress_energy_state) ? (eTtz[index]) : 0.0; @@ -228,6 +228,7 @@ void ML_BSSN_O2_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int CCTK_REAL const gt23L = gt23[index]; CCTK_REAL const gt33L = gt33[index]; CCTK_REAL const phiL = phi[index]; + CCTK_REAL const rL = r[index]; CCTK_REAL const trKL = trK[index]; CCTK_REAL const Xt1L = Xt1[index]; CCTK_REAL const Xt2L = Xt2[index]; @@ -492,6 +493,8 @@ void ML_BSSN_O2_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int 2*(PDstandardNth1beta1 + PDstandardNth2beta2 + PDstandardNth3beta3)*Xtn3 - 3*(PDstandardNth1beta3*Xtn1 + PDstandardNth2beta3*Xtn2 + PDstandardNth3beta3*Xtn3)); + CCTK_REAL const eta = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); + CCTK_REAL const beta1rhsL = (PDupwindNth1(beta1, i, j, k)*beta1L + PDupwindNth2(beta1, i, j, k)*beta2L + PDupwindNth3(beta1, i, j, k)*beta3L)*ShiftAdvectionCoeff + B1L*ShiftGammaCoeff; @@ -501,15 +504,15 @@ void ML_BSSN_O2_RHS1_Body(cGH const * restrict const cctkGH, int const dir, int CCTK_REAL const beta3rhsL = (PDupwindNth1(beta3, i, j, k)*beta1L + PDupwindNth2(beta3, i, j, k)*beta2L + PDupwindNth3(beta3, i, j, k)*beta3L)*ShiftAdvectionCoeff + B3L*ShiftGammaCoeff; - CCTK_REAL const B1rhsL = -(B1L*etaL) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + + CCTK_REAL const B1rhsL = -(B1L*eta) + ((PDupwindNth1(B1, i, j, k) - PDupwindNth1(Xt1, i, j, k))*beta1L + (PDupwindNth2(B1, i, j, k) - PDupwindNth2(Xt1, i, j, k))*beta2L + (PDupwindNth3(B1, i, j, k) - PDupwindNth3(Xt1, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt1rhsL; - CCTK_REAL const B2rhsL = -(B2L*etaL) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + + CCTK_REAL const B2rhsL = -(B2L*eta) + ((PDupwindNth1(B2, i, j, k) - PDupwindNth1(Xt2, i, j, k))*beta1L + (PDupwindNth2(B2, i, j, k) - PDupwindNth2(Xt2, i, j, k))*beta2L + (PDupwindNth3(B2, i, j, k) - PDupwindNth3(Xt2, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt2rhsL; - CCTK_REAL const B3rhsL = -(B3L*etaL) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + + CCTK_REAL const B3rhsL = -(B3L*eta) + ((PDupwindNth1(B3, i, j, k) - PDupwindNth1(Xt3, i, j, k))*beta1L + (PDupwindNth2(B3, i, j, k) - PDupwindNth2(Xt3, i, j, k))*beta2L + (PDupwindNth3(B3, i, j, k) - PDupwindNth3(Xt3, i, j, k))*beta3L)*ShiftAdvectionCoeff + Xt3rhsL; diff --git a/ML_BSSN_O2/src/ML_BSSN_O2_setBetaDriverConstant.c b/ML_BSSN_O2/src/ML_BSSN_O2_setBetaDriverConstant.c deleted file mode 100644 index b8da939..0000000 --- a/ML_BSSN_O2/src/ML_BSSN_O2_setBetaDriverConstant.c +++ /dev/null @@ -1,143 +0,0 @@ -/* File produced by Kranc */ - -#define KRANC_C - -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.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_O2_setBetaDriverConstant_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - /* Declare finite differencing variables */ - - /* Declare predefined quantities */ - // CCTK_REAL p1o2dx = INITVALUE; - // CCTK_REAL p1o2dy = INITVALUE; - // CCTK_REAL p1o2dz = INITVALUE; - // CCTK_REAL p1o4dxdy = INITVALUE; - // CCTK_REAL p1o4dxdz = INITVALUE; - // CCTK_REAL p1o4dydz = INITVALUE; - // CCTK_REAL p1odx = INITVALUE; - // CCTK_REAL p1odx2 = INITVALUE; - // CCTK_REAL p1ody = INITVALUE; - // CCTK_REAL p1ody2 = INITVALUE; - // CCTK_REAL p1odz = INITVALUE; - // CCTK_REAL p1odz2 = INITVALUE; - // CCTK_REAL pm1o2dx = INITVALUE; - // CCTK_REAL pm1o2dy = INITVALUE; - // CCTK_REAL pm1o2dz = INITVALUE; - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_O2_setBetaDriverConstant_Body"); - } - - if (cctk_iteration % ML_BSSN_O2_setBetaDriverConstant_calc_every != ML_BSSN_O2_setBetaDriverConstant_calc_offset) - { - return; - } - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - CCTK_REAL const dx = CCTK_DELTA_SPACE(0); - CCTK_REAL const dy = CCTK_DELTA_SPACE(1); - CCTK_REAL const dz = CCTK_DELTA_SPACE(2); - int const di = 1; - int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); - int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); - CCTK_REAL const dxi = 1.0 / dx; - CCTK_REAL const dyi = 1.0 / dy; - CCTK_REAL const dzi = 1.0 / dz; - CCTK_REAL const khalf = 0.5; - CCTK_REAL const kthird = 1/3.0; - CCTK_REAL const ktwothird = 2.0/3.0; - CCTK_REAL const kfourthird = 4.0/3.0; - CCTK_REAL const keightthird = 8.0/3.0; - CCTK_REAL const hdxi = 0.5 * dxi; - CCTK_REAL const hdyi = 0.5 * dyi; - CCTK_REAL const hdzi = 0.5 * dzi; - - /* Initialize predefined quantities */ - CCTK_REAL const p1o2dx = khalf*INV(dx); - CCTK_REAL const p1o2dy = khalf*INV(dy); - CCTK_REAL const p1o2dz = khalf*INV(dz); - CCTK_REAL const p1o4dxdy = (INV(dx)*INV(dy))/4.; - CCTK_REAL const p1o4dxdz = (INV(dx)*INV(dz))/4.; - CCTK_REAL const p1o4dydz = (INV(dy)*INV(dz))/4.; - CCTK_REAL const p1odx = INV(dx); - CCTK_REAL const p1odx2 = pow(dx,-2); - CCTK_REAL const p1ody = INV(dy); - CCTK_REAL const p1ody2 = pow(dy,-2); - CCTK_REAL const p1odz = INV(dz); - CCTK_REAL const p1odz2 = pow(dz,-2); - CCTK_REAL const pm1o2dx = -(khalf*INV(dx)); - CCTK_REAL const pm1o2dy = -(khalf*INV(dy)); - CCTK_REAL const pm1o2dz = -(khalf*INV(dz)); - - /* Loop over the grid points */ - #pragma omp parallel - LC_LOOP3 (ML_BSSN_O2_setBetaDriverConstant, - 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; - int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); - int const 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 etaL = 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 */ - CCTK_REAL const etaL = BetaDriver; - - - /* Copy local copies back to grid functions */ - eta[index] = etaL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (ML_BSSN_O2_setBetaDriverConstant); -} - -void ML_BSSN_O2_setBetaDriverConstant(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - GenericFD_LoopOverEverything(cctkGH, &ML_BSSN_O2_setBetaDriverConstant_Body); -} diff --git a/ML_BSSN_O2/src/ML_BSSN_O2_setBetaDriverSpatial.c b/ML_BSSN_O2/src/ML_BSSN_O2_setBetaDriverSpatial.c deleted file mode 100644 index 902f80e..0000000 --- a/ML_BSSN_O2/src/ML_BSSN_O2_setBetaDriverSpatial.c +++ /dev/null @@ -1,145 +0,0 @@ -/* File produced by Kranc */ - -#define KRANC_C - -#include <assert.h> -#include <math.h> -#include <stdio.h> -#include <stdlib.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_O2_setBetaDriverSpatial_Body(cGH const * restrict const cctkGH, int const dir, int const face, CCTK_REAL const normal[3], CCTK_REAL const tangentA[3], CCTK_REAL const tangentB[3], int const min[3], int const max[3], int const n_subblock_gfs, CCTK_REAL * restrict const subblock_gfs[]) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - - /* Declare finite differencing variables */ - - /* Declare predefined quantities */ - // CCTK_REAL p1o2dx = INITVALUE; - // CCTK_REAL p1o2dy = INITVALUE; - // CCTK_REAL p1o2dz = INITVALUE; - // CCTK_REAL p1o4dxdy = INITVALUE; - // CCTK_REAL p1o4dxdz = INITVALUE; - // CCTK_REAL p1o4dydz = INITVALUE; - // CCTK_REAL p1odx = INITVALUE; - // CCTK_REAL p1odx2 = INITVALUE; - // CCTK_REAL p1ody = INITVALUE; - // CCTK_REAL p1ody2 = INITVALUE; - // CCTK_REAL p1odz = INITVALUE; - // CCTK_REAL p1odz2 = INITVALUE; - // CCTK_REAL pm1o2dx = INITVALUE; - // CCTK_REAL pm1o2dy = INITVALUE; - // CCTK_REAL pm1o2dz = INITVALUE; - - if (verbose > 1) - { - CCTK_VInfo(CCTK_THORNSTRING,"Entering ML_BSSN_O2_setBetaDriverSpatial_Body"); - } - - if (cctk_iteration % ML_BSSN_O2_setBetaDriverSpatial_calc_every != ML_BSSN_O2_setBetaDriverSpatial_calc_offset) - { - return; - } - - /* Include user-supplied include files */ - - /* Initialise finite differencing variables */ - CCTK_REAL const dx = CCTK_DELTA_SPACE(0); - CCTK_REAL const dy = CCTK_DELTA_SPACE(1); - CCTK_REAL const dz = CCTK_DELTA_SPACE(2); - int const di = 1; - int const dj = CCTK_GFINDEX3D(cctkGH,0,1,0) - CCTK_GFINDEX3D(cctkGH,0,0,0); - int const dk = CCTK_GFINDEX3D(cctkGH,0,0,1) - CCTK_GFINDEX3D(cctkGH,0,0,0); - CCTK_REAL const dxi = 1.0 / dx; - CCTK_REAL const dyi = 1.0 / dy; - CCTK_REAL const dzi = 1.0 / dz; - CCTK_REAL const khalf = 0.5; - CCTK_REAL const kthird = 1/3.0; - CCTK_REAL const ktwothird = 2.0/3.0; - CCTK_REAL const kfourthird = 4.0/3.0; - CCTK_REAL const keightthird = 8.0/3.0; - CCTK_REAL const hdxi = 0.5 * dxi; - CCTK_REAL const hdyi = 0.5 * dyi; - CCTK_REAL const hdzi = 0.5 * dzi; - - /* Initialize predefined quantities */ - CCTK_REAL const p1o2dx = khalf*INV(dx); - CCTK_REAL const p1o2dy = khalf*INV(dy); - CCTK_REAL const p1o2dz = khalf*INV(dz); - CCTK_REAL const p1o4dxdy = (INV(dx)*INV(dy))/4.; - CCTK_REAL const p1o4dxdz = (INV(dx)*INV(dz))/4.; - CCTK_REAL const p1o4dydz = (INV(dy)*INV(dz))/4.; - CCTK_REAL const p1odx = INV(dx); - CCTK_REAL const p1odx2 = pow(dx,-2); - CCTK_REAL const p1ody = INV(dy); - CCTK_REAL const p1ody2 = pow(dy,-2); - CCTK_REAL const p1odz = INV(dz); - CCTK_REAL const p1odz2 = pow(dz,-2); - CCTK_REAL const pm1o2dx = -(khalf*INV(dx)); - CCTK_REAL const pm1o2dy = -(khalf*INV(dy)); - CCTK_REAL const pm1o2dz = -(khalf*INV(dz)); - - /* Loop over the grid points */ - #pragma omp parallel - LC_LOOP3 (ML_BSSN_O2_setBetaDriverSpatial, - 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; - int const index = CCTK_GFINDEX3D(cctkGH,i,j,k); - int const 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 etaL = INITVALUE; - // CCTK_REAL rL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - CCTK_REAL const rL = r[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 */ - CCTK_REAL const etaL = BetaDriver*IfThen(rL > SpatialBetaDriverRadius,SpatialBetaDriverRadius*INV(rL),1); - - - /* Copy local copies back to grid functions */ - eta[index] = etaL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (ML_BSSN_O2_setBetaDriverSpatial); -} - -void ML_BSSN_O2_setBetaDriverSpatial(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - GenericFD_LoopOverEverything(cctkGH, &ML_BSSN_O2_setBetaDriverSpatial_Body); -} diff --git a/ML_BSSN_O2/src/RegisterSymmetries.c b/ML_BSSN_O2/src/RegisterSymmetries.c index ce4e6be..0a85ef2 100644 --- a/ML_BSSN_O2/src/RegisterSymmetries.c +++ b/ML_BSSN_O2/src/RegisterSymmetries.c @@ -144,11 +144,6 @@ void ML_BSSN_O2_RegisterSymmetries(CCTK_ARGUMENTS) sym[0] = 1; sym[1] = 1; sym[2] = 1; - SetCartSymVN(cctkGH, sym, "ML_BSSN_O2::eta"); - - sym[0] = 1; - sym[1] = 1; - sym[2] = 1; SetCartSymVN(cctkGH, sym, "ML_BSSN_O2::cS"); sym[0] = -1; diff --git a/ML_BSSN_O2/src/make.code.defn b/ML_BSSN_O2/src/make.code.defn index 70c9a0d..0ef0d8b 100644 --- a/ML_BSSN_O2/src/make.code.defn +++ b/ML_BSSN_O2/src/make.code.defn @@ -1,3 +1,3 @@ # File produced by Kranc -SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_O2_Minkowski.c ML_BSSN_O2_convertFromADMBase.c ML_BSSN_O2_convertFromADMBaseGamma.c ML_BSSN_O2_setBetaDriverConstant.c ML_BSSN_O2_setBetaDriverSpatial.c ML_BSSN_O2_RHS.c ML_BSSN_O2_RHS1.c ML_BSSN_O2_RHS2.c ML_BSSN_O2_RHSStaticBoundary.c ML_BSSN_O2_RHSRadiativeBoundary.c ML_BSSN_O2_enforce.c ML_BSSN_O2_enforce2.c ML_BSSN_O2_boundary.c ML_BSSN_O2_convertToADMBase.c ML_BSSN_O2_convertToADMBaseDtLapseShift.c ML_BSSN_O2_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_O2_convertToADMBaseFakeDtLapseShift.c ML_BSSN_O2_constraints.c ML_BSSN_O2_constraints_boundary.c Boundaries.c +SRCS = Startup.c RegisterMoL.c RegisterSymmetries.c ML_BSSN_O2_Minkowski.c ML_BSSN_O2_convertFromADMBase.c ML_BSSN_O2_convertFromADMBaseGamma.c ML_BSSN_O2_RHS.c ML_BSSN_O2_RHS1.c ML_BSSN_O2_RHS2.c ML_BSSN_O2_RHSStaticBoundary.c ML_BSSN_O2_RHSRadiativeBoundary.c ML_BSSN_O2_enforce.c ML_BSSN_O2_enforce2.c ML_BSSN_O2_boundary.c ML_BSSN_O2_convertToADMBase.c ML_BSSN_O2_convertToADMBaseDtLapseShift.c ML_BSSN_O2_convertToADMBaseDtLapseShiftBoundary.c ML_BSSN_O2_convertToADMBaseFakeDtLapseShift.c ML_BSSN_O2_constraints.c ML_BSSN_O2_constraints_boundary.c Boundaries.c diff --git a/ML_BSSN_O2_Helper/src/SetGroupTags.c b/ML_BSSN_O2_Helper/src/SetGroupTags.c index cad3293..54b9df3 100644 --- a/ML_BSSN_O2_Helper/src/SetGroupTags.c +++ b/ML_BSSN_O2_Helper/src/SetGroupTags.c @@ -28,8 +28,6 @@ ML_BSSN_O2_SetGroupTags (void) set_group_tags (0, 0, 0, "ML_BSSN_O2::ML_Ham"); set_group_tags (0, 0, 0, "ML_BSSN_O2::ML_mom"); - set_group_tags (0, 1, 0, "ML_BSSN_O2::ML_BetaDriver"); - int const checkpoint = rhs_timelevels > 1; set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN_O2::ML_dtlapserhs"); set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN_O2::ML_dtshiftrhs"); diff --git a/m/McLachlan_BSSN.m b/m/McLachlan_BSSN.m index bb1f094..99fa876 100644 --- a/m/McLachlan_BSSN.m +++ b/m/McLachlan_BSSN.m @@ -140,7 +140,7 @@ Map [DefineTensor, e4phi, em4phi, ddetg, detgt, gtu, ddetgt, dgtu, ddgtu, Gtl, Gtlu, Gt, Rt, Rphi, gK, T00, T0, T, rho, S, - eta, x, y, z, r, + x, y, z, r, Psi0re, Psi0im, Psi1re, Psi1im, Psi2re, Psi2im, Psi3re, Psi3im, Psi4re, Psi4im, er, eth, eph, mm1A, mm1L, mm1, mm2A, mm2B, mm2L, mm2, @@ -220,8 +220,7 @@ evolvedGroups = SetGroupName [CreateGroupFromTensor [alpha ], prefix <> "lapse" ], SetGroupName [CreateGroupFromTensor [A ], prefix <> "dtlapse" ], SetGroupName [CreateGroupFromTensor [beta[ua] ], prefix <> "shift" ], - SetGroupName [CreateGroupFromTensor [B[ua] ], prefix <> "dtshift" ], - SetGroupName [CreateGroupFromTensor [eta ], prefix <> "BetaDriver"]}; + SetGroupName [CreateGroupFromTensor [B[ua] ], prefix <> "dtshift" ]}; evaluatedGroups = {SetGroupName [CreateGroupFromTensor [H ], prefix <> "Ham"], SetGroupName [CreateGroupFromTensor [M[la] ], prefix <> "mom"], @@ -361,29 +360,6 @@ convertFromADMBaseGammaCalc = } }; -setBetaDriverConstantCalc = -{ - Name -> BSSN <> "_setBetaDriverConstant", - Schedule -> {"IN "<> BSSN <> "_InitEta"}, - ConditionalOnKeyword -> {"UseSpatialBetaDriver", "no"}, - Equations -> - { - eta -> BetaDriver - } -}; - -setBetaDriverSpatialCalc = -{ - Name -> BSSN <> "_setBetaDriverSpatial", - Schedule -> {"IN "<> BSSN <> "_InitEta"}, - ConditionalOnKeyword -> {"UseSpatialBetaDriver", "yes"}, - Equations -> - { - eta -> BetaDriver - IfThen [r > SpatialBetaDriverRadius, SpatialBetaDriverRadius / r, 1] - } -}; - (******************************************************************************) (* Convert to ADMBase *) (******************************************************************************) @@ -477,7 +453,7 @@ evolCalc = Gt[ua,lb,lc], Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb], Atm[ua,lb], Atu[ua,ub], e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg, - gu[ua,ub], G[ua,lb,lc], Ats[la,lb], trAts, + gu[ua,ub], G[ua,lb,lc], Ats[la,lb], trAts, eta, rho, S[la], trS, fac1, fac2}, Equations -> { @@ -596,7 +572,10 @@ evolCalc = + LapseAdvectionCoeff beta[ua] PDu[alpha,la], dot[A] -> (1 - LapseAdvectionCoeff) (dot[trK] - AlphaDriver A), - + + eta -> BetaDriver + IfThen [r > SpatialBetaDriverRadius, SpatialBetaDriverRadius / r, 1], + (* dot[beta[ua]] -> eta Xt[ua], *) (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *) dot[beta[ua]] -> + ShiftGammaCoeff B[ua] @@ -623,7 +602,7 @@ evol1Calc = detgt, gtu[ua,ub], Gt[ua,lb,lc], Xtn[ua], Rt[la,lb], Rphi[la,lb], R[la,lb], Atm[ua,lb], Atu[ua,ub], - e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg, + e4phi, em4phi, cdphi[la], cdphi2[la,lb], g[la,lb], detg, eta, rho, S[la], trS, fac1, fac2}, Equations -> { @@ -678,6 +657,9 @@ evol1Calc = (* Equation (4.28) in Baumgarte & Shapiro (Phys. Rept. 376 (2003) 41-131) *) + addMatter (- 16 pi alpha gtu[ui,uj] S[lj]), + eta -> BetaDriver + IfThen [r > SpatialBetaDriverRadius, SpatialBetaDriverRadius / r, 1], + (* dot[beta[ua]] -> eta Xt[ua], *) (* dot[beta[ua]] -> ShiftGammaCoeff alpha^ShiftAlphaPower B[ua], *) dot[beta[ua]] -> + ShiftGammaCoeff B[ua] @@ -1217,8 +1199,6 @@ calculations = initialCalc, convertFromADMBaseCalc, convertFromADMBaseGammaCalc, - setBetaDriverConstantCalc, - setBetaDriverSpatialCalc, evolCalc, evol1Calc, evol2Calc, RHSStaticBoundaryCalc, diff --git a/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c b/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c index 3113637..0666c1b 100644 --- a/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c +++ b/m/prototype/ML_BSSN_Helper/src/SetGroupTags.c @@ -28,8 +28,6 @@ ML_BSSN_SetGroupTags (void) set_group_tags (0, 0, 0, "ML_BSSN::ML_Ham"); set_group_tags (0, 0, 0, "ML_BSSN::ML_mom"); - set_group_tags (0, 1, 0, "ML_BSSN::ML_BetaDriver"); - int const checkpoint = rhs_timelevels > 1; set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_dtlapserhs"); set_group_tags (checkpoint, checkpoint, 0, "ML_BSSN::ML_dtshiftrhs"); |