diff options
| author | Erik Schnetter <eschnett@redshift.lan> | 2007-11-04 19:29:11 -0600 |
|---|---|---|
| committer | Erik Schnetter <eschnett@redshift.lan> | 2007-11-04 19:29:11 -0600 |
| commit | bb691c9b17b5cce63e6ecbe084f9d32902da2a5d (patch) | |
| tree | 90d3d34a472fc4c18e7671615666f2f3319808c8 | |
| parent | 23f33ea59a04e5cb6c663821c6a6ecbb8d2f438c (diff) | |
Remove old files
| -rw-r--r-- | ML_ADM/src/ADM_Minkowski.c | 180 | ||||
| -rw-r--r-- | ML_ADM/src/ADM_RHS.c | 536 | ||||
| -rw-r--r-- | ML_ADM/src/ADM_constraints.c | 445 | ||||
| -rw-r--r-- | ML_ADM/src/ADM_convertFromADMBase.c | 212 | ||||
| -rw-r--r-- | ML_ADM/src/ADM_convertToADMBase.c | 228 | ||||
| -rw-r--r-- | ML_BSSN/src/BSSN_Minkowski.c | 212 | ||||
| -rw-r--r-- | ML_BSSN/src/BSSN_constraints.c | 980 | ||||
| -rw-r--r-- | ML_BSSN/src/BSSN_convertFromADMBase.c | 295 | ||||
| -rw-r--r-- | ML_BSSN/src/BSSN_convertFromADMBaseGamma.c | 254 | ||||
| -rw-r--r-- | ML_BSSN/src/BSSN_convertToADMBase.c | 267 | ||||
| -rw-r--r-- | ML_BSSN/src/BSSN_enforce.c | 179 | ||||
| -rw-r--r-- | McLachlan.m | 544 | ||||
| -rwxr-xr-x | McLachlan.sh | 6 | ||||
| -rw-r--r-- | m/McLachlan.out | 755 |
14 files changed, 0 insertions, 5093 deletions
diff --git a/ML_ADM/src/ADM_Minkowski.c b/ML_ADM/src/ADM_Minkowski.c deleted file mode 100644 index c21c98e..0000000 --- a/ML_ADM/src/ADM_Minkowski.c +++ /dev/null @@ -1,180 +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 ADM_Minkowski_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 ADM_Minkowski_Body"); - } - - if (cctk_iteration % ADM_Minkowski_calc_every != ADM_Minkowski_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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, g21L = INITVALUE, g22L = INITVALUE, g31L = INITVALUE, g32L = INITVALUE, g33L = INITVALUE; - CCTK_REAL K11L = INITVALUE, K21L = INITVALUE, K22L = INITVALUE, K31L = INITVALUE, K32L = 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; - - g21L = 0; - - g31L = 0; - - g22L = 1; - - g32L = 0; - - g33L = 1; - - K11L = 0; - - K21L = 0; - - K31L = 0; - - K22L = 0; - - K32L = 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; - g21[index] = g21L; - g22[index] = g22L; - g31[index] = g31L; - g32[index] = g32L; - g33[index] = g33L; - K11[index] = K11L; - K21[index] = K21L; - K22[index] = K22L; - K31[index] = K31L; - K32[index] = K32L; - K33[index] = K33L; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void ADM_Minkowski(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &ADM_Minkowski_Body); -} diff --git a/ML_ADM/src/ADM_RHS.c b/ML_ADM/src/ADM_RHS.c deleted file mode 100644 index 1a1a7ed..0000000 --- a/ML_ADM/src/ADM_RHS.c +++ /dev/null @@ -1,536 +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 ADM_RHS_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 ADM_RHS_Body"); - } - - if (cctk_iteration % ADM_RHS_calc_every != ADM_RHS_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 detg = INITVALUE; - CCTK_REAL G111 = INITVALUE, G121 = INITVALUE, G122 = INITVALUE, G131 = INITVALUE, G132 = INITVALUE, G133 = INITVALUE; - CCTK_REAL G211 = INITVALUE, G221 = INITVALUE, G222 = INITVALUE, G231 = INITVALUE, G232 = INITVALUE, G233 = INITVALUE; - CCTK_REAL G311 = INITVALUE, G321 = INITVALUE, G322 = INITVALUE, G331 = INITVALUE, G332 = INITVALUE, G333 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE, gu33 = INITVALUE; - CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE; - CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE; - CCTK_REAL R11 = INITVALUE, R21 = INITVALUE, R22 = INITVALUE, R31 = INITVALUE, R32 = INITVALUE, R33 = INITVALUE; - CCTK_REAL trK = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL alphaL = INITVALUE, alpharhsL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta1rhsL = INITVALUE, beta2L = INITVALUE, beta2rhsL = INITVALUE, beta3L = INITVALUE, beta3rhsL = INITVALUE; - CCTK_REAL g11L = INITVALUE, g11rhsL = INITVALUE, g21L = INITVALUE, g21rhsL = INITVALUE, g22L = INITVALUE, g22rhsL = INITVALUE; - CCTK_REAL g31L = INITVALUE, g31rhsL = INITVALUE, g32L = INITVALUE, g32rhsL = INITVALUE, g33L = INITVALUE, g33rhsL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K11rhsL = INITVALUE, K21L = INITVALUE, K21rhsL = INITVALUE, K22L = INITVALUE, K22rhsL = INITVALUE; - CCTK_REAL K31L = INITVALUE, K31rhsL = INITVALUE, K32L = INITVALUE, K32rhsL = INITVALUE, K33L = INITVALUE, K33rhsL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - CCTK_REAL PDstandard4th1alpha = INITVALUE; - CCTK_REAL PDstandard4th2alpha = INITVALUE; - CCTK_REAL PDstandard4th3alpha = INITVALUE; - CCTK_REAL PDstandard4th11alpha = INITVALUE; - CCTK_REAL PDstandard4th22alpha = INITVALUE; - CCTK_REAL PDstandard4th33alpha = INITVALUE; - CCTK_REAL PDstandard4th12alpha = INITVALUE; - CCTK_REAL PDstandard4th13alpha = INITVALUE; - CCTK_REAL PDstandard4th23alpha = INITVALUE; - CCTK_REAL PDstandard4th1beta1 = INITVALUE; - CCTK_REAL PDstandard4th2beta1 = INITVALUE; - CCTK_REAL PDstandard4th3beta1 = INITVALUE; - CCTK_REAL PDstandard4th1beta2 = INITVALUE; - CCTK_REAL PDstandard4th2beta2 = INITVALUE; - CCTK_REAL PDstandard4th3beta2 = INITVALUE; - CCTK_REAL PDstandard4th1beta3 = INITVALUE; - CCTK_REAL PDstandard4th2beta3 = INITVALUE; - CCTK_REAL PDstandard4th3beta3 = INITVALUE; - CCTK_REAL PDstandard4th1g11 = INITVALUE; - CCTK_REAL PDstandard4th2g11 = INITVALUE; - CCTK_REAL PDstandard4th3g11 = INITVALUE; - CCTK_REAL PDstandard4th22g11 = INITVALUE; - CCTK_REAL PDstandard4th33g11 = INITVALUE; - CCTK_REAL PDstandard4th12g11 = INITVALUE; - CCTK_REAL PDstandard4th13g11 = INITVALUE; - CCTK_REAL PDstandard4th21g11 = INITVALUE; - CCTK_REAL PDstandard4th23g11 = INITVALUE; - CCTK_REAL PDstandard4th31g11 = INITVALUE; - CCTK_REAL PDstandard4th32g11 = INITVALUE; - CCTK_REAL PDstandard4th1g21 = INITVALUE; - CCTK_REAL PDstandard4th2g21 = INITVALUE; - CCTK_REAL PDstandard4th3g21 = INITVALUE; - CCTK_REAL PDstandard4th33g21 = INITVALUE; - CCTK_REAL PDstandard4th12g21 = INITVALUE; - CCTK_REAL PDstandard4th13g21 = INITVALUE; - CCTK_REAL PDstandard4th21g21 = INITVALUE; - CCTK_REAL PDstandard4th23g21 = INITVALUE; - CCTK_REAL PDstandard4th31g21 = INITVALUE; - CCTK_REAL PDstandard4th32g21 = INITVALUE; - CCTK_REAL PDstandard4th1g22 = INITVALUE; - CCTK_REAL PDstandard4th2g22 = INITVALUE; - CCTK_REAL PDstandard4th3g22 = INITVALUE; - CCTK_REAL PDstandard4th11g22 = INITVALUE; - CCTK_REAL PDstandard4th33g22 = INITVALUE; - CCTK_REAL PDstandard4th12g22 = INITVALUE; - CCTK_REAL PDstandard4th13g22 = INITVALUE; - CCTK_REAL PDstandard4th21g22 = INITVALUE; - CCTK_REAL PDstandard4th23g22 = INITVALUE; - CCTK_REAL PDstandard4th31g22 = INITVALUE; - CCTK_REAL PDstandard4th32g22 = INITVALUE; - CCTK_REAL PDstandard4th1g31 = INITVALUE; - CCTK_REAL PDstandard4th2g31 = INITVALUE; - CCTK_REAL PDstandard4th3g31 = INITVALUE; - CCTK_REAL PDstandard4th22g31 = INITVALUE; - CCTK_REAL PDstandard4th12g31 = INITVALUE; - CCTK_REAL PDstandard4th13g31 = INITVALUE; - CCTK_REAL PDstandard4th21g31 = INITVALUE; - CCTK_REAL PDstandard4th23g31 = INITVALUE; - CCTK_REAL PDstandard4th31g31 = INITVALUE; - CCTK_REAL PDstandard4th32g31 = INITVALUE; - CCTK_REAL PDstandard4th1g32 = INITVALUE; - CCTK_REAL PDstandard4th2g32 = INITVALUE; - CCTK_REAL PDstandard4th3g32 = INITVALUE; - CCTK_REAL PDstandard4th11g32 = INITVALUE; - CCTK_REAL PDstandard4th12g32 = INITVALUE; - CCTK_REAL PDstandard4th13g32 = INITVALUE; - CCTK_REAL PDstandard4th21g32 = INITVALUE; - CCTK_REAL PDstandard4th23g32 = INITVALUE; - CCTK_REAL PDstandard4th31g32 = INITVALUE; - CCTK_REAL PDstandard4th32g32 = INITVALUE; - CCTK_REAL PDstandard4th1g33 = INITVALUE; - CCTK_REAL PDstandard4th2g33 = INITVALUE; - CCTK_REAL PDstandard4th3g33 = INITVALUE; - CCTK_REAL PDstandard4th11g33 = INITVALUE; - CCTK_REAL PDstandard4th22g33 = INITVALUE; - CCTK_REAL PDstandard4th12g33 = INITVALUE; - CCTK_REAL PDstandard4th13g33 = INITVALUE; - CCTK_REAL PDstandard4th21g33 = INITVALUE; - CCTK_REAL PDstandard4th23g33 = INITVALUE; - CCTK_REAL PDstandard4th31g33 = INITVALUE; - CCTK_REAL PDstandard4th32g33 = INITVALUE; - CCTK_REAL PDstandard4th1K11 = INITVALUE; - CCTK_REAL PDstandard4th2K11 = INITVALUE; - CCTK_REAL PDstandard4th3K11 = INITVALUE; - CCTK_REAL PDstandard4th1K21 = INITVALUE; - CCTK_REAL PDstandard4th2K21 = INITVALUE; - CCTK_REAL PDstandard4th3K21 = INITVALUE; - CCTK_REAL PDstandard4th1K22 = INITVALUE; - CCTK_REAL PDstandard4th2K22 = INITVALUE; - CCTK_REAL PDstandard4th3K22 = INITVALUE; - CCTK_REAL PDstandard4th1K31 = INITVALUE; - CCTK_REAL PDstandard4th2K31 = INITVALUE; - CCTK_REAL PDstandard4th3K31 = INITVALUE; - CCTK_REAL PDstandard4th1K32 = INITVALUE; - CCTK_REAL PDstandard4th2K32 = INITVALUE; - CCTK_REAL PDstandard4th3K32 = INITVALUE; - CCTK_REAL PDstandard4th1K33 = INITVALUE; - CCTK_REAL PDstandard4th2K33 = INITVALUE; - CCTK_REAL PDstandard4th3K33 = INITVALUE; - - /* Assign local copies of grid functions */ - alphaL = alpha[index]; - beta1L = beta1[index]; - beta2L = beta2[index]; - beta3L = beta3[index]; - g11L = g11[index]; - g21L = g21[index]; - g22L = g22[index]; - g31L = g31[index]; - g32L = g32[index]; - g33L = g33[index]; - K11L = K11[index]; - K21L = K21[index]; - K22L = K22[index]; - K31L = K31[index]; - K32L = K32[index]; - K33L = K33[index]; - - /* Assign local copies of subblock grid functions */ - - /* Include user supplied include files */ - - /* Precompute derivatives (new style) */ - PDstandard4th1alpha = PDstandard4th1(alpha, i, j, k); - PDstandard4th2alpha = PDstandard4th2(alpha, i, j, k); - PDstandard4th3alpha = PDstandard4th3(alpha, i, j, k); - PDstandard4th11alpha = PDstandard4th11(alpha, i, j, k); - PDstandard4th22alpha = PDstandard4th22(alpha, i, j, k); - PDstandard4th33alpha = PDstandard4th33(alpha, i, j, k); - PDstandard4th12alpha = PDstandard4th12(alpha, i, j, k); - PDstandard4th13alpha = PDstandard4th13(alpha, i, j, k); - PDstandard4th23alpha = PDstandard4th23(alpha, i, j, k); - PDstandard4th1beta1 = PDstandard4th1(beta1, i, j, k); - PDstandard4th2beta1 = PDstandard4th2(beta1, i, j, k); - PDstandard4th3beta1 = PDstandard4th3(beta1, i, j, k); - PDstandard4th1beta2 = PDstandard4th1(beta2, i, j, k); - PDstandard4th2beta2 = PDstandard4th2(beta2, i, j, k); - PDstandard4th3beta2 = PDstandard4th3(beta2, i, j, k); - PDstandard4th1beta3 = PDstandard4th1(beta3, i, j, k); - PDstandard4th2beta3 = PDstandard4th2(beta3, i, j, k); - PDstandard4th3beta3 = PDstandard4th3(beta3, i, j, k); - PDstandard4th1g11 = PDstandard4th1(g11, i, j, k); - PDstandard4th2g11 = PDstandard4th2(g11, i, j, k); - PDstandard4th3g11 = PDstandard4th3(g11, i, j, k); - PDstandard4th22g11 = PDstandard4th22(g11, i, j, k); - PDstandard4th33g11 = PDstandard4th33(g11, i, j, k); - PDstandard4th23g11 = PDstandard4th23(g11, i, j, k); - PDstandard4th1g21 = PDstandard4th1(g21, i, j, k); - PDstandard4th2g21 = PDstandard4th2(g21, i, j, k); - PDstandard4th3g21 = PDstandard4th3(g21, i, j, k); - PDstandard4th33g21 = PDstandard4th33(g21, i, j, k); - PDstandard4th12g21 = PDstandard4th12(g21, i, j, k); - PDstandard4th13g21 = PDstandard4th13(g21, i, j, k); - PDstandard4th23g21 = PDstandard4th23(g21, i, j, k); - PDstandard4th1g22 = PDstandard4th1(g22, i, j, k); - PDstandard4th2g22 = PDstandard4th2(g22, i, j, k); - PDstandard4th3g22 = PDstandard4th3(g22, i, j, k); - PDstandard4th11g22 = PDstandard4th11(g22, i, j, k); - PDstandard4th33g22 = PDstandard4th33(g22, i, j, k); - PDstandard4th13g22 = PDstandard4th13(g22, i, j, k); - PDstandard4th1g31 = PDstandard4th1(g31, i, j, k); - PDstandard4th2g31 = PDstandard4th2(g31, i, j, k); - PDstandard4th3g31 = PDstandard4th3(g31, i, j, k); - PDstandard4th22g31 = PDstandard4th22(g31, i, j, k); - PDstandard4th12g31 = PDstandard4th12(g31, i, j, k); - PDstandard4th13g31 = PDstandard4th13(g31, i, j, k); - PDstandard4th23g31 = PDstandard4th23(g31, i, j, k); - PDstandard4th1g32 = PDstandard4th1(g32, i, j, k); - PDstandard4th2g32 = PDstandard4th2(g32, i, j, k); - PDstandard4th3g32 = PDstandard4th3(g32, i, j, k); - PDstandard4th11g32 = PDstandard4th11(g32, i, j, k); - PDstandard4th12g32 = PDstandard4th12(g32, i, j, k); - PDstandard4th13g32 = PDstandard4th13(g32, i, j, k); - PDstandard4th23g32 = PDstandard4th23(g32, i, j, k); - PDstandard4th1g33 = PDstandard4th1(g33, i, j, k); - PDstandard4th2g33 = PDstandard4th2(g33, i, j, k); - PDstandard4th3g33 = PDstandard4th3(g33, i, j, k); - PDstandard4th11g33 = PDstandard4th11(g33, i, j, k); - PDstandard4th22g33 = PDstandard4th22(g33, i, j, k); - PDstandard4th12g33 = PDstandard4th12(g33, i, j, k); - PDstandard4th1K11 = PDstandard4th1(K11, i, j, k); - PDstandard4th2K11 = PDstandard4th2(K11, i, j, k); - PDstandard4th3K11 = PDstandard4th3(K11, i, j, k); - PDstandard4th1K21 = PDstandard4th1(K21, i, j, k); - PDstandard4th2K21 = PDstandard4th2(K21, i, j, k); - PDstandard4th3K21 = PDstandard4th3(K21, i, j, k); - PDstandard4th1K22 = PDstandard4th1(K22, i, j, k); - PDstandard4th2K22 = PDstandard4th2(K22, i, j, k); - PDstandard4th3K22 = PDstandard4th3(K22, i, j, k); - PDstandard4th1K31 = PDstandard4th1(K31, i, j, k); - PDstandard4th2K31 = PDstandard4th2(K31, i, j, k); - PDstandard4th3K31 = PDstandard4th3(K31, i, j, k); - PDstandard4th1K32 = PDstandard4th1(K32, i, j, k); - PDstandard4th2K32 = PDstandard4th2(K32, i, j, k); - PDstandard4th3K32 = PDstandard4th3(K32, i, j, k); - PDstandard4th1K33 = PDstandard4th1(K33, i, j, k); - PDstandard4th2K33 = PDstandard4th2(K33, i, j, k); - PDstandard4th3K33 = PDstandard4th3(K33, i, j, k); - - /* Precompute derivatives (old style) */ - - /* Calculate temporaries and grid functions */ - detg = 2*g21L*g31L*g32L + g33L*(g11L*g22L - SQR(g21L)) - g22L*SQR(g31L) - g11L*SQR(g32L); - - gu11 = INV(detg)*(g22L*g33L - SQR(g32L)); - - gu12 = (g31L*g32L - g21L*g33L)*INV(detg); - - gu13 = (-(g22L*g31L) + g21L*g32L)*INV(detg); - - gu22 = INV(detg)*(g11L*g33L - SQR(g31L)); - - gu23 = (g21L*g31L - g11L*g32L)*INV(detg); - - gu33 = INV(detg)*(g11L*g22L - SQR(g21L)); - - G111 = khalf*(gu11*PDstandard4th1g11 + 2*(gu12*PDstandard4th1g21 + gu13*PDstandard4th1g31) - gu12*PDstandard4th2g11 - - gu13*PDstandard4th3g11); - - G211 = khalf*(gu12*PDstandard4th1g11 + 2*(gu22*PDstandard4th1g21 + gu23*PDstandard4th1g31) - gu22*PDstandard4th2g11 - - gu23*PDstandard4th3g11); - - G311 = khalf*(gu13*PDstandard4th1g11 + 2*(gu23*PDstandard4th1g21 + gu33*PDstandard4th1g31) - gu23*PDstandard4th2g11 - - gu33*PDstandard4th3g11); - - G121 = khalf*(gu12*PDstandard4th1g22 + gu11*PDstandard4th2g11 + - gu13*(PDstandard4th1g32 + PDstandard4th2g31 - PDstandard4th3g21)); - - G221 = khalf*(gu22*PDstandard4th1g22 + gu12*PDstandard4th2g11 + - gu23*(PDstandard4th1g32 + PDstandard4th2g31 - PDstandard4th3g21)); - - G321 = khalf*(gu23*PDstandard4th1g22 + gu13*PDstandard4th2g11 + - gu33*(PDstandard4th1g32 + PDstandard4th2g31 - PDstandard4th3g21)); - - G131 = khalf*(gu13*PDstandard4th1g33 + gu11*PDstandard4th3g11 + - gu12*(PDstandard4th1g32 - PDstandard4th2g31 + PDstandard4th3g21)); - - G231 = khalf*(gu23*PDstandard4th1g33 + gu12*PDstandard4th3g11 + - gu22*(PDstandard4th1g32 - PDstandard4th2g31 + PDstandard4th3g21)); - - G331 = khalf*(gu33*PDstandard4th1g33 + gu13*PDstandard4th3g11 + - gu23*(PDstandard4th1g32 - PDstandard4th2g31 + PDstandard4th3g21)); - - G122 = khalf*(gu11*(-PDstandard4th1g22 + 2*PDstandard4th2g21) + gu12*PDstandard4th2g22 + - gu13*(2*PDstandard4th2g32 - PDstandard4th3g22)); - - G222 = khalf*(gu12*(-PDstandard4th1g22 + 2*PDstandard4th2g21) + gu22*PDstandard4th2g22 + - gu23*(2*PDstandard4th2g32 - PDstandard4th3g22)); - - G322 = khalf*(gu13*(-PDstandard4th1g22 + 2*PDstandard4th2g21) + gu23*PDstandard4th2g22 + - gu33*(2*PDstandard4th2g32 - PDstandard4th3g22)); - - G132 = khalf*(gu13*PDstandard4th2g33 + gu11*(-PDstandard4th1g32 + PDstandard4th2g31 + PDstandard4th3g21) + - gu12*PDstandard4th3g22); - - G232 = khalf*(gu23*PDstandard4th2g33 + gu12*(-PDstandard4th1g32 + PDstandard4th2g31 + PDstandard4th3g21) + - gu22*PDstandard4th3g22); - - G332 = khalf*(gu33*PDstandard4th2g33 + gu13*(-PDstandard4th1g32 + PDstandard4th2g31 + PDstandard4th3g21) + - gu23*PDstandard4th3g22); - - G133 = khalf*(-(gu11*PDstandard4th1g33) - gu12*PDstandard4th2g33 + 2*gu11*PDstandard4th3g31 + - 2*gu12*PDstandard4th3g32 + gu13*PDstandard4th3g33); - - G233 = khalf*(-(gu12*PDstandard4th1g33) - gu22*PDstandard4th2g33 + 2*gu12*PDstandard4th3g31 + - 2*gu22*PDstandard4th3g32 + gu23*PDstandard4th3g33); - - G333 = khalf*(-(gu13*PDstandard4th1g33) - gu23*PDstandard4th2g33 + 2*gu13*PDstandard4th3g31 + - 2*gu23*PDstandard4th3g32 + gu33*PDstandard4th3g33); - - R11 = khalf*(4*G231*G321 + G211*(2*G121 - 2*G222 - 2*G332) + G311*(2*G131 - 2*G333) - gu22*PDstandard4th11g22 - - 2*(G111*G221 + G232*G311 + G111*G331 + gu23*PDstandard4th11g32) - gu33*PDstandard4th11g33 + - 2*gu22*PDstandard4th12g21 + 2*gu23*PDstandard4th12g31 + 2*gu23*PDstandard4th13g21 + 2*gu33*PDstandard4th13g31 - - gu22*PDstandard4th22g11 - 2*gu23*PDstandard4th23g11 - gu33*PDstandard4th33g11 + 2*SQR(G221) + 2*SQR(G331)); - - R21 = khalf*(2*(G122*G211 + G132*G311 + G231*G322 + G331*G332) - - 2*(G121*G221 + G121*G331 + G221*G332 + G321*G333 + gu12*PDstandard4th12g21) - gu23*PDstandard4th12g32 - - gu33*PDstandard4th12g33 + gu13*(PDstandard4th11g32 - PDstandard4th12g31 - PDstandard4th13g21) + - gu23*PDstandard4th13g22 + gu33*PDstandard4th13g32 + gu12*(PDstandard4th11g22 + PDstandard4th22g11) + - gu23*PDstandard4th22g31 + gu13*PDstandard4th23g11 - gu23*PDstandard4th23g21 + gu33*PDstandard4th23g31 - - gu33*PDstandard4th33g21); - - R31 = khalf*(2*(G132*G211 + G221*G232 + G133*G311 + G233*G321) - - 2*(G222*G231 + G232*G331 + G131*(G221 + G331) + gu13*PDstandard4th13g31) + - gu12*(PDstandard4th11g32 - PDstandard4th12g31 - PDstandard4th13g21 + PDstandard4th23g11) + - gu22*(PDstandard4th12g32 - PDstandard4th13g22 - PDstandard4th22g31 + PDstandard4th23g21) + - gu13*(PDstandard4th11g33 + PDstandard4th33g11) + - gu23*(PDstandard4th12g33 - PDstandard4th13g32 - PDstandard4th23g31 + PDstandard4th33g21)); - - R22 = khalf*(4*G132*G321 + G122*(2*G221 - 2*(G111 + G331)) + G322*(2*G232 - 2*G333) - - 2*(G131*G322 + G222*(G121 + G332) + gu13*PDstandard4th13g22) + - gu11*(-PDstandard4th11g22 + 2*PDstandard4th12g21 - PDstandard4th22g11) + - gu13*(-2*PDstandard4th22g31 + 2*(PDstandard4th12g32 + PDstandard4th23g21)) + - gu33*(-PDstandard4th22g33 + 2*PDstandard4th23g32 - PDstandard4th33g22) + 2*(SQR(G121) + SQR(G332))); - - R32 = khalf*(2*(G121*G131 + G122*G231 + G133*G321 + G233*G322) + - gu11*(-PDstandard4th11g32 + PDstandard4th12g31 + PDstandard4th13g21 - PDstandard4th23g11) + - gu12*(-PDstandard4th12g32 + PDstandard4th13g22 + PDstandard4th22g31 - PDstandard4th23g21) - - 2*(G111*G132 + G131*G332 + G232*(G121 + G332) + gu23*PDstandard4th23g32) + - gu13*(PDstandard4th12g33 - PDstandard4th13g32 - PDstandard4th23g31 + PDstandard4th33g21) + - gu23*(PDstandard4th22g33 + PDstandard4th33g22)); - - R33 = khalf*(4*G132*G231 - gu11*PDstandard4th11g33 - - 2*(G111*G133 + G133*G221 + G121*G233 + G222*G233 + G131*G333 + G232*G333 + gu12*PDstandard4th12g33) + - 2*(G133*G331 + G233*G332 + gu11*PDstandard4th13g31) + 2*gu12*PDstandard4th13g32 - gu22*PDstandard4th22g33 + - 2*gu12*PDstandard4th23g31 + 2*gu22*PDstandard4th23g32 - gu11*PDstandard4th33g11 - 2*gu12*PDstandard4th33g21 - - gu22*PDstandard4th33g22 + 2*SQR(G131) + 2*SQR(G232)); - - Km11 = gu11*K11L + gu12*K21L + gu13*K31L; - - Km21 = gu12*K11L + gu22*K21L + gu23*K31L; - - Km31 = gu13*K11L + gu23*K21L + gu33*K31L; - - Km12 = gu11*K21L + gu12*K22L + gu13*K32L; - - Km22 = gu12*K21L + gu22*K22L + gu23*K32L; - - Km32 = gu13*K21L + gu23*K22L + gu33*K32L; - - Km13 = gu11*K31L + gu12*K32L + gu13*K33L; - - Km23 = gu12*K31L + gu22*K32L + gu23*K33L; - - Km33 = gu13*K31L + gu23*K32L + gu33*K33L; - - trK = Km11 + Km22 + Km33; - - g11rhsL = -2*alphaL*K11L + 2*(g11L*PDstandard4th1beta1 + g21L*PDstandard4th1beta2 + g31L*PDstandard4th1beta3) + - beta1L*PDstandard4th1g11 + beta2L*PDstandard4th2g11 + beta3L*PDstandard4th3g11; - - g21rhsL = -2*alphaL*K21L + g22L*PDstandard4th1beta2 + g32L*PDstandard4th1beta3 + beta1L*PDstandard4th1g21 + - g11L*PDstandard4th2beta1 + g21L*(PDstandard4th1beta1 + PDstandard4th2beta2) + g31L*PDstandard4th2beta3 + - beta2L*PDstandard4th2g21 + beta3L*PDstandard4th3g21; - - g31rhsL = -2*alphaL*K31L + g32L*PDstandard4th1beta2 + g33L*PDstandard4th1beta3 + beta1L*PDstandard4th1g31 + - beta2L*PDstandard4th2g31 + g11L*PDstandard4th3beta1 + g21L*PDstandard4th3beta2 + - g31L*(PDstandard4th1beta1 + PDstandard4th3beta3) + beta3L*PDstandard4th3g31; - - g22rhsL = -2*alphaL*K22L + beta1L*PDstandard4th1g22 + - 2*(g21L*PDstandard4th2beta1 + g22L*PDstandard4th2beta2 + g32L*PDstandard4th2beta3) + beta2L*PDstandard4th2g22 + - beta3L*PDstandard4th3g22; - - g32rhsL = -2*alphaL*K32L + beta1L*PDstandard4th1g32 + g31L*PDstandard4th2beta1 + g33L*PDstandard4th2beta3 + - beta2L*PDstandard4th2g32 + g21L*PDstandard4th3beta1 + g22L*PDstandard4th3beta2 + - g32L*(PDstandard4th2beta2 + PDstandard4th3beta3) + beta3L*PDstandard4th3g32; - - g33rhsL = -2*alphaL*K33L + beta1L*PDstandard4th1g33 + beta2L*PDstandard4th2g33 + - 2*(g31L*PDstandard4th3beta1 + g32L*PDstandard4th3beta2 + g33L*PDstandard4th3beta3) + beta3L*PDstandard4th3g33; - - K11rhsL = -PDstandard4th11alpha + G111*PDstandard4th1alpha + - 2*(K11L*PDstandard4th1beta1 + K21L*PDstandard4th1beta2 + K31L*PDstandard4th1beta3) + beta1L*PDstandard4th1K11 + - G211*PDstandard4th2alpha + beta2L*PDstandard4th2K11 + G311*PDstandard4th3alpha + beta3L*PDstandard4th3K11 + - alphaL*(-2*(K11L*Km11 + K21L*Km21 + K31L*Km31) + R11 + K11L*trK); - - K21rhsL = -PDstandard4th12alpha + G121*PDstandard4th1alpha + K22L*PDstandard4th1beta2 + K32L*PDstandard4th1beta3 + - beta1L*PDstandard4th1K21 + G221*PDstandard4th2alpha + K11L*PDstandard4th2beta1 + - K21L*(PDstandard4th1beta1 + PDstandard4th2beta2) + K31L*PDstandard4th2beta3 + beta2L*PDstandard4th2K21 + - G321*PDstandard4th3alpha + beta3L*PDstandard4th3K21 + - alphaL*(-2*(K11L*Km12 + K21L*Km22 + K31L*Km32) + R21 + K21L*trK); - - K31rhsL = -PDstandard4th13alpha + G131*PDstandard4th1alpha + K32L*PDstandard4th1beta2 + K33L*PDstandard4th1beta3 + - beta1L*PDstandard4th1K31 + G231*PDstandard4th2alpha + beta2L*PDstandard4th2K31 + G331*PDstandard4th3alpha + - K11L*PDstandard4th3beta1 + K21L*PDstandard4th3beta2 + K31L*(PDstandard4th1beta1 + PDstandard4th3beta3) + - beta3L*PDstandard4th3K31 + alphaL*(-2*(K11L*Km13 + K21L*Km23 + K31L*Km33) + R31 + K31L*trK); - - K22rhsL = G122*PDstandard4th1alpha + beta1L*PDstandard4th1K22 - PDstandard4th22alpha + G222*PDstandard4th2alpha + - 2*(K21L*PDstandard4th2beta1 + K22L*PDstandard4th2beta2 + K32L*PDstandard4th2beta3) + beta2L*PDstandard4th2K22 + - G322*PDstandard4th3alpha + beta3L*PDstandard4th3K22 + - alphaL*(-2*(K21L*Km12 + K22L*Km22 + K32L*Km32) + R22 + K22L*trK); - - K32rhsL = G132*PDstandard4th1alpha + beta1L*PDstandard4th1K32 - PDstandard4th23alpha + G232*PDstandard4th2alpha + - K31L*PDstandard4th2beta1 + K33L*PDstandard4th2beta3 + beta2L*PDstandard4th2K32 + G332*PDstandard4th3alpha + - K21L*PDstandard4th3beta1 + K22L*PDstandard4th3beta2 + K32L*(PDstandard4th2beta2 + PDstandard4th3beta3) + - beta3L*PDstandard4th3K32 + alphaL*(-2*(K21L*Km13 + K22L*Km23 + K32L*Km33) + R32 + K32L*trK); - - K33rhsL = G133*PDstandard4th1alpha + beta1L*PDstandard4th1K33 + G233*PDstandard4th2alpha + beta2L*PDstandard4th2K33 - - PDstandard4th33alpha + G333*PDstandard4th3alpha + - 2*(K31L*PDstandard4th3beta1 + K32L*PDstandard4th3beta2 + K33L*PDstandard4th3beta3) + beta3L*PDstandard4th3K33 + - alphaL*(-2*(K31L*Km13 + K32L*Km23 + K33L*Km33) + R33 + K33L*trK); - - alpharhsL = 0; - - beta1rhsL = 0; - - beta2rhsL = 0; - - beta3rhsL = 0; - - - /* Copy local copies back to grid functions */ - alpharhs[index] = alpharhsL; - beta1rhs[index] = beta1rhsL; - beta2rhs[index] = beta2rhsL; - beta3rhs[index] = beta3rhsL; - g11rhs[index] = g11rhsL; - g21rhs[index] = g21rhsL; - g22rhs[index] = g22rhsL; - g31rhs[index] = g31rhsL; - g32rhs[index] = g32rhsL; - g33rhs[index] = g33rhsL; - K11rhs[index] = K11rhsL; - K21rhs[index] = K21rhsL; - K22rhs[index] = K22rhsL; - K31rhs[index] = K31rhsL; - K32rhs[index] = K32rhsL; - K33rhs[index] = K33rhsL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void ADM_RHS(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &ADM_RHS_Body); -} diff --git a/ML_ADM/src/ADM_constraints.c b/ML_ADM/src/ADM_constraints.c deleted file mode 100644 index aadfe22..0000000 --- a/ML_ADM/src/ADM_constraints.c +++ /dev/null @@ -1,445 +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 ADM_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[]) -{ - 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 ADM_constraints_Body"); - } - - if (cctk_iteration % ADM_constraints_calc_every != ADM_constraints_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 detg = INITVALUE; - CCTK_REAL G111 = INITVALUE, G121 = INITVALUE, G122 = INITVALUE, G131 = INITVALUE, G132 = INITVALUE, G133 = INITVALUE; - CCTK_REAL G211 = INITVALUE, G221 = INITVALUE, G222 = INITVALUE, G231 = INITVALUE, G232 = INITVALUE, G233 = INITVALUE; - CCTK_REAL G311 = INITVALUE, G321 = INITVALUE, G322 = INITVALUE, G331 = INITVALUE, G332 = INITVALUE, G333 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE, gu33 = INITVALUE; - CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE; - CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE; - CCTK_REAL R11 = INITVALUE, R21 = INITVALUE, R22 = INITVALUE, R31 = INITVALUE, R32 = INITVALUE, R33 = INITVALUE; - CCTK_REAL trK = INITVALUE; - CCTK_REAL trR = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL g11L = INITVALUE, g21L = INITVALUE, g22L = INITVALUE, g31L = INITVALUE, g32L = INITVALUE, g33L = INITVALUE; - CCTK_REAL HL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K21L = INITVALUE, K22L = INITVALUE, K31L = INITVALUE, K32L = INITVALUE, K33L = INITVALUE; - CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - CCTK_REAL PDstandard4th1g11 = INITVALUE; - CCTK_REAL PDstandard4th2g11 = INITVALUE; - CCTK_REAL PDstandard4th3g11 = INITVALUE; - CCTK_REAL PDstandard4th22g11 = INITVALUE; - CCTK_REAL PDstandard4th33g11 = INITVALUE; - CCTK_REAL PDstandard4th12g11 = INITVALUE; - CCTK_REAL PDstandard4th13g11 = INITVALUE; - CCTK_REAL PDstandard4th21g11 = INITVALUE; - CCTK_REAL PDstandard4th23g11 = INITVALUE; - CCTK_REAL PDstandard4th31g11 = INITVALUE; - CCTK_REAL PDstandard4th32g11 = INITVALUE; - CCTK_REAL PDstandard4th1g21 = INITVALUE; - CCTK_REAL PDstandard4th2g21 = INITVALUE; - CCTK_REAL PDstandard4th3g21 = INITVALUE; - CCTK_REAL PDstandard4th33g21 = INITVALUE; - CCTK_REAL PDstandard4th12g21 = INITVALUE; - CCTK_REAL PDstandard4th13g21 = INITVALUE; - CCTK_REAL PDstandard4th21g21 = INITVALUE; - CCTK_REAL PDstandard4th23g21 = INITVALUE; - CCTK_REAL PDstandard4th31g21 = INITVALUE; - CCTK_REAL PDstandard4th32g21 = INITVALUE; - CCTK_REAL PDstandard4th1g22 = INITVALUE; - CCTK_REAL PDstandard4th2g22 = INITVALUE; - CCTK_REAL PDstandard4th3g22 = INITVALUE; - CCTK_REAL PDstandard4th11g22 = INITVALUE; - CCTK_REAL PDstandard4th33g22 = INITVALUE; - CCTK_REAL PDstandard4th12g22 = INITVALUE; - CCTK_REAL PDstandard4th13g22 = INITVALUE; - CCTK_REAL PDstandard4th21g22 = INITVALUE; - CCTK_REAL PDstandard4th23g22 = INITVALUE; - CCTK_REAL PDstandard4th31g22 = INITVALUE; - CCTK_REAL PDstandard4th32g22 = INITVALUE; - CCTK_REAL PDstandard4th1g31 = INITVALUE; - CCTK_REAL PDstandard4th2g31 = INITVALUE; - CCTK_REAL PDstandard4th3g31 = INITVALUE; - CCTK_REAL PDstandard4th22g31 = INITVALUE; - CCTK_REAL PDstandard4th12g31 = INITVALUE; - CCTK_REAL PDstandard4th13g31 = INITVALUE; - CCTK_REAL PDstandard4th21g31 = INITVALUE; - CCTK_REAL PDstandard4th23g31 = INITVALUE; - CCTK_REAL PDstandard4th31g31 = INITVALUE; - CCTK_REAL PDstandard4th32g31 = INITVALUE; - CCTK_REAL PDstandard4th1g32 = INITVALUE; - CCTK_REAL PDstandard4th2g32 = INITVALUE; - CCTK_REAL PDstandard4th3g32 = INITVALUE; - CCTK_REAL PDstandard4th11g32 = INITVALUE; - CCTK_REAL PDstandard4th12g32 = INITVALUE; - CCTK_REAL PDstandard4th13g32 = INITVALUE; - CCTK_REAL PDstandard4th21g32 = INITVALUE; - CCTK_REAL PDstandard4th23g32 = INITVALUE; - CCTK_REAL PDstandard4th31g32 = INITVALUE; - CCTK_REAL PDstandard4th32g32 = INITVALUE; - CCTK_REAL PDstandard4th1g33 = INITVALUE; - CCTK_REAL PDstandard4th2g33 = INITVALUE; - CCTK_REAL PDstandard4th3g33 = INITVALUE; - CCTK_REAL PDstandard4th11g33 = INITVALUE; - CCTK_REAL PDstandard4th22g33 = INITVALUE; - CCTK_REAL PDstandard4th12g33 = INITVALUE; - CCTK_REAL PDstandard4th13g33 = INITVALUE; - CCTK_REAL PDstandard4th21g33 = INITVALUE; - CCTK_REAL PDstandard4th23g33 = INITVALUE; - CCTK_REAL PDstandard4th31g33 = INITVALUE; - CCTK_REAL PDstandard4th32g33 = INITVALUE; - CCTK_REAL PDstandard4th2K11 = INITVALUE; - CCTK_REAL PDstandard4th3K11 = INITVALUE; - CCTK_REAL PDstandard4th1K21 = INITVALUE; - CCTK_REAL PDstandard4th2K21 = INITVALUE; - CCTK_REAL PDstandard4th3K21 = INITVALUE; - CCTK_REAL PDstandard4th1K22 = INITVALUE; - CCTK_REAL PDstandard4th3K22 = INITVALUE; - CCTK_REAL PDstandard4th1K31 = INITVALUE; - CCTK_REAL PDstandard4th2K31 = INITVALUE; - CCTK_REAL PDstandard4th3K31 = INITVALUE; - CCTK_REAL PDstandard4th1K32 = INITVALUE; - CCTK_REAL PDstandard4th2K32 = INITVALUE; - CCTK_REAL PDstandard4th3K32 = INITVALUE; - CCTK_REAL PDstandard4th1K33 = INITVALUE; - CCTK_REAL PDstandard4th2K33 = INITVALUE; - - /* Assign local copies of grid functions */ - g11L = g11[index]; - g21L = g21[index]; - g22L = g22[index]; - g31L = g31[index]; - g32L = g32[index]; - g33L = g33[index]; - K11L = K11[index]; - K21L = K21[index]; - K22L = K22[index]; - K31L = K31[index]; - K32L = K32[index]; - K33L = K33[index]; - - /* Assign local copies of subblock grid functions */ - - /* Include user supplied include files */ - - /* Precompute derivatives (new style) */ - PDstandard4th1g11 = PDstandard4th1(g11, i, j, k); - PDstandard4th2g11 = PDstandard4th2(g11, i, j, k); - PDstandard4th3g11 = PDstandard4th3(g11, i, j, k); - PDstandard4th22g11 = PDstandard4th22(g11, i, j, k); - PDstandard4th33g11 = PDstandard4th33(g11, i, j, k); - PDstandard4th23g11 = PDstandard4th23(g11, i, j, k); - PDstandard4th1g21 = PDstandard4th1(g21, i, j, k); - PDstandard4th2g21 = PDstandard4th2(g21, i, j, k); - PDstandard4th3g21 = PDstandard4th3(g21, i, j, k); - PDstandard4th33g21 = PDstandard4th33(g21, i, j, k); - PDstandard4th12g21 = PDstandard4th12(g21, i, j, k); - PDstandard4th13g21 = PDstandard4th13(g21, i, j, k); - PDstandard4th23g21 = PDstandard4th23(g21, i, j, k); - PDstandard4th1g22 = PDstandard4th1(g22, i, j, k); - PDstandard4th2g22 = PDstandard4th2(g22, i, j, k); - PDstandard4th3g22 = PDstandard4th3(g22, i, j, k); - PDstandard4th11g22 = PDstandard4th11(g22, i, j, k); - PDstandard4th33g22 = PDstandard4th33(g22, i, j, k); - PDstandard4th13g22 = PDstandard4th13(g22, i, j, k); - PDstandard4th1g31 = PDstandard4th1(g31, i, j, k); - PDstandard4th2g31 = PDstandard4th2(g31, i, j, k); - PDstandard4th3g31 = PDstandard4th3(g31, i, j, k); - PDstandard4th22g31 = PDstandard4th22(g31, i, j, k); - PDstandard4th12g31 = PDstandard4th12(g31, i, j, k); - PDstandard4th13g31 = PDstandard4th13(g31, i, j, k); - PDstandard4th23g31 = PDstandard4th23(g31, i, j, k); - PDstandard4th1g32 = PDstandard4th1(g32, i, j, k); - PDstandard4th2g32 = PDstandard4th2(g32, i, j, k); - PDstandard4th3g32 = PDstandard4th3(g32, i, j, k); - PDstandard4th11g32 = PDstandard4th11(g32, i, j, k); - PDstandard4th12g32 = PDstandard4th12(g32, i, j, k); - PDstandard4th13g32 = PDstandard4th13(g32, i, j, k); - PDstandard4th23g32 = PDstandard4th23(g32, i, j, k); - PDstandard4th1g33 = PDstandard4th1(g33, i, j, k); - PDstandard4th2g33 = PDstandard4th2(g33, i, j, k); - PDstandard4th3g33 = PDstandard4th3(g33, i, j, k); - PDstandard4th11g33 = PDstandard4th11(g33, i, j, k); - PDstandard4th22g33 = PDstandard4th22(g33, i, j, k); - PDstandard4th12g33 = PDstandard4th12(g33, i, j, k); - PDstandard4th2K11 = PDstandard4th2(K11, i, j, k); - PDstandard4th3K11 = PDstandard4th3(K11, i, j, k); - PDstandard4th1K21 = PDstandard4th1(K21, i, j, k); - PDstandard4th2K21 = PDstandard4th2(K21, i, j, k); - PDstandard4th3K21 = PDstandard4th3(K21, i, j, k); - PDstandard4th1K22 = PDstandard4th1(K22, i, j, k); - PDstandard4th3K22 = PDstandard4th3(K22, i, j, k); - PDstandard4th1K31 = PDstandard4th1(K31, i, j, k); - PDstandard4th2K31 = PDstandard4th2(K31, i, j, k); - PDstandard4th3K31 = PDstandard4th3(K31, i, j, k); - PDstandard4th1K32 = PDstandard4th1(K32, i, j, k); - PDstandard4th2K32 = PDstandard4th2(K32, i, j, k); - PDstandard4th3K32 = PDstandard4th3(K32, i, j, k); - PDstandard4th1K33 = PDstandard4th1(K33, i, j, k); - PDstandard4th2K33 = PDstandard4th2(K33, i, j, k); - - /* Precompute derivatives (old style) */ - - /* Calculate temporaries and grid functions */ - detg = 2*g21L*g31L*g32L + g33L*(g11L*g22L - SQR(g21L)) - g22L*SQR(g31L) - g11L*SQR(g32L); - - gu11 = INV(detg)*(g22L*g33L - SQR(g32L)); - - gu12 = (g31L*g32L - g21L*g33L)*INV(detg); - - gu13 = (-(g22L*g31L) + g21L*g32L)*INV(detg); - - gu22 = INV(detg)*(g11L*g33L - SQR(g31L)); - - gu23 = (g21L*g31L - g11L*g32L)*INV(detg); - - gu33 = INV(detg)*(g11L*g22L - SQR(g21L)); - - G111 = khalf*(gu11*PDstandard4th1g11 + 2*(gu12*PDstandard4th1g21 + gu13*PDstandard4th1g31) - gu12*PDstandard4th2g11 - - gu13*PDstandard4th3g11); - - G211 = khalf*(gu12*PDstandard4th1g11 + 2*(gu22*PDstandard4th1g21 + gu23*PDstandard4th1g31) - gu22*PDstandard4th2g11 - - gu23*PDstandard4th3g11); - - G311 = khalf*(gu13*PDstandard4th1g11 + 2*(gu23*PDstandard4th1g21 + gu33*PDstandard4th1g31) - gu23*PDstandard4th2g11 - - gu33*PDstandard4th3g11); - - G121 = khalf*(gu12*PDstandard4th1g22 + gu11*PDstandard4th2g11 + - gu13*(PDstandard4th1g32 + PDstandard4th2g31 - PDstandard4th3g21)); - - G221 = khalf*(gu22*PDstandard4th1g22 + gu12*PDstandard4th2g11 + - gu23*(PDstandard4th1g32 + PDstandard4th2g31 - PDstandard4th3g21)); - - G321 = khalf*(gu23*PDstandard4th1g22 + gu13*PDstandard4th2g11 + - gu33*(PDstandard4th1g32 + PDstandard4th2g31 - PDstandard4th3g21)); - - G131 = khalf*(gu13*PDstandard4th1g33 + gu11*PDstandard4th3g11 + - gu12*(PDstandard4th1g32 - PDstandard4th2g31 + PDstandard4th3g21)); - - G231 = khalf*(gu23*PDstandard4th1g33 + gu12*PDstandard4th3g11 + - gu22*(PDstandard4th1g32 - PDstandard4th2g31 + PDstandard4th3g21)); - - G331 = khalf*(gu33*PDstandard4th1g33 + gu13*PDstandard4th3g11 + - gu23*(PDstandard4th1g32 - PDstandard4th2g31 + PDstandard4th3g21)); - - G122 = khalf*(gu11*(-PDstandard4th1g22 + 2*PDstandard4th2g21) + gu12*PDstandard4th2g22 + - gu13*(2*PDstandard4th2g32 - PDstandard4th3g22)); - - G222 = khalf*(gu12*(-PDstandard4th1g22 + 2*PDstandard4th2g21) + gu22*PDstandard4th2g22 + - gu23*(2*PDstandard4th2g32 - PDstandard4th3g22)); - - G322 = khalf*(gu13*(-PDstandard4th1g22 + 2*PDstandard4th2g21) + gu23*PDstandard4th2g22 + - gu33*(2*PDstandard4th2g32 - PDstandard4th3g22)); - - G132 = khalf*(gu13*PDstandard4th2g33 + gu11*(-PDstandard4th1g32 + PDstandard4th2g31 + PDstandard4th3g21) + - gu12*PDstandard4th3g22); - - G232 = khalf*(gu23*PDstandard4th2g33 + gu12*(-PDstandard4th1g32 + PDstandard4th2g31 + PDstandard4th3g21) + - gu22*PDstandard4th3g22); - - G332 = khalf*(gu33*PDstandard4th2g33 + gu13*(-PDstandard4th1g32 + PDstandard4th2g31 + PDstandard4th3g21) + - gu23*PDstandard4th3g22); - - G133 = khalf*(-(gu11*PDstandard4th1g33) - gu12*PDstandard4th2g33 + 2*gu11*PDstandard4th3g31 + - 2*gu12*PDstandard4th3g32 + gu13*PDstandard4th3g33); - - G233 = khalf*(-(gu12*PDstandard4th1g33) - gu22*PDstandard4th2g33 + 2*gu12*PDstandard4th3g31 + - 2*gu22*PDstandard4th3g32 + gu23*PDstandard4th3g33); - - G333 = khalf*(-(gu13*PDstandard4th1g33) - gu23*PDstandard4th2g33 + 2*gu13*PDstandard4th3g31 + - 2*gu23*PDstandard4th3g32 + gu33*PDstandard4th3g33); - - R11 = khalf*(-(gu22*PDstandard4th11g22) - 2*(G111*G122 + G111*G133 + G211*G222 + G211*G233 + G311*G322 + G311*G333 + - gu23*PDstandard4th11g32) - gu33*PDstandard4th11g33 + 2*gu22*PDstandard4th12g21 + 2*gu23*PDstandard4th12g31 + - 2*gu23*PDstandard4th13g21 + 2*gu33*PDstandard4th13g31 - gu22*PDstandard4th22g11 - 2*gu23*PDstandard4th23g11 - - gu33*PDstandard4th33g11 + 2*SQR(G121) + 2*SQR(G131) + 2*SQR(G221) + 2*SQR(G231) + 2*SQR(G321) + 2*SQR(G331)); - - R21 = khalf*(2*(G131*G132 + G231*G232 + G331*G332) - 2*(G121*G133 + G221*G233 + G321*G333 + gu12*PDstandard4th12g21) - - gu23*PDstandard4th12g32 - gu33*PDstandard4th12g33 + - gu13*(PDstandard4th11g32 - PDstandard4th12g31 - PDstandard4th13g21) + gu23*PDstandard4th13g22 + - gu33*PDstandard4th13g32 + gu12*(PDstandard4th11g22 + PDstandard4th22g11) + gu23*PDstandard4th22g31 + - gu13*PDstandard4th23g11 - gu23*PDstandard4th23g21 + gu33*PDstandard4th23g31 - gu33*PDstandard4th33g21); - - R31 = khalf*(2*(G121*G132 + G221*G232 + G321*G332) - 2*(G122*G131 + G222*G231 + G322*G331 + gu13*PDstandard4th13g31) + - gu12*(PDstandard4th11g32 - PDstandard4th12g31 - PDstandard4th13g21 + PDstandard4th23g11) + - gu22*(PDstandard4th12g32 - PDstandard4th13g22 - PDstandard4th22g31 + PDstandard4th23g21) + - gu13*(PDstandard4th11g33 + PDstandard4th33g11) + - gu23*(PDstandard4th12g33 - PDstandard4th13g32 - PDstandard4th23g31 + PDstandard4th33g21)); - - R22 = khalf*(-2*(G122*(G111 + G133) + G222*(G211 + G233) + G322*(G311 + G333) + gu13*PDstandard4th13g22) + - gu11*(-PDstandard4th11g22 + 2*PDstandard4th12g21 - PDstandard4th22g11) + - gu13*(-2*PDstandard4th22g31 + 2*(PDstandard4th12g32 + PDstandard4th23g21)) + - gu33*(-PDstandard4th22g33 + 2*PDstandard4th23g32 - PDstandard4th33g22) + - 2*(SQR(G121) + SQR(G132) + SQR(G221) + SQR(G232) + SQR(G321) + SQR(G332))); - - R32 = khalf*(2*(G121*G131 + G221*G231 + G321*G331) + - gu11*(-PDstandard4th11g32 + PDstandard4th12g31 + PDstandard4th13g21 - PDstandard4th23g11) + - gu12*(-PDstandard4th12g32 + PDstandard4th13g22 + PDstandard4th22g31 - PDstandard4th23g21) - - 2*(G111*G132 + G211*G232 + G311*G332 + gu23*PDstandard4th23g32) + - gu13*(PDstandard4th12g33 - PDstandard4th13g32 - PDstandard4th23g31 + PDstandard4th33g21) + - gu23*(PDstandard4th22g33 + PDstandard4th33g22)); - - R33 = khalf*(gu11*(-PDstandard4th11g33 + 2*PDstandard4th13g31 - PDstandard4th33g11) - - 2*((G111 + G122)*G133 + (G211 + G222)*G233 + (G311 + G322)*G333 + - gu12*(PDstandard4th12g33 + PDstandard4th33g21)) + - gu22*(-PDstandard4th22g33 + 2*PDstandard4th23g32 - PDstandard4th33g22) + - 2*(gu12*(PDstandard4th13g32 + PDstandard4th23g31) + SQR(G131) + SQR(G132) + SQR(G231) + SQR(G232) + SQR(G331) + - SQR(G332))); - - trR = gu11*R11 + gu22*R22 + 2*(gu12*R21 + gu13*R31 + gu23*R32) + gu33*R33; - - Km11 = gu11*K11L + gu12*K21L + gu13*K31L; - - Km21 = gu12*K11L + gu22*K21L + gu23*K31L; - - Km31 = gu13*K11L + gu23*K21L + gu33*K31L; - - Km12 = gu11*K21L + gu12*K22L + gu13*K32L; - - Km22 = gu12*K21L + gu22*K22L + gu23*K32L; - - Km32 = gu13*K21L + gu23*K22L + gu33*K32L; - - Km13 = gu11*K31L + gu12*K32L + gu13*K33L; - - Km23 = gu12*K31L + gu22*K32L + gu23*K33L; - - Km33 = gu13*K31L + gu23*K32L + gu33*K33L; - - trK = Km11 + Km22 + Km33; - - HL = -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - SQR(Km33) + SQR(trK); - - M1L = gu12*(-(G121*K11L) + G111*K21L - G221*K21L + G211*K22L - G321*K31L + G311*K32L - PDstandard4th1K21 + - PDstandard4th2K11) + gu22*(-(G122*K11L) + G121*K21L - G222*K21L + G221*K22L - G322*K31L + G321*K32L - - PDstandard4th1K22 + PDstandard4th2K21) + gu13* - (-(G131*K11L) - G231*K21L + G111*K31L - G331*K31L + G211*K32L + G311*K33L - PDstandard4th1K31 + PDstandard4th3K11)\ - + gu23*(G131*K21L + G231*K22L + G121*K31L + (G221 + G331)*K32L + G321*K33L - - 2*(G132*K11L + G232*K21L + G332*K31L + PDstandard4th1K32) + PDstandard4th2K31 + PDstandard4th3K21) + - gu33*(-(G133*K11L) - G233*K21L + G131*K31L - G333*K31L + G231*K32L + G331*K33L - PDstandard4th1K33 + - PDstandard4th3K31); - - M2L = gu11*(G121*K11L + (-G111 + G221)*K21L - G211*K22L + G321*K31L - G311*K32L + PDstandard4th1K21 - - PDstandard4th2K11) + gu12*(G122*K11L + (-G121 + G222)*K21L - G221*K22L + G322*K31L - G321*K32L + - PDstandard4th1K22 - PDstandard4th2K21) + gu13* - (G132*K11L + (-2*G131 + G232)*K21L + (G121 + G332)*K31L + G221*K32L + G321*K33L + PDstandard4th1K32 - - 2*(G231*K22L + G331*K32L + PDstandard4th2K31) + PDstandard4th3K21) + - gu23*(-(G132*K21L) - G232*K22L + G122*K31L + G222*K32L - G332*K32L + G322*K33L - PDstandard4th2K32 + - PDstandard4th3K22) + gu33*(-(G133*K21L) - G233*K22L + G132*K31L + G232*K32L - G333*K32L + G332*K33L - - PDstandard4th2K33 + PDstandard4th3K32); - - M3L = gu11*(G131*K11L + G231*K21L + (-G111 + G331)*K31L - G211*K32L - G311*K33L + PDstandard4th1K31 - - PDstandard4th3K11) + gu12*(G132*K11L + (G131 + G232)*K21L + G231*K22L + (-2*G121 + G332)*K31L + - (-2*G221 + G331)*K32L + PDstandard4th1K32 + PDstandard4th2K31 - 2*(G321*K33L + PDstandard4th3K21)) + - gu22*(G132*K21L + G232*K22L - G122*K31L - G222*K32L + G332*K32L - G322*K33L + PDstandard4th2K32 - - PDstandard4th3K22) + gu13*(G133*K11L + G233*K21L + (-G131 + G333)*K31L - G231*K32L - G331*K33L + - PDstandard4th1K33 - PDstandard4th3K31) + gu23* - (G133*K21L + G233*K22L - G132*K31L - G232*K32L + G333*K32L - G332*K33L + PDstandard4th2K33 - PDstandard4th3K32); - - - /* 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 (somename); -} - -void ADM_constraints(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &ADM_constraints_Body); -} diff --git a/ML_ADM/src/ADM_convertFromADMBase.c b/ML_ADM/src/ADM_convertFromADMBase.c deleted file mode 100644 index c3275a4..0000000 --- a/ML_ADM/src/ADM_convertFromADMBase.c +++ /dev/null @@ -1,212 +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 ADM_convertFromADMBase_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 ADM_convertFromADMBase_Body"); - } - - if (cctk_iteration % ADM_convertFromADMBase_calc_every != ADM_convertFromADMBase_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 alpL = INITVALUE; - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL betaxL = INITVALUE; - CCTK_REAL betayL = INITVALUE; - CCTK_REAL betazL = INITVALUE; - CCTK_REAL g11L = INITVALUE, g21L = INITVALUE, g22L = INITVALUE, g31L = INITVALUE, g32L = INITVALUE, g33L = INITVALUE; - CCTK_REAL gxxL = INITVALUE; - CCTK_REAL gxyL = INITVALUE; - CCTK_REAL gxzL = INITVALUE; - CCTK_REAL gyyL = INITVALUE; - CCTK_REAL gyzL = INITVALUE; - CCTK_REAL gzzL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K21L = INITVALUE, K22L = INITVALUE, K31L = INITVALUE, K32L = INITVALUE, K33L = INITVALUE; - CCTK_REAL kxxL = INITVALUE; - CCTK_REAL kxyL = INITVALUE; - CCTK_REAL kxzL = INITVALUE; - CCTK_REAL kyyL = INITVALUE; - CCTK_REAL kyzL = INITVALUE; - CCTK_REAL kzzL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - alpL = alp[index]; - betaxL = betax[index]; - betayL = betay[index]; - betazL = betaz[index]; - gxxL = gxx[index]; - gxyL = gxy[index]; - gxzL = gxz[index]; - gyyL = gyy[index]; - gyzL = gyz[index]; - gzzL = gzz[index]; - kxxL = kxx[index]; - kxyL = kxy[index]; - kxzL = kxz[index]; - kyyL = kyy[index]; - kyzL = kyz[index]; - kzzL = kzz[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 */ - g11L = gxxL; - - g21L = gxyL; - - g31L = gxzL; - - g22L = gyyL; - - g32L = gyzL; - - g33L = gzzL; - - K11L = kxxL; - - K21L = kxyL; - - K31L = kxzL; - - K22L = kyyL; - - K32L = kyzL; - - K33L = kzzL; - - alphaL = alpL; - - beta1L = betaxL; - - beta2L = betayL; - - beta3L = betazL; - - - /* Copy local copies back to grid functions */ - alpha[index] = alphaL; - beta1[index] = beta1L; - beta2[index] = beta2L; - beta3[index] = beta3L; - g11[index] = g11L; - g21[index] = g21L; - g22[index] = g22L; - g31[index] = g31L; - g32[index] = g32L; - g33[index] = g33L; - K11[index] = K11L; - K21[index] = K21L; - K22[index] = K22L; - K31[index] = K31L; - K32[index] = K32L; - K33[index] = K33L; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void ADM_convertFromADMBase(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &ADM_convertFromADMBase_Body); -} diff --git a/ML_ADM/src/ADM_convertToADMBase.c b/ML_ADM/src/ADM_convertToADMBase.c deleted file mode 100644 index 77d9852..0000000 --- a/ML_ADM/src/ADM_convertToADMBase.c +++ /dev/null @@ -1,228 +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 ADM_convertToADMBase_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 ADM_convertToADMBase_Body"); - } - - if (cctk_iteration % ADM_convertToADMBase_calc_every != ADM_convertToADMBase_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 alpL = INITVALUE; - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL betaxL = INITVALUE; - CCTK_REAL betayL = INITVALUE; - CCTK_REAL betazL = INITVALUE; - CCTK_REAL dtalpL = INITVALUE; - CCTK_REAL dtbetaxL = INITVALUE; - CCTK_REAL dtbetayL = INITVALUE; - CCTK_REAL dtbetazL = INITVALUE; - CCTK_REAL g11L = INITVALUE, g21L = INITVALUE, g22L = INITVALUE, g31L = INITVALUE, g32L = INITVALUE, g33L = INITVALUE; - CCTK_REAL gxxL = INITVALUE; - CCTK_REAL gxyL = INITVALUE; - CCTK_REAL gxzL = INITVALUE; - CCTK_REAL gyyL = INITVALUE; - CCTK_REAL gyzL = INITVALUE; - CCTK_REAL gzzL = INITVALUE; - CCTK_REAL K11L = INITVALUE, K21L = INITVALUE, K22L = INITVALUE, K31L = INITVALUE, K32L = INITVALUE, K33L = INITVALUE; - CCTK_REAL kxxL = INITVALUE; - CCTK_REAL kxyL = INITVALUE; - CCTK_REAL kxzL = INITVALUE; - CCTK_REAL kyyL = INITVALUE; - CCTK_REAL kyzL = INITVALUE; - CCTK_REAL kzzL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - alphaL = alpha[index]; - beta1L = beta1[index]; - beta2L = beta2[index]; - beta3L = beta3[index]; - g11L = g11[index]; - g21L = g21[index]; - g22L = g22[index]; - g31L = g31[index]; - g32L = g32[index]; - g33L = g33[index]; - K11L = K11[index]; - K21L = K21[index]; - K22L = K22[index]; - K31L = K31[index]; - K32L = K32[index]; - K33L = K33[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 */ - gxxL = g11L; - - gxyL = g21L; - - gxzL = g31L; - - gyyL = g22L; - - gyzL = g32L; - - gzzL = g33L; - - kxxL = K11L; - - kxyL = K21L; - - kxzL = K31L; - - kyyL = K22L; - - kyzL = K32L; - - kzzL = K33L; - - alpL = alphaL; - - dtalpL = 0; - - betaxL = beta1L; - - betayL = beta2L; - - betazL = beta3L; - - dtbetaxL = 0; - - dtbetayL = 0; - - dtbetazL = 0; - - - /* Copy local copies back to grid functions */ - alp[index] = alpL; - betax[index] = betaxL; - betay[index] = betayL; - betaz[index] = betazL; - dtalp[index] = dtalpL; - dtbetax[index] = dtbetaxL; - dtbetay[index] = dtbetayL; - dtbetaz[index] = dtbetazL; - gxx[index] = gxxL; - gxy[index] = gxyL; - gxz[index] = gxzL; - gyy[index] = gyyL; - gyz[index] = gyzL; - gzz[index] = gzzL; - kxx[index] = kxxL; - kxy[index] = kxyL; - kxz[index] = kxzL; - kyy[index] = kyyL; - kyz[index] = kyzL; - kzz[index] = kzzL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void ADM_convertToADMBase(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &ADM_convertToADMBase_Body); -} diff --git a/ML_BSSN/src/BSSN_Minkowski.c b/ML_BSSN/src/BSSN_Minkowski.c deleted file mode 100644 index 85b1bc6..0000000 --- a/ML_BSSN/src/BSSN_Minkowski.c +++ /dev/null @@ -1,212 +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 BSSN_Minkowski_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 BSSN_Minkowski_Body"); - } - - if (cctk_iteration % BSSN_Minkowski_calc_every != BSSN_Minkowski_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 At11L = INITVALUE, At21L = INITVALUE, At22L = INITVALUE, At31L = INITVALUE, At32L = INITVALUE, At33L = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL dtalphaL = INITVALUE; - CCTK_REAL dtbeta1L = INITVALUE, dtbeta2L = INITVALUE, dtbeta3L = INITVALUE; - CCTK_REAL gt11L = INITVALUE, gt21L = INITVALUE, gt22L = INITVALUE, gt31L = INITVALUE, gt32L = 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 */ - - /* 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 */ - phiL = 0; - - gt11L = 1; - - gt21L = 0; - - gt31L = 0; - - gt22L = 1; - - gt32L = 0; - - gt33L = 1; - - trKL = 0; - - At11L = 0; - - At21L = 0; - - At31L = 0; - - At22L = 0; - - At32L = 0; - - At33L = 0; - - Xt1L = 0; - - Xt2L = 0; - - Xt3L = 0; - - alphaL = 1; - - dtalphaL = 0; - - beta1L = 0; - - beta2L = 0; - - beta3L = 0; - - dtbeta1L = 0; - - dtbeta2L = 0; - - dtbeta3L = 0; - - - /* Copy local copies back to grid functions */ - alpha[index] = alphaL; - At11[index] = At11L; - At21[index] = At21L; - At22[index] = At22L; - At31[index] = At31L; - At32[index] = At32L; - At33[index] = At33L; - beta1[index] = beta1L; - beta2[index] = beta2L; - beta3[index] = beta3L; - dtalpha[index] = dtalphaL; - dtbeta1[index] = dtbeta1L; - dtbeta2[index] = dtbeta2L; - dtbeta3[index] = dtbeta3L; - gt11[index] = gt11L; - gt21[index] = gt21L; - gt22[index] = gt22L; - gt31[index] = gt31L; - gt32[index] = gt32L; - 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 (somename); -} - -void BSSN_Minkowski(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &BSSN_Minkowski_Body); -} diff --git a/ML_BSSN/src/BSSN_constraints.c b/ML_BSSN/src/BSSN_constraints.c deleted file mode 100644 index 9378886..0000000 --- a/ML_BSSN/src/BSSN_constraints.c +++ /dev/null @@ -1,980 +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 BSSN_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[]) -{ - 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 BSSN_constraints_Body"); - } - - if (cctk_iteration % BSSN_constraints_calc_every != BSSN_constraints_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 ddetg1 = INITVALUE, ddetg2 = INITVALUE, ddetg3 = INITVALUE; - CCTK_REAL ddetgt1 = INITVALUE, ddetgt2 = INITVALUE, ddetgt3 = INITVALUE; - CCTK_REAL detg = INITVALUE; - CCTK_REAL detgt = INITVALUE; - CCTK_REAL e4phi = INITVALUE; - CCTK_REAL g11 = INITVALUE; - CCTK_REAL G111 = INITVALUE, G121 = INITVALUE, G122 = INITVALUE, G131 = INITVALUE, G132 = INITVALUE, G133 = INITVALUE; - CCTK_REAL g21 = INITVALUE; - CCTK_REAL G211 = INITVALUE; - CCTK_REAL g22 = INITVALUE; - CCTK_REAL G221 = INITVALUE, G222 = INITVALUE, G231 = INITVALUE, G232 = INITVALUE, G233 = INITVALUE; - CCTK_REAL g31 = INITVALUE; - CCTK_REAL G311 = INITVALUE; - CCTK_REAL g32 = INITVALUE; - CCTK_REAL G321 = INITVALUE, G322 = INITVALUE; - CCTK_REAL g33 = INITVALUE; - CCTK_REAL G331 = INITVALUE, G332 = INITVALUE, G333 = INITVALUE; - CCTK_REAL gK112 = INITVALUE, gK113 = INITVALUE, gK211 = INITVALUE, gK212 = INITVALUE, gK213 = INITVALUE, gK221 = INITVALUE; - CCTK_REAL gK223 = INITVALUE, gK311 = INITVALUE, gK312 = INITVALUE, gK313 = INITVALUE, gK321 = INITVALUE, gK322 = INITVALUE; - CCTK_REAL gK323 = INITVALUE, gK331 = INITVALUE, gK332 = INITVALUE; - CCTK_REAL Gt111 = INITVALUE, Gt121 = INITVALUE, Gt122 = INITVALUE, Gt131 = INITVALUE, Gt132 = INITVALUE, Gt133 = INITVALUE; - CCTK_REAL Gt211 = INITVALUE, Gt221 = INITVALUE, Gt222 = INITVALUE, Gt231 = INITVALUE, Gt232 = INITVALUE, Gt233 = INITVALUE; - CCTK_REAL Gt311 = INITVALUE, Gt321 = INITVALUE, Gt322 = INITVALUE, Gt331 = INITVALUE, Gt332 = INITVALUE, Gt333 = INITVALUE; - CCTK_REAL gtu11 = INITVALUE, gtu12 = INITVALUE, gtu13 = INITVALUE, gtu22 = INITVALUE, gtu23 = INITVALUE, gtu33 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE, gu33 = INITVALUE; - CCTK_REAL K11 = INITVALUE, K21 = INITVALUE, K22 = INITVALUE, K31 = INITVALUE, K32 = INITVALUE, K33 = INITVALUE; - CCTK_REAL Km11 = INITVALUE, Km12 = INITVALUE, Km13 = INITVALUE, Km21 = INITVALUE, Km22 = INITVALUE, Km23 = INITVALUE; - CCTK_REAL Km31 = INITVALUE, Km32 = INITVALUE, Km33 = INITVALUE; - CCTK_REAL R11 = INITVALUE, R21 = INITVALUE, R22 = INITVALUE, R31 = INITVALUE, R32 = INITVALUE, R33 = INITVALUE; - CCTK_REAL Rphi11 = INITVALUE, Rphi21 = INITVALUE, Rphi22 = INITVALUE, Rphi31 = INITVALUE, Rphi32 = INITVALUE, Rphi33 = INITVALUE; - CCTK_REAL Rt11 = INITVALUE, Rt21 = INITVALUE, Rt22 = INITVALUE, Rt31 = INITVALUE, Rt32 = INITVALUE, Rt33 = INITVALUE; - CCTK_REAL trR = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL At11L = INITVALUE, At21L = INITVALUE, At22L = INITVALUE, At31L = INITVALUE, At32L = INITVALUE, At33L = INITVALUE; - CCTK_REAL cAL = INITVALUE; - CCTK_REAL cSL = INITVALUE; - CCTK_REAL cXt1L = INITVALUE, cXt2L = INITVALUE, cXt3L = INITVALUE; - CCTK_REAL gt11L = INITVALUE, gt21L = INITVALUE, gt22L = INITVALUE, gt31L = INITVALUE, gt32L = INITVALUE, gt33L = INITVALUE; - CCTK_REAL HL = INITVALUE; - CCTK_REAL M1L = INITVALUE, M2L = INITVALUE, M3L = INITVALUE; - CCTK_REAL phiL = INITVALUE; - CCTK_REAL trKL = INITVALUE; - CCTK_REAL Xt1L = INITVALUE, Xt2L = INITVALUE, Xt3L = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - CCTK_REAL PDstandard4th2At11 = INITVALUE; - CCTK_REAL PDstandard4th3At11 = INITVALUE; - CCTK_REAL PDstandard4th1At21 = INITVALUE; - CCTK_REAL PDstandard4th2At21 = INITVALUE; - CCTK_REAL PDstandard4th3At21 = INITVALUE; - CCTK_REAL PDstandard4th1At22 = INITVALUE; - CCTK_REAL PDstandard4th3At22 = INITVALUE; - CCTK_REAL PDstandard4th1At31 = INITVALUE; - CCTK_REAL PDstandard4th2At31 = INITVALUE; - CCTK_REAL PDstandard4th3At31 = INITVALUE; - CCTK_REAL PDstandard4th1At32 = INITVALUE; - CCTK_REAL PDstandard4th2At32 = INITVALUE; - CCTK_REAL PDstandard4th3At32 = INITVALUE; - CCTK_REAL PDstandard4th1At33 = INITVALUE; - CCTK_REAL PDstandard4th2At33 = INITVALUE; - CCTK_REAL PDstandard4th1gt11 = INITVALUE; - CCTK_REAL PDstandard4th2gt11 = INITVALUE; - CCTK_REAL PDstandard4th3gt11 = INITVALUE; - CCTK_REAL PDstandard4th11gt11 = INITVALUE; - CCTK_REAL PDstandard4th22gt11 = INITVALUE; - CCTK_REAL PDstandard4th33gt11 = INITVALUE; - CCTK_REAL PDstandard4th12gt11 = INITVALUE; - CCTK_REAL PDstandard4th13gt11 = INITVALUE; - CCTK_REAL PDstandard4th21gt11 = INITVALUE; - CCTK_REAL PDstandard4th23gt11 = INITVALUE; - CCTK_REAL PDstandard4th31gt11 = INITVALUE; - CCTK_REAL PDstandard4th32gt11 = INITVALUE; - CCTK_REAL PDstandard4th1gt21 = INITVALUE; - CCTK_REAL PDstandard4th2gt21 = INITVALUE; - CCTK_REAL PDstandard4th3gt21 = INITVALUE; - CCTK_REAL PDstandard4th11gt21 = INITVALUE; - CCTK_REAL PDstandard4th22gt21 = INITVALUE; - CCTK_REAL PDstandard4th33gt21 = INITVALUE; - CCTK_REAL PDstandard4th12gt21 = INITVALUE; - CCTK_REAL PDstandard4th13gt21 = INITVALUE; - CCTK_REAL PDstandard4th21gt21 = INITVALUE; - CCTK_REAL PDstandard4th23gt21 = INITVALUE; - CCTK_REAL PDstandard4th31gt21 = INITVALUE; - CCTK_REAL PDstandard4th32gt21 = INITVALUE; - CCTK_REAL PDstandard4th1gt22 = INITVALUE; - CCTK_REAL PDstandard4th2gt22 = INITVALUE; - CCTK_REAL PDstandard4th3gt22 = INITVALUE; - CCTK_REAL PDstandard4th11gt22 = INITVALUE; - CCTK_REAL PDstandard4th22gt22 = INITVALUE; - CCTK_REAL PDstandard4th33gt22 = INITVALUE; - CCTK_REAL PDstandard4th12gt22 = INITVALUE; - CCTK_REAL PDstandard4th13gt22 = INITVALUE; - CCTK_REAL PDstandard4th21gt22 = INITVALUE; - CCTK_REAL PDstandard4th23gt22 = INITVALUE; - CCTK_REAL PDstandard4th31gt22 = INITVALUE; - CCTK_REAL PDstandard4th32gt22 = INITVALUE; - CCTK_REAL PDstandard4th1gt31 = INITVALUE; - CCTK_REAL PDstandard4th2gt31 = INITVALUE; - CCTK_REAL PDstandard4th3gt31 = INITVALUE; - CCTK_REAL PDstandard4th11gt31 = INITVALUE; - CCTK_REAL PDstandard4th22gt31 = INITVALUE; - CCTK_REAL PDstandard4th33gt31 = INITVALUE; - CCTK_REAL PDstandard4th12gt31 = INITVALUE; - CCTK_REAL PDstandard4th13gt31 = INITVALUE; - CCTK_REAL PDstandard4th21gt31 = INITVALUE; - CCTK_REAL PDstandard4th23gt31 = INITVALUE; - CCTK_REAL PDstandard4th31gt31 = INITVALUE; - CCTK_REAL PDstandard4th32gt31 = INITVALUE; - CCTK_REAL PDstandard4th1gt32 = INITVALUE; - CCTK_REAL PDstandard4th2gt32 = INITVALUE; - CCTK_REAL PDstandard4th3gt32 = INITVALUE; - CCTK_REAL PDstandard4th11gt32 = INITVALUE; - CCTK_REAL PDstandard4th22gt32 = INITVALUE; - CCTK_REAL PDstandard4th33gt32 = INITVALUE; - CCTK_REAL PDstandard4th12gt32 = INITVALUE; - CCTK_REAL PDstandard4th13gt32 = INITVALUE; - CCTK_REAL PDstandard4th21gt32 = INITVALUE; - CCTK_REAL PDstandard4th23gt32 = INITVALUE; - CCTK_REAL PDstandard4th31gt32 = INITVALUE; - CCTK_REAL PDstandard4th32gt32 = INITVALUE; - CCTK_REAL PDstandard4th1gt33 = INITVALUE; - CCTK_REAL PDstandard4th2gt33 = INITVALUE; - CCTK_REAL PDstandard4th3gt33 = INITVALUE; - CCTK_REAL PDstandard4th11gt33 = INITVALUE; - CCTK_REAL PDstandard4th22gt33 = INITVALUE; - CCTK_REAL PDstandard4th33gt33 = INITVALUE; - CCTK_REAL PDstandard4th12gt33 = INITVALUE; - CCTK_REAL PDstandard4th13gt33 = INITVALUE; - CCTK_REAL PDstandard4th21gt33 = INITVALUE; - CCTK_REAL PDstandard4th23gt33 = INITVALUE; - CCTK_REAL PDstandard4th31gt33 = INITVALUE; - CCTK_REAL PDstandard4th32gt33 = INITVALUE; - CCTK_REAL PDstandard4th1phi = INITVALUE; - CCTK_REAL PDstandard4th2phi = INITVALUE; - CCTK_REAL PDstandard4th3phi = INITVALUE; - CCTK_REAL PDstandard4th11phi = INITVALUE; - CCTK_REAL PDstandard4th22phi = INITVALUE; - CCTK_REAL PDstandard4th33phi = INITVALUE; - CCTK_REAL PDstandard4th12phi = INITVALUE; - CCTK_REAL PDstandard4th13phi = INITVALUE; - CCTK_REAL PDstandard4th21phi = INITVALUE; - CCTK_REAL PDstandard4th23phi = INITVALUE; - CCTK_REAL PDstandard4th31phi = INITVALUE; - CCTK_REAL PDstandard4th32phi = INITVALUE; - CCTK_REAL PDstandard4th1trK = INITVALUE; - CCTK_REAL PDstandard4th2trK = INITVALUE; - CCTK_REAL PDstandard4th3trK = INITVALUE; - CCTK_REAL PDstandard4th1Xt1 = INITVALUE; - CCTK_REAL PDstandard4th2Xt1 = INITVALUE; - CCTK_REAL PDstandard4th3Xt1 = INITVALUE; - CCTK_REAL PDstandard4th1Xt2 = INITVALUE; - CCTK_REAL PDstandard4th2Xt2 = INITVALUE; - CCTK_REAL PDstandard4th3Xt2 = INITVALUE; - CCTK_REAL PDstandard4th1Xt3 = INITVALUE; - CCTK_REAL PDstandard4th2Xt3 = INITVALUE; - CCTK_REAL PDstandard4th3Xt3 = INITVALUE; - - /* Assign local copies of grid functions */ - At11L = At11[index]; - At21L = At21[index]; - At22L = At22[index]; - At31L = At31[index]; - At32L = At32[index]; - At33L = At33[index]; - gt11L = gt11[index]; - gt21L = gt21[index]; - gt22L = gt22[index]; - gt31L = gt31[index]; - gt32L = gt32[index]; - gt33L = gt33[index]; - phiL = phi[index]; - trKL = trK[index]; - Xt1L = Xt1[index]; - Xt2L = Xt2[index]; - Xt3L = Xt3[index]; - - /* Assign local copies of subblock grid functions */ - - /* Include user supplied include files */ - - /* Precompute derivatives (new style) */ - PDstandard4th2At11 = PDstandard4th2(At11, i, j, k); - PDstandard4th3At11 = PDstandard4th3(At11, i, j, k); - PDstandard4th1At21 = PDstandard4th1(At21, i, j, k); - PDstandard4th2At21 = PDstandard4th2(At21, i, j, k); - PDstandard4th3At21 = PDstandard4th3(At21, i, j, k); - PDstandard4th1At22 = PDstandard4th1(At22, i, j, k); - PDstandard4th3At22 = PDstandard4th3(At22, i, j, k); - PDstandard4th1At31 = PDstandard4th1(At31, i, j, k); - PDstandard4th2At31 = PDstandard4th2(At31, i, j, k); - PDstandard4th3At31 = PDstandard4th3(At31, i, j, k); - PDstandard4th1At32 = PDstandard4th1(At32, i, j, k); - PDstandard4th2At32 = PDstandard4th2(At32, i, j, k); - PDstandard4th3At32 = PDstandard4th3(At32, i, j, k); - PDstandard4th1At33 = PDstandard4th1(At33, i, j, k); - PDstandard4th2At33 = PDstandard4th2(At33, i, j, k); - PDstandard4th1gt11 = PDstandard4th1(gt11, i, j, k); - PDstandard4th2gt11 = PDstandard4th2(gt11, i, j, k); - PDstandard4th3gt11 = PDstandard4th3(gt11, i, j, k); - PDstandard4th11gt11 = PDstandard4th11(gt11, i, j, k); - PDstandard4th22gt11 = PDstandard4th22(gt11, i, j, k); - PDstandard4th33gt11 = PDstandard4th33(gt11, i, j, k); - PDstandard4th12gt11 = PDstandard4th12(gt11, i, j, k); - PDstandard4th13gt11 = PDstandard4th13(gt11, i, j, k); - PDstandard4th23gt11 = PDstandard4th23(gt11, i, j, k); - PDstandard4th1gt21 = PDstandard4th1(gt21, i, j, k); - PDstandard4th2gt21 = PDstandard4th2(gt21, i, j, k); - PDstandard4th3gt21 = PDstandard4th3(gt21, i, j, k); - PDstandard4th11gt21 = PDstandard4th11(gt21, i, j, k); - PDstandard4th22gt21 = PDstandard4th22(gt21, i, j, k); - PDstandard4th33gt21 = PDstandard4th33(gt21, i, j, k); - PDstandard4th12gt21 = PDstandard4th12(gt21, i, j, k); - PDstandard4th13gt21 = PDstandard4th13(gt21, i, j, k); - PDstandard4th23gt21 = PDstandard4th23(gt21, i, j, k); - PDstandard4th1gt22 = PDstandard4th1(gt22, i, j, k); - PDstandard4th2gt22 = PDstandard4th2(gt22, i, j, k); - PDstandard4th3gt22 = PDstandard4th3(gt22, i, j, k); - PDstandard4th11gt22 = PDstandard4th11(gt22, i, j, k); - PDstandard4th22gt22 = PDstandard4th22(gt22, i, j, k); - PDstandard4th33gt22 = PDstandard4th33(gt22, i, j, k); - PDstandard4th12gt22 = PDstandard4th12(gt22, i, j, k); - PDstandard4th13gt22 = PDstandard4th13(gt22, i, j, k); - PDstandard4th23gt22 = PDstandard4th23(gt22, i, j, k); - PDstandard4th1gt31 = PDstandard4th1(gt31, i, j, k); - PDstandard4th2gt31 = PDstandard4th2(gt31, i, j, k); - PDstandard4th3gt31 = PDstandard4th3(gt31, i, j, k); - PDstandard4th11gt31 = PDstandard4th11(gt31, i, j, k); - PDstandard4th22gt31 = PDstandard4th22(gt31, i, j, k); - PDstandard4th33gt31 = PDstandard4th33(gt31, i, j, k); - PDstandard4th12gt31 = PDstandard4th12(gt31, i, j, k); - PDstandard4th13gt31 = PDstandard4th13(gt31, i, j, k); - PDstandard4th23gt31 = PDstandard4th23(gt31, i, j, k); - PDstandard4th1gt32 = PDstandard4th1(gt32, i, j, k); - PDstandard4th2gt32 = PDstandard4th2(gt32, i, j, k); - PDstandard4th3gt32 = PDstandard4th3(gt32, i, j, k); - PDstandard4th11gt32 = PDstandard4th11(gt32, i, j, k); - PDstandard4th22gt32 = PDstandard4th22(gt32, i, j, k); - PDstandard4th33gt32 = PDstandard4th33(gt32, i, j, k); - PDstandard4th12gt32 = PDstandard4th12(gt32, i, j, k); - PDstandard4th13gt32 = PDstandard4th13(gt32, i, j, k); - PDstandard4th23gt32 = PDstandard4th23(gt32, i, j, k); - PDstandard4th1gt33 = PDstandard4th1(gt33, i, j, k); - PDstandard4th2gt33 = PDstandard4th2(gt33, i, j, k); - PDstandard4th3gt33 = PDstandard4th3(gt33, i, j, k); - PDstandard4th11gt33 = PDstandard4th11(gt33, i, j, k); - PDstandard4th22gt33 = PDstandard4th22(gt33, i, j, k); - PDstandard4th33gt33 = PDstandard4th33(gt33, i, j, k); - PDstandard4th12gt33 = PDstandard4th12(gt33, i, j, k); - PDstandard4th13gt33 = PDstandard4th13(gt33, i, j, k); - PDstandard4th23gt33 = PDstandard4th23(gt33, i, j, k); - PDstandard4th1phi = PDstandard4th1(phi, i, j, k); - PDstandard4th2phi = PDstandard4th2(phi, i, j, k); - PDstandard4th3phi = PDstandard4th3(phi, i, j, k); - PDstandard4th11phi = PDstandard4th11(phi, i, j, k); - PDstandard4th22phi = PDstandard4th22(phi, i, j, k); - PDstandard4th33phi = PDstandard4th33(phi, i, j, k); - PDstandard4th12phi = PDstandard4th12(phi, i, j, k); - PDstandard4th13phi = PDstandard4th13(phi, i, j, k); - PDstandard4th23phi = PDstandard4th23(phi, i, j, k); - PDstandard4th1trK = PDstandard4th1(trK, i, j, k); - PDstandard4th2trK = PDstandard4th2(trK, i, j, k); - PDstandard4th3trK = PDstandard4th3(trK, i, j, k); - PDstandard4th1Xt1 = PDstandard4th1(Xt1, i, j, k); - PDstandard4th2Xt1 = PDstandard4th2(Xt1, i, j, k); - PDstandard4th3Xt1 = PDstandard4th3(Xt1, i, j, k); - PDstandard4th1Xt2 = PDstandard4th1(Xt2, i, j, k); - PDstandard4th2Xt2 = PDstandard4th2(Xt2, i, j, k); - PDstandard4th3Xt2 = PDstandard4th3(Xt2, i, j, k); - PDstandard4th1Xt3 = PDstandard4th1(Xt3, i, j, k); - PDstandard4th2Xt3 = PDstandard4th2(Xt3, i, j, k); - PDstandard4th3Xt3 = PDstandard4th3(Xt3, i, j, k); - - /* Precompute derivatives (old style) */ - - /* Calculate temporaries and grid functions */ - detgt = 2*gt21L*gt31L*gt32L + gt33L*(gt11L*gt22L - SQR(gt21L)) - gt22L*SQR(gt31L) - gt11L*SQR(gt32L); - - ddetgt1 = (-2*gt22L*gt31L + 2*gt21L*gt32L)*PDstandard4th1gt31 + - 2*((gt31L*gt32L - gt21L*gt33L)*PDstandard4th1gt21 + (gt21L*gt31L - gt11L*gt32L)*PDstandard4th1gt32) + - PDstandard4th1gt33*(gt11L*gt22L - SQR(gt21L)) + PDstandard4th1gt22*(gt11L*gt33L - SQR(gt31L)) + - PDstandard4th1gt11*(gt22L*gt33L - SQR(gt32L)); - - ddetgt2 = (-2*gt22L*gt31L + 2*gt21L*gt32L)*PDstandard4th2gt31 + - 2*((gt31L*gt32L - gt21L*gt33L)*PDstandard4th2gt21 + (gt21L*gt31L - gt11L*gt32L)*PDstandard4th2gt32) + - PDstandard4th2gt33*(gt11L*gt22L - SQR(gt21L)) + PDstandard4th2gt22*(gt11L*gt33L - SQR(gt31L)) + - PDstandard4th2gt11*(gt22L*gt33L - SQR(gt32L)); - - ddetgt3 = (-2*gt22L*gt31L + 2*gt21L*gt32L)*PDstandard4th3gt31 + - 2*((gt31L*gt32L - gt21L*gt33L)*PDstandard4th3gt21 + (gt21L*gt31L - gt11L*gt32L)*PDstandard4th3gt32) + - PDstandard4th3gt33*(gt11L*gt22L - SQR(gt21L)) + PDstandard4th3gt22*(gt11L*gt33L - SQR(gt31L)) + - PDstandard4th3gt11*(gt22L*gt33L - SQR(gt32L)); - - gtu11 = INV(detgt)*(gt22L*gt33L - SQR(gt32L)); - - gtu12 = (gt31L*gt32L - gt21L*gt33L)*INV(detgt); - - gtu13 = (-(gt22L*gt31L) + gt21L*gt32L)*INV(detgt); - - gtu22 = INV(detgt)*(gt11L*gt33L - SQR(gt31L)); - - gtu23 = (gt21L*gt31L - gt11L*gt32L)*INV(detgt); - - gtu33 = INV(detgt)*(gt11L*gt22L - SQR(gt21L)); - - Gt111 = khalf*(gtu11*PDstandard4th1gt11 + 2*(gtu12*PDstandard4th1gt21 + gtu13*PDstandard4th1gt31) - - gtu12*PDstandard4th2gt11 - gtu13*PDstandard4th3gt11); - - Gt211 = khalf*(gtu12*PDstandard4th1gt11 + 2*(gtu22*PDstandard4th1gt21 + gtu23*PDstandard4th1gt31) - - gtu22*PDstandard4th2gt11 - gtu23*PDstandard4th3gt11); - - Gt311 = khalf*(gtu13*PDstandard4th1gt11 + 2*(gtu23*PDstandard4th1gt21 + gtu33*PDstandard4th1gt31) - - gtu23*PDstandard4th2gt11 - gtu33*PDstandard4th3gt11); - - Gt121 = khalf*(gtu12*PDstandard4th1gt22 + gtu11*PDstandard4th2gt11 + - gtu13*(PDstandard4th1gt32 + PDstandard4th2gt31 - PDstandard4th3gt21)); - - Gt221 = khalf*(gtu22*PDstandard4th1gt22 + gtu12*PDstandard4th2gt11 + - gtu23*(PDstandard4th1gt32 + PDstandard4th2gt31 - PDstandard4th3gt21)); - - Gt321 = khalf*(gtu23*PDstandard4th1gt22 + gtu13*PDstandard4th2gt11 + - gtu33*(PDstandard4th1gt32 + PDstandard4th2gt31 - PDstandard4th3gt21)); - - Gt131 = khalf*(gtu13*PDstandard4th1gt33 + gtu11*PDstandard4th3gt11 + - gtu12*(PDstandard4th1gt32 - PDstandard4th2gt31 + PDstandard4th3gt21)); - - Gt231 = khalf*(gtu23*PDstandard4th1gt33 + gtu12*PDstandard4th3gt11 + - gtu22*(PDstandard4th1gt32 - PDstandard4th2gt31 + PDstandard4th3gt21)); - - Gt331 = khalf*(gtu33*PDstandard4th1gt33 + gtu13*PDstandard4th3gt11 + - gtu23*(PDstandard4th1gt32 - PDstandard4th2gt31 + PDstandard4th3gt21)); - - Gt122 = khalf*(gtu11*(-PDstandard4th1gt22 + 2*PDstandard4th2gt21) + gtu12*PDstandard4th2gt22 + - gtu13*(2*PDstandard4th2gt32 - PDstandard4th3gt22)); - - Gt222 = khalf*(gtu12*(-PDstandard4th1gt22 + 2*PDstandard4th2gt21) + gtu22*PDstandard4th2gt22 + - gtu23*(2*PDstandard4th2gt32 - PDstandard4th3gt22)); - - Gt322 = khalf*(gtu13*(-PDstandard4th1gt22 + 2*PDstandard4th2gt21) + gtu23*PDstandard4th2gt22 + - gtu33*(2*PDstandard4th2gt32 - PDstandard4th3gt22)); - - Gt132 = khalf*(gtu13*PDstandard4th2gt33 + gtu11*(-PDstandard4th1gt32 + PDstandard4th2gt31 + PDstandard4th3gt21) + - gtu12*PDstandard4th3gt22); - - Gt232 = khalf*(gtu23*PDstandard4th2gt33 + gtu12*(-PDstandard4th1gt32 + PDstandard4th2gt31 + PDstandard4th3gt21) + - gtu22*PDstandard4th3gt22); - - Gt332 = khalf*(gtu33*PDstandard4th2gt33 + gtu13*(-PDstandard4th1gt32 + PDstandard4th2gt31 + PDstandard4th3gt21) + - gtu23*PDstandard4th3gt22); - - Gt133 = khalf*(-(gtu11*PDstandard4th1gt33) - gtu12*PDstandard4th2gt33 + 2*gtu11*PDstandard4th3gt31 + - 2*gtu12*PDstandard4th3gt32 + gtu13*PDstandard4th3gt33); - - Gt233 = khalf*(-(gtu12*PDstandard4th1gt33) - gtu22*PDstandard4th2gt33 + 2*gtu12*PDstandard4th3gt31 + - 2*gtu22*PDstandard4th3gt32 + gtu23*PDstandard4th3gt33); - - Gt333 = khalf*(-(gtu13*PDstandard4th1gt33) - gtu23*PDstandard4th2gt33 + 2*gtu13*PDstandard4th3gt31 + - 2*gtu23*PDstandard4th3gt32 + gtu33*PDstandard4th3gt33); - - e4phi = exp(4*phiL); - - g11 = e4phi*gt11L; - - g21 = e4phi*gt21L; - - g31 = e4phi*gt31L; - - g22 = e4phi*gt22L; - - g32 = e4phi*gt32L; - - g33 = e4phi*gt33L; - - detg = 2*g21*g31*g32 + g33*(g11*g22 - SQR(g21)) - g22*SQR(g31) - g11*SQR(g32); - - gu11 = INV(detg)*(g22*g33 - SQR(g32)); - - gu12 = (g31*g32 - g21*g33)*INV(detg); - - gu13 = (-(g22*g31) + g21*g32)*INV(detg); - - gu22 = INV(detg)*(g11*g33 - SQR(g31)); - - gu23 = (g21*g31 - g11*g32)*INV(detg); - - gu33 = INV(detg)*(g11*g22 - SQR(g21)); - - ddetg1 = e4phi*(ddetgt1 + 4*detgt*PDstandard4th1phi); - - ddetg2 = e4phi*(ddetgt2 + 4*detgt*PDstandard4th2phi); - - ddetg3 = e4phi*(ddetgt3 + 4*detgt*PDstandard4th3phi); - - G111 = -((-6*detg*Gt111 + ddetg1*(-6 + g11*gu11) + g11*(ddetg2*gu12 + ddetg3*gu13))*INV(detg))/6.; - - G211 = ((6*detg*Gt211 - g11*(ddetg1*gu12 + ddetg2*gu22 + ddetg3*gu23))*INV(detg))/6.; - - G311 = ((6*detg*Gt311 - g11*(ddetg1*gu13 + ddetg2*gu23 + ddetg3*gu33))*INV(detg))/6.; - - G121 = -((-6*detg*Gt121 + ddetg2*(-3 + g21*gu12) + g21*(ddetg1*gu11 + ddetg3*gu13))*INV(detg))/6.; - - G221 = -((-6*detg*Gt221 + ddetg1*(-3 + g21*gu12) + g21*(ddetg2*gu22 + ddetg3*gu23))*INV(detg))/6.; - - G321 = ((6*detg*Gt321 - g21*(ddetg1*gu13 + ddetg2*gu23 + ddetg3*gu33))*INV(detg))/6.; - - G131 = -((-6*detg*Gt131 + g31*(ddetg1*gu11 + ddetg2*gu12) + ddetg3*(-3 + g31*gu13))*INV(detg))/6.; - - G231 = ((6*detg*Gt231 - g31*(ddetg1*gu12 + ddetg2*gu22 + ddetg3*gu23))*INV(detg))/6.; - - G331 = -((-6*detg*Gt331 + ddetg1*(-3 + g31*gu13) + g31*(ddetg2*gu23 + ddetg3*gu33))*INV(detg))/6.; - - G122 = ((6*detg*Gt122 - g22*(ddetg1*gu11 + ddetg2*gu12 + ddetg3*gu13))*INV(detg))/6.; - - G222 = -((-6*detg*Gt222 + ddetg2*(-6 + g22*gu22) + g22*(ddetg1*gu12 + ddetg3*gu23))*INV(detg))/6.; - - G322 = ((6*detg*Gt322 - g22*(ddetg1*gu13 + ddetg2*gu23 + ddetg3*gu33))*INV(detg))/6.; - - G132 = ((6*detg*Gt132 - g32*(ddetg1*gu11 + ddetg2*gu12 + ddetg3*gu13))*INV(detg))/6.; - - G232 = -((-6*detg*Gt232 + g32*(ddetg1*gu12 + ddetg2*gu22) + ddetg3*(-3 + g32*gu23))*INV(detg))/6.; - - G332 = -((-6*detg*Gt332 + ddetg2*(-3 + g32*gu23) + g32*(ddetg1*gu13 + ddetg3*gu33))*INV(detg))/6.; - - G133 = ((6*detg*Gt133 - g33*(ddetg1*gu11 + ddetg2*gu12 + ddetg3*gu13))*INV(detg))/6.; - - G233 = ((6*detg*Gt233 - g33*(ddetg1*gu12 + ddetg2*gu22 + ddetg3*gu23))*INV(detg))/6.; - - G333 = -((-6*detg*Gt333 + g33*(ddetg1*gu13 + ddetg2*gu23) + ddetg3*(-6 + g33*gu33))*INV(detg))/6.; - - Rt11 = -(gtu11*khalf*PDstandard4th11gt11) + gtu12* - (2*Gt211*Gt221*gt22L + 4*gt21L*Gt232*Gt311 + 6*Gt121*Gt311*gt31L + 4*gt11L*Gt131*Gt321 + 4*Gt221*gt31L*Gt321 + - 4*Gt211*gt31L*Gt322 + 2*Gt221*Gt311*gt32L + 2*Gt211*Gt321*gt32L + 4*gt31L*Gt321*Gt331 + 4*Gt311*gt31L*Gt332 + - 2*Gt311*Gt321*gt33L - PDstandard4th12gt11) + - gtu13*(2*Gt231*Gt311*gt32L + 4*gt11L*Gt131*Gt331 + 2*Gt211*gt32L*Gt331 + 4*Gt211*gt31L*Gt332 + - 4*Gt311*gt31L*Gt333 + 2*Gt311*Gt331*gt33L - PDstandard4th13gt11) + 2*gt11L*PDstandard4th1Xt1 + - gt21L*(6*Gt111*Gt221*gtu12 + 4*Gt211*Gt222*gtu12 + 6*Gt131*Gt211*gtu13 + 6*Gt121*Gt221*gtu22 + - 6*Gt131*Gt221*gtu23 + 6*Gt131*Gt231*gtu33 + 2*PDstandard4th1Xt2) + 2*gt31L*PDstandard4th1Xt3 - - gtu22*khalf*PDstandard4th22gt11 - gtu23*PDstandard4th23gt11 - gtu33*khalf*PDstandard4th33gt11 + - Gt111*(10*gt11L*Gt131*gtu13 + 6*gt21L*Gt231*gtu13 + 2*gt11L*Xt1L) + - Gt211*(4*gt11L*Gt121*gtu11 + 6*Gt111*gt21L*gtu11 + 4*gt31L*Gt321*gtu11 + 2*Gt311*gt32L*gtu11 + - 4*gt11L*Gt122*gtu12 + 4*gt11L*Gt132*gtu13 + 2*gt22L*Gt231*gtu13 + 4*gt21L*Gt232*gtu13 + 2*gt21L*Xt1L) + - Gt311*(4*gt21L*Gt231*gtu11 + 6*Gt111*gt31L*gtu11 + 4*gt11L*Gt132*gtu12 + 4*gt11L*Gt133*gtu13 + - 4*gt21L*Gt233*gtu13 + 6*Gt131*gt31L*gtu13 + 2*gt31L*Xt1L) + - Gt121*(10*Gt111*gt11L*gtu12 + 6*Gt211*gt21L*gtu12 + 4*gt11L*Gt231*gtu13 + 6*gt31L*Gt321*gtu22 + - 10*gt11L*Gt131*gtu23 + 2*gt11L*Xt2L) + Gt221* - (4*gt31L*Gt322*gtu22 + 2*Gt321*gt32L*gtu22 + 2*gt22L*Gt231*gtu23 + 2*gt32L*Gt331*gtu23 + 4*gt31L*Gt332*gtu23 + - 2*gt21L*Xt2L) + Gt321*(4*gt21L*Gt231*gtu12 + 6*Gt111*gt31L*gtu12 + 4*Gt231*gt31L*gtu13 + 4*gt11L*Gt132*gtu22 + - 4*gt31L*Gt332*gtu22 + 4*gt11L*Gt133*gtu23 + 4*gt21L*Gt233*gtu23 + 6*Gt131*gt31L*gtu23 + 2*gt31L*Xt2L) + - 2*gt11L*Gt131*Xt3L + Gt231*(4*gt11L*Gt122*gtu23 + 6*Gt121*gt21L*gtu23 + 4*gt31L*Gt322*gtu23 + 2*Gt321*gt32L*gtu23 + - 4*gt11L*Gt132*gtu33 + 2*gt32L*Gt331*gtu33 + 4*gt31L*Gt332*gtu33 + 2*gt21L*Xt3L) + - Gt331*(4*gt21L*Gt231*gtu13 + 6*Gt111*gt31L*gtu13 + 4*gt11L*Gt132*gtu23 + 4*gt21L*Gt232*gtu23 + - 6*Gt121*gt31L*gtu23 + 2*Gt321*gt33L*gtu23 + 4*gt11L*Gt133*gtu33 + 4*gt21L*Gt233*gtu33 + 6*Gt131*gt31L*gtu33 + - 4*gt31L*Gt333*gtu33 + 2*gt31L*Xt3L) + 5*gt11L*gtu11*SQR(Gt111) + 5*gt11L*gtu22*SQR(Gt121) + - 5*gt11L*gtu33*SQR(Gt131) + gt22L*gtu11*SQR(Gt211) + gt22L*gtu22*SQR(Gt221) + - 4*(Gt211*gt21L*Gt221*gtu11 + gt11L*Gt131*Gt311*gtu11 + Gt311*gt31L*Gt331*gtu11 + gt11L*Gt121*Gt221*gtu12 + - gt21L*Gt221*Gt231*gtu13 + gt11L*Gt122*Gt221*gtu22 + gt21L*Gt221*Gt222*gtu22 + gt21L*Gt232*Gt321*gtu22 + - gt11L*Gt132*Gt221*gtu23 + gt21L*Gt222*Gt231*gtu23 + gt21L*Gt221*Gt232*gtu23 + gt31L*Gt331*Gt332*gtu23 + - gt31L*Gt321*Gt333*gtu23 + gt21L*Gt231*Gt232*gtu33 + gt21L*gtu12*SQR(Gt221)) + gt22L*gtu33*SQR(Gt231) + - gt33L*gtu11*SQR(Gt311) + gt33L*gtu22*SQR(Gt321) + 4*gt31L*gtu13*SQR(Gt331) + gt33L*gtu33*SQR(Gt331); - - Rt21 = gt22L*PDstandard4th1Xt2 + gt32L*PDstandard4th1Xt3 + gt11L*PDstandard4th2Xt1 + - gt21L*(PDstandard4th1Xt1 + PDstandard4th2Xt2) + gt31L*PDstandard4th2Xt3 + - khalf*(-(gtu11*PDstandard4th11gt21) - 2*gtu12*PDstandard4th12gt21 - 2*gtu13*PDstandard4th13gt21 - - gtu22*PDstandard4th22gt21 - 2*gtu23*PDstandard4th23gt21 - gtu33*PDstandard4th33gt21) + - (gt11L*Gt121 + gt21L*(Gt111 + Gt221) + Gt211*gt22L + gt31L*Gt321 + Gt311*gt32L)*Xt1L + - (gt11L*Gt122 + gt21L*(Gt121 + Gt222) + Gt221*gt22L + gt31L*Gt322 + Gt321*gt32L)*Xt2L + - (gt11L*Gt132 + gt22L*Gt231 + gt21L*(Gt131 + Gt232) + gt32L*Gt331 + gt31L*Gt332)*Xt3L + - gtu11*(gt11L*(Gt121*(3*Gt111 + 2*Gt221) + 2*Gt131*Gt321) + - Gt111*(gt21L*Gt221 + gt31L*Gt321 + 2*(Gt211*gt22L + Gt311*gt32L)) + - Gt211*(5*Gt121*gt21L + 3*(Gt221*gt22L + Gt321*gt32L)) + gt31L*(3*Gt121*Gt311 + 2*Gt321*(Gt221 + Gt331)) + - Gt311*(Gt221*gt32L + Gt321*gt33L) + 2*(Gt311*(Gt131*gt21L + gt22L*Gt231 + gt32L*Gt331) + - gt21L*(Gt231*Gt321 + SQR(Gt111) + SQR(Gt221)))) + - gtu22*(gt11L*(Gt122*(3*Gt121 + 2*Gt222) + 2*Gt132*Gt322) + - Gt121*(gt21L*Gt222 + gt31L*Gt322 + 2*(Gt221*gt22L + Gt321*gt32L)) + - Gt221*(5*Gt122*gt21L + 3*(Gt222*gt22L + Gt322*gt32L)) + gt31L*(3*Gt122*Gt321 + 2*Gt322*(Gt222 + Gt332)) + - Gt321*(Gt222*gt32L + Gt322*gt33L) + 2*(Gt321*(Gt132*gt21L + gt22L*Gt232 + gt32L*Gt332) + - gt21L*(Gt232*Gt322 + SQR(Gt121) + SQR(Gt222)))) + - gtu33*(gt11L*(Gt132*(3*Gt131 + 2*Gt232) + 2*Gt133*Gt332) + - Gt131*(gt21L*Gt232 + 2*(gt22L*Gt231 + gt32L*Gt331) + gt31L*Gt332) + - Gt231*(5*Gt132*gt21L + 3*(gt22L*Gt232 + gt32L*Gt332)) + - Gt331*(3*Gt132*gt31L + gt32L*(Gt232 + 2*Gt333) + Gt332*gt33L) + - 2*((Gt133*gt21L + gt22L*Gt233)*Gt331 + Gt332*(gt21L*Gt233 + gt31L*(Gt232 + Gt333)) + - gt21L*(SQR(Gt131) + SQR(Gt232)))) + gtu12* - (Gt122*(Gt111*gt11L + 3*Gt211*gt21L + Gt311*gt31L) + - gt21L*(Gt111*Gt222 + Gt221*(4*Gt121 + 2*Gt222) + 2*(Gt132*Gt311 + Gt232*Gt321)) + - gt31L*((Gt111 + 2*Gt221)*Gt322 + Gt321*(4*Gt121 + 2*Gt332)) + - gt32L*(2*Gt221*Gt321 + 3*Gt211*Gt322 + Gt311*(Gt222 + 2*Gt332)) + - gt11L*(2*(Gt122*Gt221 + Gt132*Gt321) + 3*SQR(Gt121)) + gt22L*(3*Gt211*Gt222 + 2*Gt232*Gt311 + SQR(Gt221)) + - 2*(gt11L*(Gt111*Gt122 + Gt121*Gt222 + Gt131*Gt322) + - gt21L*(Gt122*Gt211 + Gt121*(Gt111 + Gt221) + Gt221*Gt222 + Gt131*Gt321 + Gt231*Gt322) + - Gt111*(Gt121*gt21L + Gt221*gt22L + Gt321*gt32L) + gt31L*(Gt122*Gt311 + Gt222*Gt321 + Gt322*Gt331) + - gt32L*(Gt121*Gt311 + Gt321*(Gt221 + Gt331)) + gt22L*(Gt121*Gt211 + Gt231*Gt321 + SQR(Gt221))) + - gt33L*(Gt311*Gt322 + SQR(Gt321))) + gtu13* - (gt22L*(3*Gt211*Gt232 + 2*Gt233*Gt311) + Gt132*(Gt111*gt11L + 3*Gt211*gt21L + Gt311*gt31L) + - gt11L*(3*Gt121*Gt131 + 2*(Gt132*Gt221 + Gt133*Gt321)) + - gt21L*(Gt131*Gt221 + 3*Gt121*Gt231 + Gt111*Gt232 + 2*(Gt133*Gt311 + Gt233*Gt321)) + - Gt221*(gt22L*Gt231 + 2*gt21L*Gt232 + gt32L*Gt331) + gt31L*(3*Gt121*Gt331 + Gt111*Gt332) + - gt32L*(3*Gt211*Gt332 + Gt311*(Gt232 + 2*Gt333)) + Gt311*Gt332*gt33L + - Gt321*(Gt231*gt32L + gt31L*(Gt131 + 2*Gt333) + Gt331*gt33L) + - 2*(Gt111*(Gt131*gt21L + gt22L*Gt231 + gt32L*Gt331) + gt22L*(Gt131*Gt211 + Gt231*(Gt221 + Gt331)) + - gt11L*(Gt111*Gt132 + Gt121*Gt232 + Gt131*Gt332) + gt31L*(Gt132*Gt311 + Gt232*Gt321 + (Gt221 + Gt331)*Gt332) + - gt21L*(Gt132*Gt211 + Gt221*Gt232 + Gt131*(Gt111 + Gt331) + Gt231*(Gt121 + Gt332)) + - gt32L*(Gt131*Gt311 + Gt231*Gt321 + SQR(Gt331)))) + - gtu23*(gt22L*(3*Gt221*Gt232 + 2*Gt233*Gt321) + Gt132*(gt11L*Gt121 + 3*gt21L*Gt221 + gt31L*Gt321) + - gt11L*(3*Gt122*Gt131 + 2*(Gt132*Gt222 + Gt133*Gt322)) + - gt21L*(Gt131*Gt222 + 3*Gt122*Gt231 + Gt121*Gt232 + 2*(Gt133*Gt321 + Gt233*Gt322)) + - Gt222*(gt22L*Gt231 + 2*gt21L*Gt232 + gt32L*Gt331) + gt31L*(3*Gt122*Gt331 + Gt121*Gt332) + - gt32L*(3*Gt221*Gt332 + Gt321*(Gt232 + 2*Gt333)) + Gt321*Gt332*gt33L + - Gt322*(Gt231*gt32L + gt31L*(Gt131 + 2*Gt333) + Gt331*gt33L) + - 2*(Gt131*(Gt121*gt21L + Gt221*gt22L + Gt321*gt32L) + Gt231*(Gt122*gt21L + Gt222*gt22L + Gt322*gt32L) + - Gt121*(Gt131*gt21L + gt22L*Gt231 + gt32L*Gt331) + Gt331*(Gt132*gt21L + gt22L*Gt232 + gt32L*Gt332) + - gt11L*(Gt122*Gt232 + Gt132*(Gt121 + Gt332)) + gt21L*(Gt132*Gt221 + Gt232*(Gt222 + Gt332)) + - gt31L*(Gt132*Gt321 + Gt232*Gt322 + Gt222*Gt332 + SQR(Gt332)))); - - Rt31 = gt32L*PDstandard4th1Xt2 + gt33L*PDstandard4th1Xt3 + - khalf*(-(gtu11*PDstandard4th11gt31) - 2*gtu12*PDstandard4th12gt31 - 2*gtu13*PDstandard4th13gt31 - - gtu22*PDstandard4th22gt31 - 2*gtu23*PDstandard4th23gt31 - gtu33*PDstandard4th33gt31) + gt11L*PDstandard4th3Xt1 + - gt21L*PDstandard4th3Xt2 + gt31L*(PDstandard4th1Xt1 + PDstandard4th3Xt3) + - (gt11L*Gt131 + gt21L*Gt231 + Gt211*gt32L + gt31L*(Gt111 + Gt331) + Gt311*gt33L)*Xt1L + - (gt11L*Gt132 + gt21L*Gt232 + Gt221*gt32L + gt31L*(Gt121 + Gt332) + Gt321*gt33L)*Xt2L + - (gt11L*Gt133 + gt21L*Gt233 + Gt231*gt32L + gt31L*(Gt131 + Gt333) + Gt331*gt33L)*Xt3L + - gtu12*(gt22L*(Gt221*Gt231 + Gt211*Gt232) + Gt132*(Gt111*gt11L + Gt211*gt21L + 3*Gt311*gt31L) + - gt21L*((Gt121 + 2*Gt222)*Gt231 + Gt232*(Gt111 + 2*Gt331)) + - gt11L*(3*Gt121*Gt131 + 2*(Gt122*Gt231 + Gt132*Gt331)) + gt31L*(Gt121*Gt331 + Gt111*Gt332) + - gt32L*(Gt231*Gt321 + Gt221*Gt331 + Gt211*(2*Gt222 + Gt332)) + (2*Gt211*Gt322 + Gt321*Gt331)*gt33L + - 3*(Gt131*(gt21L*Gt221 + gt31L*Gt321) + Gt311*(Gt232*gt32L + Gt332*gt33L)) + - 2*(gt11L*(Gt111*Gt132 + Gt121*Gt232 + Gt131*Gt332) + gt21L*(Gt132*Gt211 + Gt221*Gt232 + Gt231*Gt332) + - gt31L*(Gt122*Gt211 + Gt121*(2*Gt111 + Gt221) + Gt132*Gt311 + (Gt131 + Gt232)*Gt321 + Gt231*Gt322 + - 2*Gt331*Gt332) + (Gt121*Gt311 + Gt321*(Gt221 + Gt331))*gt33L + Gt111*(Gt221*gt32L + Gt321*gt33L) + - gt32L*(Gt121*Gt211 + Gt231*Gt321 + SQR(Gt221)))) + - gtu23*(Gt233*(Gt221*gt22L + 2*gt31L*Gt322 + Gt321*gt32L) + Gt132*(gt11L*Gt131 + gt21L*Gt231 + 3*gt31L*Gt331) + - Gt232*(Gt131*gt21L + gt22L*Gt231 + 3*gt32L*Gt331) + (Gt131*gt31L + Gt231*gt32L)*Gt332 + - gt21L*((Gt121 + 2*Gt222)*Gt233 + 2*Gt232*Gt333) + gt11L*(3*Gt121*Gt133 + 2*(Gt122*Gt233 + Gt132*Gt333)) + - Gt333*(Gt121*gt31L + Gt221*gt32L + Gt321*gt33L) + 3*(Gt133*(gt21L*Gt221 + gt31L*Gt321) + Gt331*Gt332*gt33L) + - 2*((Gt222*Gt231 + Gt221*(Gt131 + Gt232) + Gt233*Gt321)*gt32L + gt11L*(Gt132*(Gt131 + Gt232) + Gt133*Gt332) + - gt31L*(Gt121*Gt131 + Gt122*Gt231 + Gt133*Gt321 + Gt132*(Gt221 + Gt331) + Gt332*(Gt232 + 2*Gt333)) + - (Gt231*Gt322 + Gt221*Gt332 + Gt321*(Gt131 + Gt333))*gt33L + Gt121*(Gt131*gt31L + Gt231*gt32L + Gt331*gt33L) + - gt21L*(Gt132*Gt231 + Gt233*Gt332 + SQR(Gt232)))) + - gtu11*(Gt231*(Gt211*gt22L + 2*(gt21L*Gt221 + gt31L*Gt321) + 3*Gt311*gt32L) + - gt21L*(3*Gt131*Gt211 + Gt231*(Gt111 + 2*Gt331)) + gt11L*(2*Gt121*Gt231 + Gt131*(3*Gt111 + 2*Gt331)) + - 3*Gt311*Gt331*gt33L + Gt211*(gt32L*Gt331 + 2*Gt321*gt33L) + - gt31L*(5*Gt131*Gt311 + Gt111*Gt331 + 2*(Gt121*Gt211 + SQR(Gt111))) + - 2*(Gt211*(Gt111 + Gt221)*gt32L + Gt111*Gt311*gt33L + gt31L*SQR(Gt331))) + - gtu13*(Gt133*(Gt111*gt11L + Gt211*gt21L + 3*Gt311*gt31L) + gt32L*(3*Gt233*Gt311 + 2*Gt231*(Gt111 + Gt331)) + - gt21L*(Gt111*Gt233 + 2*(Gt231*Gt232 + Gt233*Gt331)) + Gt111*(gt31L*Gt333 + 2*Gt331*gt33L) + - Gt211*(gt22L*Gt233 + gt32L*Gt333 + 2*(Gt132*gt31L + Gt232*gt32L + Gt332*gt33L)) + - gt11L*(2*(Gt132*Gt231 + Gt133*Gt331) + 3*SQR(Gt131)) + gt22L*SQR(Gt231) + gt33L*(3*Gt311*Gt333 + SQR(Gt331)) + - 2*(gt32L*(Gt131*Gt211 + Gt231*(Gt221 + Gt331)) + Gt131*(Gt111*gt31L + 2*(gt21L*Gt231 + gt31L*Gt331)) + - gt11L*(Gt111*Gt133 + Gt121*Gt233 + Gt131*Gt333) + gt21L*(Gt133*Gt211 + Gt221*Gt233 + Gt231*Gt333) + - gt31L*(Gt133*Gt311 + Gt233*Gt321 + Gt131*(Gt111 + Gt331) + Gt231*(Gt121 + Gt332) + 2*Gt331*Gt333) + - gt33L*(Gt131*Gt311 + Gt231*Gt321 + SQR(Gt331)))) + - gtu22*(Gt232*(Gt221*gt22L + 2*(gt21L*Gt222 + gt31L*Gt322) + 3*Gt321*gt32L) + - gt21L*(3*Gt132*Gt221 + Gt232*(Gt121 + 2*Gt332)) + gt11L*(2*Gt122*Gt232 + Gt132*(3*Gt121 + 2*Gt332)) + - 3*Gt321*Gt332*gt33L + Gt221*(gt32L*Gt332 + 2*Gt322*gt33L) + - gt31L*(5*Gt132*Gt321 + Gt121*Gt332 + 2*(Gt122*Gt221 + SQR(Gt121))) + - 2*(Gt221*(Gt121 + Gt222)*gt32L + Gt121*Gt321*gt33L + gt31L*SQR(Gt332))) + - gtu33*(Gt233*(gt22L*Gt231 + 2*gt21L*Gt232 + 3*gt32L*Gt331) + gt21L*(3*Gt133*Gt231 + Gt233*(Gt131 + 2*Gt333)) + - gt11L*(2*Gt132*Gt233 + Gt133*(3*Gt131 + 2*Gt333)) + 3*Gt331*Gt333*gt33L + Gt231*(gt32L*Gt333 + 2*Gt332*gt33L) + - gt31L*(5*Gt133*Gt331 + Gt131*Gt333 + 2*SQR(Gt131)) + - 2*(Gt231*(Gt132*gt31L + (Gt131 + Gt232)*gt32L) + Gt131*Gt331*gt33L + gt31L*(Gt233*Gt332 + SQR(Gt333)))); - - Rt22 = 6*(Gt122*gt21L*Gt221*gtu12 + Gt121*gt21L*Gt222*gtu12 + Gt222*Gt321*gt32L*gtu12 + Gt221*Gt322*gt32L*gtu12 + - Gt132*gt21L*Gt221*gtu13 + Gt122*gt21L*Gt222*gtu22 + Gt132*gt21L*Gt222*gtu23 + Gt232*gt32L*Gt332*gtu33) - - gtu11*khalf*PDstandard4th11gt22 + gtu12*(2*Gt122*gt31L*Gt321 + 4*Gt131*gt21L*Gt322 + 4*Gt121*Gt321*gt32L + - 4*Gt322*gt32L*Gt331 + 4*Gt321*gt32L*Gt332 + 2*Gt321*Gt322*gt33L - PDstandard4th12gt22) - - gtu13*PDstandard4th13gt22 - gtu22*khalf*PDstandard4th22gt22 + - gtu23*(2*Gt132*gt31L*Gt322 + 4*Gt122*gt32L*Gt331 + 4*gt22L*Gt232*Gt332 + 4*Gt322*gt32L*Gt333 + - 2*Gt322*Gt332*gt33L - PDstandard4th23gt22) + 2*gt21L*PDstandard4th2Xt1 + - gt22L*(10*Gt221*Gt222*gtu12 + 4*Gt232*Gt321*gtu12 + 4*Gt132*Gt211*gtu13 + 4*Gt122*Gt221*gtu22 + - 4*Gt132*Gt221*gtu23 + 10*Gt222*Gt232*gtu23 + 4*Gt233*Gt322*gtu23 + 4*Gt132*Gt231*gtu33 + 2*PDstandard4th2Xt2) + - gt32L*(4*Gt121*Gt311*gtu11 + 6*Gt221*Gt321*gtu11 + 4*Gt122*Gt311*gtu12 + 4*Gt132*Gt311*gtu13 + - 6*Gt232*Gt321*gtu13 + 4*Gt121*Gt331*gtu13 + 4*Gt122*Gt321*gtu22 + 6*Gt232*Gt322*gtu23 + 4*Gt132*Gt331*gtu33 + - 2*PDstandard4th2Xt3) - gtu33*khalf*PDstandard4th33gt22 + 2*Gt221*gt22L*Xt1L + - Gt121*(4*Gt111*gt21L*gtu11 + 6*gt21L*Gt221*gtu11 + 2*gt11L*Gt122*gtu12 + 2*gt11L*Gt132*gtu13 + - 4*Gt122*gt21L*gtu22 + 4*Gt132*gt21L*gtu23 + 2*gt21L*Xt1L) + - Gt321*(4*gt22L*Gt231*gtu11 + 2*Gt121*gt31L*gtu11 + 4*gt32L*Gt331*gtu11 + 4*Gt132*gt21L*gtu12 + - 4*Gt133*gt21L*gtu13 + 4*gt22L*Gt233*gtu13 + 2*Gt132*gt31L*gtu13 + 2*Gt332*gt33L*gtu13 + 4*Gt132*gt32L*gtu23 + - 2*gt32L*Xt1L) + 2*Gt222*gt22L*Xt2L + Gt122*(2*gt11L*Gt132*gtu23 + 4*gt22L*Gt231*gtu23 + 2*gt21L*Xt2L) + - Gt322*(4*gt22L*Gt231*gtu12 + 2*Gt121*gt31L*gtu12 + 4*Gt132*gt21L*gtu22 + 2*Gt122*gt31L*gtu22 + - 6*Gt222*gt32L*gtu22 + 4*Gt133*gt21L*gtu23 + 2*gt32L*Xt2L) + 2*Gt132*gt21L*Xt3L + - Gt232*(6*Gt121*gt21L*gtu13 + 10*Gt221*gt22L*gtu13 + 4*gt22L*Gt322*gtu22 + 6*Gt122*gt21L*gtu23 + - 6*Gt132*gt21L*gtu33 + 2*gt22L*Xt3L) + Gt332* - (2*Gt121*gt31L*gtu13 + 6*Gt221*gt32L*gtu13 + 4*Gt322*gt32L*gtu22 + 4*Gt132*gt21L*gtu23 + 2*Gt122*gt31L*gtu23 + - 6*Gt222*gt32L*gtu23 + 4*Gt133*gt21L*gtu33 + 4*gt22L*Gt233*gtu33 + 2*Gt132*gt31L*gtu33 + 2*gt32L*Xt3L) + - gt11L*gtu11*SQR(Gt121) + 4*(Gt121*Gt211*gt22L*gtu11 + Gt131*gt21L*Gt321*gtu11 + Gt111*Gt122*gt21L*gtu12 + - Gt122*Gt211*gt22L*gtu12 + Gt121*Gt221*gt22L*gtu12 + Gt121*Gt131*gt21L*gtu13 + Gt111*Gt132*gt21L*gtu13 + - Gt121*gt22L*Gt231*gtu13 + Gt131*gt21L*Gt332*gtu13 + gt22L*Gt231*Gt332*gtu13 + gt32L*Gt331*Gt332*gtu13 + - Gt321*gt32L*Gt333*gtu13 + Gt122*Gt131*gt21L*gtu23 + Gt131*Gt132*gt21L*gtu33 + gt32L*Gt332*Gt333*gtu33 + - gt21L*gtu12*SQR(Gt121)) + gt11L*gtu22*SQR(Gt122) + gt11L*gtu33*SQR(Gt132) + 5*gt22L*gtu11*SQR(Gt221) + - 5*gt22L*gtu22*SQR(Gt222) + 5*gt22L*gtu33*SQR(Gt232) + gt33L*gtu11*SQR(Gt321) + gt33L*gtu22*SQR(Gt322) + - 4*gt32L*gtu23*SQR(Gt332) + gt33L*gtu33*SQR(Gt332); - - Rt32 = gt31L*PDstandard4th2Xt1 + gt33L*PDstandard4th2Xt3 + - khalf*(-(gtu11*PDstandard4th11gt32) - 2*gtu12*PDstandard4th12gt32 - 2*gtu13*PDstandard4th13gt32 - - gtu22*PDstandard4th22gt32 - 2*gtu23*PDstandard4th23gt32 - gtu33*PDstandard4th33gt32) + gt21L*PDstandard4th3Xt1 + - gt22L*PDstandard4th3Xt2 + gt32L*(PDstandard4th2Xt2 + PDstandard4th3Xt3) + - (Gt131*gt21L + gt22L*Gt231 + Gt121*gt31L + gt32L*(Gt221 + Gt331) + Gt321*gt33L)*Xt1L + - (Gt132*gt21L + gt22L*Gt232 + Gt122*gt31L + gt32L*(Gt222 + Gt332) + Gt322*gt33L)*Xt2L + - (Gt133*gt21L + gt22L*Gt233 + Gt132*gt31L + gt32L*(Gt232 + Gt333) + Gt332*gt33L)*Xt3L + - gtu12*(gt11L*(Gt122*Gt131 + Gt121*Gt132) + gt32L* - (2*(Gt132*Gt311 + Gt131*Gt321) + 5*(Gt232*Gt321 + Gt231*Gt322) + Gt222*(4*Gt221 + Gt331) + - (Gt221 + 4*Gt331)*Gt332) + gt21L*(3*Gt122*Gt231 + Gt132*(Gt221 + 2*Gt331) + Gt131*(Gt222 + 2*Gt332)) + - 3*(Gt121*gt21L*Gt232 + gt22L*(Gt222*Gt231 + Gt221*Gt232) + gt31L*(Gt132*Gt321 + Gt131*Gt322) + - (Gt322*Gt331 + Gt321*Gt332)*gt33L) + 2* - ((Gt122*Gt221 + Gt121*Gt222)*gt31L + Gt111*(Gt132*gt21L + Gt122*gt31L) + - gt22L*(Gt132*Gt211 + Gt131*Gt221 + Gt232*Gt331 + Gt231*Gt332) + (Gt222*Gt321 + Gt221*Gt322)*gt33L + - Gt122*(Gt211*gt32L + Gt311*gt33L) + Gt121*(Gt131*gt21L + Gt221*gt32L + Gt321*gt33L)) + - gt31L*(Gt122*Gt331 + Gt121*Gt332 + 2*SQR(Gt121))) + - gtu13*(gt11L*(Gt131*Gt132 + Gt121*Gt133) + Gt133*(gt21L*Gt221 + 3*gt31L*Gt321 + 2*(Gt211*gt22L + Gt311*gt32L)) + - gt32L*(5*Gt233*Gt321 + Gt232*(4*Gt221 + Gt331) + Gt231*(2*Gt121 + 5*Gt332) + Gt221*Gt333) + - Gt331*(2*gt22L*Gt233 + Gt132*gt31L + 4*gt32L*Gt333) + - Gt131*(2*(gt22L*Gt231 + Gt121*gt31L + gt32L*Gt331) + gt31L*Gt332 + gt21L*(Gt232 + 2*Gt333)) + - Gt121*(gt31L*Gt333 + 2*Gt331*gt33L) + 3*(Gt231*(Gt132*gt21L + gt22L*Gt232) + (Gt121*gt21L + Gt221*gt22L)*Gt233 + - (Gt331*Gt332 + Gt321*Gt333)*gt33L) + 2* - (Gt111*(Gt133*gt21L + Gt132*gt31L) + gt31L*(Gt132*Gt221 + Gt121*Gt232 + Gt131*Gt332) + gt22L*Gt231*Gt333 + - (Gt232*Gt321 + Gt221*Gt332)*gt33L + Gt132*(Gt211*gt32L + Gt311*gt33L) + gt21L*(Gt133*Gt331 + SQR(Gt131)))) + - gtu11*(Gt131*(gt11L*Gt121 + gt21L*Gt221 + 3*gt31L*Gt321 + 2*(Gt211*gt22L + Gt311*gt32L)) + - (Gt121*gt31L + Gt221*gt32L)*Gt331 + Gt231*(5*Gt321*gt32L + 2*gt22L*Gt331) + - 3*((Gt121*gt21L + Gt221*gt22L)*Gt231 + Gt321*Gt331*gt33L) + - 2*(Gt111*(Gt131*gt21L + Gt121*gt31L) + Gt131*gt21L*Gt331 + Gt221*Gt321*gt33L + - Gt121*(Gt221*gt31L + Gt211*gt32L + Gt311*gt33L) + gt32L*(SQR(Gt221) + SQR(Gt331)))) + - gtu23*(Gt233*(Gt122*gt21L + Gt222*gt22L + 3*Gt322*gt32L) + - Gt133*(gt11L*Gt122 + 3*gt31L*Gt322 + gt21L*(Gt222 + 2*Gt332)) + - Gt132*(4*gt21L*Gt232 + 2*(gt22L*Gt231 + gt32L*Gt331 + gt31L*(Gt222 + Gt332))) + - gt32L*(4*Gt232*Gt332 + Gt222*Gt333) + Gt122*(gt31L*Gt333 + 2*Gt331*gt33L) + gt11L*SQR(Gt132) + - 3*(Gt322*Gt333*gt33L + gt22L*SQR(Gt232)) + gt33L*SQR(Gt332) + - 2*(Gt131*(Gt132*gt21L + Gt122*gt31L) + Gt133*(Gt121*gt21L + Gt221*gt22L + Gt321*gt32L) + - gt31L*(Gt122*Gt232 + Gt132*(Gt121 + Gt332)) + Gt233*(Gt122*gt21L + Gt322*gt32L + gt22L*(Gt222 + Gt332)) + - (Gt132*gt21L + gt22L*Gt232)*Gt333 + gt32L* - (Gt132*Gt221 + Gt122*Gt231 + Gt232*(2*Gt222 + Gt332) + 2*Gt332*Gt333) + - gt33L*(Gt132*Gt321 + Gt232*Gt322 + Gt222*Gt332 + SQR(Gt332)))) + - gtu22*(Gt132*(gt11L*Gt122 + gt21L*Gt222 + 3*gt31L*Gt322 + 2*(Gt221*gt22L + Gt321*gt32L)) + - (Gt122*gt31L + Gt222*gt32L)*Gt332 + Gt232*(5*Gt322*gt32L + 2*gt22L*Gt332) + - 3*((Gt122*gt21L + Gt222*gt22L)*Gt232 + Gt322*Gt332*gt33L) + - 2*(Gt121*(Gt132*gt21L + Gt122*gt31L) + Gt132*gt21L*Gt332 + Gt222*Gt322*gt33L + - Gt122*(Gt222*gt31L + Gt221*gt32L + Gt321*gt33L) + gt32L*(SQR(Gt222) + SQR(Gt332)))) + - gtu33*(Gt133*(gt11L*Gt132 + gt21L*Gt232 + 2*(gt22L*Gt231 + gt32L*Gt331)) + (Gt132*gt31L + Gt232*gt32L)*Gt333 + - Gt233*(5*gt32L*Gt332 + 2*gt22L*Gt333) + 3* - ((Gt132*gt21L + gt22L*Gt232)*Gt233 + Gt332*(Gt133*gt31L + Gt333*gt33L)) + - 2*(Gt131*(Gt133*gt21L + Gt132*gt31L) + Gt133*gt21L*Gt333 + Gt232*Gt332*gt33L + - Gt132*(Gt232*gt31L + Gt231*gt32L + Gt331*gt33L) + gt32L*(SQR(Gt232) + SQR(Gt333)))); - - Rt33 = 6*(Gt133*gt31L*Gt331*gtu13 + Gt233*gt32L*Gt331*gtu13 + Gt131*gt31L*Gt333*gtu13 + Gt231*gt32L*Gt333*gtu13 + - Gt132*gt31L*Gt332*gtu22 + Gt133*gt31L*Gt332*gtu23 + Gt132*gt31L*Gt333*gtu23 + Gt232*gt32L*Gt333*gtu23 + - Gt133*gt31L*Gt333*gtu33) + gtu12*(2*gt22L*Gt231*Gt232 + 4*Gt111*Gt132*gt31L + 4*Gt132*Gt211*gt32L + - 4*Gt221*Gt232*gt32L + 6*Gt132*gt31L*Gt331 + 6*Gt232*gt32L*Gt331 + 6*Gt131*gt31L*Gt332 + 6*Gt231*gt32L*Gt332 + - 4*Gt132*Gt311*gt33L - PDstandard4th12gt33) - gtu13*PDstandard4th13gt33 - gtu22*khalf*PDstandard4th22gt33 - - gtu23*PDstandard4th23gt33 - gtu33*khalf*PDstandard4th33gt33 + - 2*(Gt132*gt21L*Gt231*gtu12 + Gt131*gt21L*Gt232*gtu12 + Gt133*gt21L*Gt231*gtu13 + Gt131*gt21L*Gt233*gtu13 + - gt22L*Gt231*Gt233*gtu13 + Gt132*gt21L*Gt232*gtu22 + gt11L*Gt132*Gt133*gtu23 + Gt132*gt21L*Gt233*gtu23 + - gt22L*Gt232*Gt233*gtu23 + Gt133*gt21L*Gt233*gtu33 + gt31L*PDstandard4th3Xt1) + 2*gt32L*PDstandard4th3Xt2 + - gt33L*(4*Gt231*Gt322*gtu12 + 10*Gt331*Gt332*gtu12 + 4*Gt133*Gt311*gtu13 + 4*Gt231*Gt332*gtu13 + - 10*Gt331*Gt333*gtu13 + 4*Gt132*Gt321*gtu22 + 4*Gt133*Gt321*gtu23 + 4*Gt133*Gt331*gtu33 + 2*PDstandard4th3Xt3) + - 2*Gt231*gt32L*Xt1L + 2*Gt331*gt33L*Xt1L + Gt131* - (2*gt21L*Gt231*gtu11 + 4*Gt111*gt31L*gtu11 + 2*gt11L*Gt132*gtu12 + 2*gt11L*Gt133*gtu13 + 4*Gt132*gt31L*gtu23 + - 4*Gt133*gt31L*gtu33 + 2*gt31L*Xt1L) + 2*Gt132*gt31L*Xt2L + - Gt232*(4*Gt222*gt32L*gtu22 + 6*gt32L*Gt332*gtu22 + 2*Gt133*gt21L*gtu23 + 4*Gt233*gt32L*gtu33 + 2*gt32L*Xt2L) + - Gt332*(4*Gt232*gt33L*gtu23 + 10*Gt333*gt33L*gtu23 + 4*Gt233*gt33L*gtu33 + 2*gt33L*Xt2L) + 2*Gt133*gt31L*Xt3L + - 2*Gt333*gt33L*Xt3L + Gt233*(4*Gt222*gt32L*gtu23 + 6*gt32L*Gt332*gtu23 + 4*Gt132*gt31L*gtu33 + 6*gt32L*Gt333*gtu33 + - 2*gt32L*Xt3L) + gtu11*(4*Gt221*Gt231*gt32L + 6*Gt131*gt31L*Gt331 + 6*Gt231*gt32L*Gt331 + 4*Gt131*Gt311*gt33L + - 4*Gt231*Gt321*gt33L - khalf*PDstandard4th11gt33 + gt11L*SQR(Gt131)) + - 4*(Gt121*Gt231*gt31L*gtu11 + Gt131*Gt211*gt32L*gtu11 + Gt121*Gt131*gt31L*gtu12 + Gt122*Gt231*gt31L*gtu12 + - Gt121*Gt232*gt31L*gtu12 + Gt131*Gt221*gt32L*gtu12 + Gt222*Gt231*gt32L*gtu12 + Gt131*Gt321*gt33L*gtu12 + - Gt232*Gt321*gt33L*gtu12 + Gt111*Gt133*gt31L*gtu13 + Gt132*Gt231*gt31L*gtu13 + Gt121*Gt233*gt31L*gtu13 + - Gt133*Gt211*gt32L*gtu13 + Gt131*Gt231*gt32L*gtu13 + Gt231*Gt232*gt32L*gtu13 + Gt221*Gt233*gt32L*gtu13 + - Gt233*Gt321*gt33L*gtu13 + Gt131*Gt331*gt33L*gtu13 + Gt121*Gt132*gt31L*gtu22 + Gt122*Gt232*gt31L*gtu22 + - Gt132*Gt221*gt32L*gtu22 + Gt232*Gt322*gt33L*gtu22 + Gt121*Gt133*gt31L*gtu23 + Gt132*Gt232*gt31L*gtu23 + - Gt122*Gt233*gt31L*gtu23 + Gt133*Gt221*gt32L*gtu23 + Gt132*Gt231*gt32L*gtu23 + Gt233*Gt322*gt33L*gtu23 + - Gt132*Gt331*gt33L*gtu23 + Gt133*Gt231*gt32L*gtu33 + gt31L*gtu13*SQR(Gt131)) + gt11L*gtu22*SQR(Gt132) + - gt11L*gtu33*SQR(Gt133) + gt22L*gtu11*SQR(Gt231) + gt22L*gtu22*SQR(Gt232) + 4*gt32L*gtu23*SQR(Gt232) + - gt22L*gtu33*SQR(Gt233) + 5*gt33L*gtu11*SQR(Gt331) + 5*gt33L*gtu22*SQR(Gt332) + 5*gt33L*gtu33*SQR(Gt333); - - Rphi11 = 2*(-PDstandard4th11phi - gt11L*gtu11*PDstandard4th11phi - 2*gt11L*gtu12*PDstandard4th12phi - - 2*gt11L*gtu13*PDstandard4th13phi - gt11L*gtu22*PDstandard4th22phi - 2*gt11L*gtu23*PDstandard4th23phi - - gt11L*gtu33*PDstandard4th33phi + Gt311*PDstandard4th3phi + gt11L*Gt311*gtu11*PDstandard4th3phi + - 2*gt11L*Gt321*gtu12*PDstandard4th3phi + 2*gt11L*Gt331*gtu13*PDstandard4th3phi + - gt11L*Gt322*gtu22*PDstandard4th3phi + 2*gt11L*Gt332*gtu23*PDstandard4th3phi + - gt11L*Gt333*gtu33*PDstandard4th3phi + PDstandard4th1phi* - (Gt111 + Gt111*gt11L*gtu11 + 2*gt11L*Gt121*gtu12 + 2*gt11L*Gt131*gtu13 + gt11L*Gt122*gtu22 + - 2*gt11L*Gt132*gtu23 + gt11L*Gt133*gtu33 - 4*gt11L*gtu12*PDstandard4th2phi - 4*gt11L*gtu13*PDstandard4th3phi) + - PDstandard4th2phi*(Gt211 + gt11L*Gt211*gtu11 + - gt11L*(2*Gt221*gtu12 + 2*Gt231*gtu13 + Gt222*gtu22 + 2*Gt232*gtu23 + Gt233*gtu33) - - 4*gt11L*gtu23*PDstandard4th3phi) + (2 - 2*gt11L*gtu11)*SQR(PDstandard4th1phi) - - 2*gt11L*gtu22*SQR(PDstandard4th2phi) - 2*gt11L*gtu33*SQR(PDstandard4th3phi)); - - Rphi21 = 2*(-(gt21L*gtu11*PDstandard4th11phi) - PDstandard4th12phi - 2*gt21L*gtu12*PDstandard4th12phi - - 2*gt21L*gtu13*PDstandard4th13phi - gt21L*gtu22*PDstandard4th22phi - 2*gt21L*gtu23*PDstandard4th23phi - - gt21L*gtu33*PDstandard4th33phi + Gt321*PDstandard4th3phi + gt21L*Gt311*gtu11*PDstandard4th3phi + - 2*gt21L*Gt321*gtu12*PDstandard4th3phi + 2*gt21L*Gt331*gtu13*PDstandard4th3phi + - gt21L*Gt322*gtu22*PDstandard4th3phi + 2*gt21L*Gt332*gtu23*PDstandard4th3phi + - gt21L*Gt333*gtu33*PDstandard4th3phi + PDstandard4th1phi* - (Gt121 + Gt111*gt21L*gtu11 + 2*Gt121*gt21L*gtu12 + 2*Gt131*gt21L*gtu13 + Gt122*gt21L*gtu22 + - 2*Gt132*gt21L*gtu23 + Gt133*gt21L*gtu33 + (2 - 4*gt21L*gtu12)*PDstandard4th2phi - - 4*gt21L*gtu13*PDstandard4th3phi) + PDstandard4th2phi* - (Gt221 + 2*gt21L*Gt221*gtu12 + gt21L*(Gt211*gtu11 + 2*Gt231*gtu13 + Gt222*gtu22 + 2*Gt232*gtu23 + Gt233*gtu33) - - 4*gt21L*gtu23*PDstandard4th3phi) - 2*gt21L*gtu11*SQR(PDstandard4th1phi) - - 2*gt21L*gtu22*SQR(PDstandard4th2phi) - 2*gt21L*gtu33*SQR(PDstandard4th3phi)); - - Rphi31 = 2*(-PDstandard4th13phi + gt31L*(-(gtu11*PDstandard4th11phi) - 2*gtu12*PDstandard4th12phi - - 2*gtu13*PDstandard4th13phi) - gt31L*gtu22*PDstandard4th22phi - 2*gt31L*gtu23*PDstandard4th23phi - - gt31L*gtu33*PDstandard4th33phi + Gt331*PDstandard4th3phi + Gt311*gt31L*gtu11*PDstandard4th3phi + - 2*gt31L*Gt321*gtu12*PDstandard4th3phi + 2*gt31L*Gt331*gtu13*PDstandard4th3phi + - gt31L*Gt322*gtu22*PDstandard4th3phi + 2*gt31L*Gt332*gtu23*PDstandard4th3phi + - gt31L*Gt333*gtu33*PDstandard4th3phi + PDstandard4th1phi* - (Gt131 + Gt111*gt31L*gtu11 + 2*Gt121*gt31L*gtu12 + 2*Gt131*gt31L*gtu13 + Gt122*gt31L*gtu22 + - 2*Gt132*gt31L*gtu23 + Gt133*gt31L*gtu33 - 4*gt31L*gtu12*PDstandard4th2phi + - (2 - 4*gt31L*gtu13)*PDstandard4th3phi) + - PDstandard4th2phi*(Gt231 + 2*Gt231*gt31L*gtu13 + - gt31L*(Gt211*gtu11 + 2*Gt221*gtu12 + Gt222*gtu22 + 2*Gt232*gtu23 + Gt233*gtu33) - - 4*gt31L*gtu23*PDstandard4th3phi) - 2*gt31L*gtu11*SQR(PDstandard4th1phi) - - 2*gt31L*gtu22*SQR(PDstandard4th2phi) - 2*gt31L*gtu33*SQR(PDstandard4th3phi)); - - Rphi22 = 2*(-PDstandard4th22phi + gt22L*(-(gtu11*PDstandard4th11phi) - 2*gtu12*PDstandard4th12phi - - 2*gtu13*PDstandard4th13phi - gtu22*PDstandard4th22phi) - 2*gt22L*gtu23*PDstandard4th23phi - - gt22L*gtu33*PDstandard4th33phi + Gt322*PDstandard4th3phi + gt22L*Gt311*gtu11*PDstandard4th3phi + - 2*gt22L*Gt321*gtu12*PDstandard4th3phi + 2*gt22L*Gt331*gtu13*PDstandard4th3phi + - gt22L*Gt322*gtu22*PDstandard4th3phi + 2*gt22L*Gt332*gtu23*PDstandard4th3phi + - gt22L*Gt333*gtu33*PDstandard4th3phi + PDstandard4th1phi* - (Gt122 + Gt111*gt22L*gtu11 + 2*Gt121*gt22L*gtu12 + 2*Gt131*gt22L*gtu13 + Gt122*gt22L*gtu22 + - 2*Gt132*gt22L*gtu23 + Gt133*gt22L*gtu33 - 4*gt22L*gtu12*PDstandard4th2phi - 4*gt22L*gtu13*PDstandard4th3phi) + - PDstandard4th2phi*(Gt222 + Gt222*gt22L*gtu22 + - gt22L*(Gt211*gtu11 + 2*Gt221*gtu12 + 2*Gt231*gtu13 + 2*Gt232*gtu23 + Gt233*gtu33) - - 4*gt22L*gtu23*PDstandard4th3phi) - 2*gt22L*gtu11*SQR(PDstandard4th1phi) + - (2 - 2*gt22L*gtu22)*SQR(PDstandard4th2phi) - 2*gt22L*gtu33*SQR(PDstandard4th3phi)); - - Rphi32 = 2*(-PDstandard4th23phi + gt32L*(-(gtu11*PDstandard4th11phi) - 2*gtu12*PDstandard4th12phi - - 2*gtu13*PDstandard4th13phi - gtu22*PDstandard4th22phi - 2*gtu23*PDstandard4th23phi) - - gt32L*gtu33*PDstandard4th33phi + Gt332*PDstandard4th3phi + Gt311*gt32L*gtu11*PDstandard4th3phi + - 2*Gt321*gt32L*gtu12*PDstandard4th3phi + 2*gt32L*Gt331*gtu13*PDstandard4th3phi + - Gt322*gt32L*gtu22*PDstandard4th3phi + 2*gt32L*Gt332*gtu23*PDstandard4th3phi + - gt32L*Gt333*gtu33*PDstandard4th3phi + PDstandard4th1phi* - (Gt132 + Gt111*gt32L*gtu11 + 2*Gt121*gt32L*gtu12 + 2*Gt131*gt32L*gtu13 + Gt122*gt32L*gtu22 + - 2*Gt132*gt32L*gtu23 + Gt133*gt32L*gtu33 - 4*gt32L*gtu12*PDstandard4th2phi - 4*gt32L*gtu13*PDstandard4th3phi) + - PDstandard4th2phi*(Gt232 + 2*Gt232*gt32L*gtu23 + - gt32L*(Gt211*gtu11 + 2*Gt221*gtu12 + 2*Gt231*gtu13 + Gt222*gtu22 + Gt233*gtu33) + - (2 - 4*gt32L*gtu23)*PDstandard4th3phi) - 2*gt32L*gtu11*SQR(PDstandard4th1phi) - - 2*gt32L*gtu22*SQR(PDstandard4th2phi) - 2*gt32L*gtu33*SQR(PDstandard4th3phi)); - - Rphi33 = 2*(-PDstandard4th33phi + (Gt333 + gt33L* - (Gt322*gtu22 + 2*(Gt321*gtu12 + Gt331*gtu13 + Gt332*gtu23) + Gt333*gtu33))*PDstandard4th3phi + - PDstandard4th2phi*(Gt233 + gt33L*(Gt211*gtu11 + Gt222*gtu22 + 2*(Gt221*gtu12 + Gt231*gtu13 + Gt232*gtu23) + - Gt233*gtu33 - 4*gtu23*PDstandard4th3phi)) + - PDstandard4th1phi*(Gt133 + gt33L*(Gt111*gtu11 + Gt122*gtu22 + 2*(Gt121*gtu12 + Gt131*gtu13 + Gt132*gtu23) + - Gt133*gtu33 - 4*(gtu12*PDstandard4th2phi + gtu13*PDstandard4th3phi))) + 2*SQR(PDstandard4th3phi) + - gt33L*(-(gtu11*PDstandard4th11phi) - 2*gtu12*PDstandard4th12phi - 2*gtu13*PDstandard4th13phi - - gtu22*PDstandard4th22phi - 2*gtu23*PDstandard4th23phi - gtu33*PDstandard4th33phi + - Gt311*gtu11*PDstandard4th3phi - 2*gtu11*SQR(PDstandard4th1phi) - 2*gtu22*SQR(PDstandard4th2phi) - - 2*gtu33*SQR(PDstandard4th3phi))); - - R11 = Rphi11 + Rt11; - - R21 = Rphi21 + Rt21; - - R31 = Rphi31 + Rt31; - - R22 = Rphi22 + Rt22; - - R32 = Rphi32 + Rt32; - - R33 = Rphi33 + Rt33; - - trR = gu11*R11 + gu22*R22 + 2*(gu12*R21 + gu13*R31 + gu23*R32) + gu33*R33; - - K11 = At11L*e4phi + g11*kthird*trKL; - - K21 = At21L*e4phi + g21*kthird*trKL; - - K31 = At31L*e4phi + g31*kthird*trKL; - - K22 = At22L*e4phi + g22*kthird*trKL; - - K32 = At32L*e4phi + g32*kthird*trKL; - - K33 = At33L*e4phi + g33*kthird*trKL; - - Km11 = gu11*K11 + gu12*K21 + gu13*K31; - - Km21 = gu12*K11 + gu22*K21 + gu23*K31; - - Km31 = gu13*K11 + gu23*K21 + gu33*K31; - - Km12 = gu11*K21 + gu12*K22 + gu13*K32; - - Km22 = gu12*K21 + gu22*K22 + gu23*K32; - - Km32 = gu13*K21 + gu23*K22 + gu33*K32; - - Km13 = gu11*K31 + gu12*K32 + gu13*K33; - - Km23 = gu12*K31 + gu22*K32 + gu23*K33; - - Km33 = gu13*K31 + gu23*K32 + gu33*K33; - - HL = -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - SQR(Km33) + SQR(trKL); - - gK112 = e4phi*(-2*(At11L*G121 + At21L*G221 + At31L*G321) + PDstandard4th2At11 + 4*At11L*PDstandard4th2phi) + - g11*kthird*PDstandard4th2trK; - - gK113 = e4phi*(-2*(At11L*G131 + At21L*G231 + At31L*G331) + PDstandard4th3At11 + 4*At11L*PDstandard4th3phi) + - g11*kthird*PDstandard4th3trK; - - gK211 = e4phi*(-(At11L*G121) - At22L*G211 - At21L*(G111 + G221) - At32L*G311 - At31L*G321 + PDstandard4th1At21 + - 4*At21L*PDstandard4th1phi) + g21*kthird*PDstandard4th1trK; - - gK212 = e4phi*(-(At11L*G122) - At22L*G221 - At21L*(G121 + G222) - At32L*G321 - At31L*G322 + PDstandard4th2At21 + - 4*At21L*PDstandard4th2phi) + g21*kthird*PDstandard4th2trK; - - gK213 = e4phi*(-(At11L*G132) - At22L*G231 - At21L*(G131 + G232) - At32L*G331 - At31L*G332 + PDstandard4th3At21 + - 4*At21L*PDstandard4th3phi) + g21*kthird*PDstandard4th3trK; - - gK311 = e4phi*(-(At11L*G131) - At32L*G211 - At21L*G231 - At33L*G311 - At31L*(G111 + G331) + PDstandard4th1At31 + - 4*At31L*PDstandard4th1phi) + g31*kthird*PDstandard4th1trK; - - gK312 = e4phi*(-(At11L*G132) - At32L*G221 - At21L*G232 - At33L*G321 - At31L*(G121 + G332) + PDstandard4th2At31 + - 4*At31L*PDstandard4th2phi) + g31*kthird*PDstandard4th2trK; - - gK313 = e4phi*(-(At11L*G133) - At32L*G231 - At21L*G233 - At33L*G331 - At31L*(G131 + G333) + PDstandard4th3At31 + - 4*At31L*PDstandard4th3phi) + g31*kthird*PDstandard4th3trK; - - gK221 = e4phi*(-2*(At21L*G121 + At22L*G221 + At32L*G321) + PDstandard4th1At22 + 4*At22L*PDstandard4th1phi) + - g22*kthird*PDstandard4th1trK; - - gK223 = e4phi*(-2*(At21L*G132 + At22L*G232 + At32L*G332) + PDstandard4th3At22 + 4*At22L*PDstandard4th3phi) + - g22*kthird*PDstandard4th3trK; - - gK321 = e4phi*(-(At31L*G121) - At21L*G131 - At32L*G221 - At22L*G231 - At33L*G321 - At32L*G331 + PDstandard4th1At32 + - 4*At32L*PDstandard4th1phi) + g32*kthird*PDstandard4th1trK; - - gK322 = e4phi*(-(At31L*G122) - At21L*G132 - At32L*G222 - At22L*G232 - At33L*G322 - At32L*G332 + PDstandard4th2At32 + - 4*At32L*PDstandard4th2phi) + g32*kthird*PDstandard4th2trK; - - gK323 = e4phi*(-(At31L*G132) - At21L*G133 - At32L*G232 - At22L*G233 - At33L*G332 - At32L*G333 + PDstandard4th3At32 + - 4*At32L*PDstandard4th3phi) + g32*kthird*PDstandard4th3trK; - - gK331 = e4phi*(-2*(At31L*G131 + At32L*G231 + At33L*G331) + PDstandard4th1At33 + 4*At33L*PDstandard4th1phi) + - g33*kthird*PDstandard4th1trK; - - gK332 = e4phi*(-2*(At31L*G132 + At32L*G232 + At33L*G332) + PDstandard4th2At33 + 4*At33L*PDstandard4th2phi) + - g33*kthird*PDstandard4th2trK; - - M1L = (gK112 - gK211)*gu12 + (gK113 - gK311)*gu13 + (gK212 - gK221)*gu22 + (gK213 + gK312 - 2*gK321)*gu23 + - (gK313 - gK331)*gu33; - - M2L = (-gK112 + gK211)*gu11 + (-gK212 + gK221)*gu12 + (gK213 - 2*gK312 + gK321)*gu13 + (gK223 - gK322)*gu23 + - (gK323 - gK332)*gu33; - - M3L = (-gK113 + gK311)*gu11 + (-2*gK213 + gK312 + gK321)*gu12 + (-gK313 + gK331)*gu13 + (-gK223 + gK322)*gu22 + - (-gK323 + gK332)*gu23; - - cSL = Log(detgt); - - cXt1L = Gt111*gtu11 + Gt122*gtu22 + 2*(Gt121*gtu12 + Gt131*gtu13 + Gt132*gtu23) + Gt133*gtu33 - Xt1L; - - cXt2L = Gt211*gtu11 + Gt222*gtu22 + 2*(Gt221*gtu12 + Gt231*gtu13 + Gt232*gtu23) + Gt233*gtu33 - Xt2L; - - cXt3L = Gt311*gtu11 + Gt322*gtu22 + 2*(Gt321*gtu12 + Gt331*gtu13 + Gt332*gtu23) + Gt333*gtu33 - Xt3L; - - cAL = At11L*gtu11 + At22L*gtu22 + 2*(At21L*gtu12 + At31L*gtu13 + At32L*gtu23) + At33L*gtu33; - - - /* Copy local copies back to grid functions */ - cA[index] = cAL; - cS[index] = cSL; - cXt1[index] = cXt1L; - cXt2[index] = cXt2L; - cXt3[index] = cXt3L; - H[index] = HL; - M1[index] = M1L; - M2[index] = M2L; - M3[index] = M3L; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void BSSN_constraints(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &BSSN_constraints_Body); -} diff --git a/ML_BSSN/src/BSSN_convertFromADMBase.c b/ML_BSSN/src/BSSN_convertFromADMBase.c deleted file mode 100644 index 8e9c834..0000000 --- a/ML_BSSN/src/BSSN_convertFromADMBase.c +++ /dev/null @@ -1,295 +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 BSSN_convertFromADMBase_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 BSSN_convertFromADMBase_Body"); - } - - if (cctk_iteration % BSSN_convertFromADMBase_calc_every != BSSN_convertFromADMBase_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 detg = INITVALUE; - CCTK_REAL em4phi = INITVALUE; - CCTK_REAL g11 = INITVALUE, g21 = INITVALUE, g22 = INITVALUE, g31 = INITVALUE, g32 = INITVALUE, g33 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE, gu33 = INITVALUE; - CCTK_REAL K11 = INITVALUE, K21 = INITVALUE, K22 = INITVALUE, K31 = INITVALUE, K32 = INITVALUE, K33 = INITVALUE; - CCTK_REAL Km11 = INITVALUE, Km22 = INITVALUE, Km33 = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL alpL = INITVALUE; - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL At11L = INITVALUE, At21L = INITVALUE, At22L = INITVALUE, At31L = INITVALUE, At32L = INITVALUE, At33L = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL betaxL = INITVALUE; - CCTK_REAL betayL = INITVALUE; - CCTK_REAL betazL = INITVALUE; - CCTK_REAL dtalpL = INITVALUE; - CCTK_REAL dtalphaL = INITVALUE; - CCTK_REAL dtbeta1L = INITVALUE, dtbeta2L = INITVALUE, dtbeta3L = INITVALUE; - CCTK_REAL dtbetaxL = INITVALUE; - CCTK_REAL dtbetayL = INITVALUE; - CCTK_REAL dtbetazL = INITVALUE; - CCTK_REAL gt11L = INITVALUE, gt21L = INITVALUE, gt22L = INITVALUE, gt31L = INITVALUE, gt32L = INITVALUE, gt33L = INITVALUE; - CCTK_REAL gxxL = INITVALUE; - CCTK_REAL gxyL = INITVALUE; - CCTK_REAL gxzL = INITVALUE; - CCTK_REAL gyyL = INITVALUE; - CCTK_REAL gyzL = INITVALUE; - CCTK_REAL gzzL = INITVALUE; - CCTK_REAL kxxL = INITVALUE; - CCTK_REAL kxyL = INITVALUE; - CCTK_REAL kxzL = INITVALUE; - CCTK_REAL kyyL = INITVALUE; - CCTK_REAL kyzL = INITVALUE; - CCTK_REAL kzzL = INITVALUE; - CCTK_REAL phiL = INITVALUE; - CCTK_REAL trKL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - alpL = alp[index]; - betaxL = betax[index]; - betayL = betay[index]; - betazL = betaz[index]; - dtalpL = dtalp[index]; - dtbetaxL = dtbetax[index]; - dtbetayL = dtbetay[index]; - dtbetazL = dtbetaz[index]; - gxxL = gxx[index]; - gxyL = gxy[index]; - gxzL = gxz[index]; - gyyL = gyy[index]; - gyzL = gyz[index]; - gzzL = gzz[index]; - kxxL = kxx[index]; - kxyL = kxy[index]; - kxzL = kxz[index]; - kyyL = kyy[index]; - kyzL = kyz[index]; - kzzL = kzz[index]; - trKL = trK[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 */ - g11 = gxxL; - - g21 = gxyL; - - g31 = gxzL; - - g22 = gyyL; - - g32 = gyzL; - - g33 = gzzL; - - detg = 2*g21*g31*g32 + g33*(g11*g22 - SQR(g21)) - g22*SQR(g31) - g11*SQR(g32); - - gu11 = INV(detg)*(g22*g33 - SQR(g32)); - - gu12 = (g31*g32 - g21*g33)*INV(detg); - - gu13 = (-(g22*g31) + g21*g32)*INV(detg); - - gu22 = INV(detg)*(g11*g33 - SQR(g31)); - - gu23 = (g21*g31 - g11*g32)*INV(detg); - - gu33 = INV(detg)*(g11*g22 - SQR(g21)); - - em4phi = pow(detg,-3); - - phiL = Log(detg)/12.; - - gt11L = em4phi*g11; - - gt21L = em4phi*g21; - - gt31L = em4phi*g31; - - gt22L = em4phi*g22; - - gt32L = em4phi*g32; - - gt33L = em4phi*g33; - - K11 = kxxL; - - K21 = kxyL; - - K31 = kxzL; - - K22 = kyyL; - - K32 = kyzL; - - K33 = kzzL; - - Km11 = gu11*K11 + gu12*K21 + gu13*K31; - - Km22 = gu12*K21 + gu22*K22 + gu23*K32; - - Km33 = gu13*K31 + gu23*K32 + gu33*K33; - - trKL = Km11 + Km22 + Km33; - - At11L = em4phi*(K11 - g11*kthird*trKL); - - At21L = em4phi*(K21 - g21*kthird*trKL); - - At31L = em4phi*(K31 - g31*kthird*trKL); - - At22L = em4phi*(K22 - g22*kthird*trKL); - - At32L = em4phi*(K32 - g32*kthird*trKL); - - At33L = em4phi*(K33 - g33*kthird*trKL); - - alphaL = alpL; - - dtalphaL = dtalpL; - - beta1L = betaxL; - - beta2L = betayL; - - beta3L = betazL; - - dtbeta1L = dtbetaxL; - - dtbeta2L = dtbetayL; - - dtbeta3L = dtbetazL; - - - /* Copy local copies back to grid functions */ - alpha[index] = alphaL; - At11[index] = At11L; - At21[index] = At21L; - At22[index] = At22L; - At31[index] = At31L; - At32[index] = At32L; - At33[index] = At33L; - beta1[index] = beta1L; - beta2[index] = beta2L; - beta3[index] = beta3L; - dtalpha[index] = dtalphaL; - dtbeta1[index] = dtbeta1L; - dtbeta2[index] = dtbeta2L; - dtbeta3[index] = dtbeta3L; - gt11[index] = gt11L; - gt21[index] = gt21L; - gt22[index] = gt22L; - gt31[index] = gt31L; - gt32[index] = gt32L; - gt33[index] = gt33L; - phi[index] = phiL; - trK[index] = trKL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void BSSN_convertFromADMBase(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &BSSN_convertFromADMBase_Body); -} diff --git a/ML_BSSN/src/BSSN_convertFromADMBaseGamma.c b/ML_BSSN/src/BSSN_convertFromADMBaseGamma.c deleted file mode 100644 index 1cd3ffb..0000000 --- a/ML_BSSN/src/BSSN_convertFromADMBaseGamma.c +++ /dev/null @@ -1,254 +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 BSSN_convertFromADMBaseGamma_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 BSSN_convertFromADMBaseGamma_Body"); - } - - if (cctk_iteration % BSSN_convertFromADMBaseGamma_calc_every != BSSN_convertFromADMBaseGamma_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 Gt111 = INITVALUE, Gt121 = INITVALUE, Gt122 = INITVALUE, Gt131 = INITVALUE, Gt132 = INITVALUE, Gt133 = INITVALUE; - CCTK_REAL Gt211 = INITVALUE, Gt221 = INITVALUE, Gt222 = INITVALUE, Gt231 = INITVALUE, Gt232 = INITVALUE, Gt233 = INITVALUE; - CCTK_REAL Gt311 = INITVALUE, Gt321 = INITVALUE, Gt322 = INITVALUE, Gt331 = INITVALUE, Gt332 = INITVALUE, Gt333 = INITVALUE; - CCTK_REAL gtu11 = INITVALUE, gtu12 = INITVALUE, gtu13 = INITVALUE, gtu22 = INITVALUE, gtu23 = INITVALUE, gtu33 = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL gt11L = INITVALUE, gt21L = INITVALUE, gt22L = INITVALUE, gt31L = INITVALUE, gt32L = INITVALUE, gt33L = INITVALUE; - CCTK_REAL Xt1L = INITVALUE, Xt2L = INITVALUE, Xt3L = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - CCTK_REAL PDstandard4th1gt11 = INITVALUE; - CCTK_REAL PDstandard4th2gt11 = INITVALUE; - CCTK_REAL PDstandard4th3gt11 = INITVALUE; - CCTK_REAL PDstandard4th1gt21 = INITVALUE; - CCTK_REAL PDstandard4th2gt21 = INITVALUE; - CCTK_REAL PDstandard4th3gt21 = INITVALUE; - CCTK_REAL PDstandard4th1gt22 = INITVALUE; - CCTK_REAL PDstandard4th2gt22 = INITVALUE; - CCTK_REAL PDstandard4th3gt22 = INITVALUE; - CCTK_REAL PDstandard4th1gt31 = INITVALUE; - CCTK_REAL PDstandard4th2gt31 = INITVALUE; - CCTK_REAL PDstandard4th3gt31 = INITVALUE; - CCTK_REAL PDstandard4th1gt32 = INITVALUE; - CCTK_REAL PDstandard4th2gt32 = INITVALUE; - CCTK_REAL PDstandard4th3gt32 = INITVALUE; - CCTK_REAL PDstandard4th1gt33 = INITVALUE; - CCTK_REAL PDstandard4th2gt33 = INITVALUE; - CCTK_REAL PDstandard4th3gt33 = INITVALUE; - - /* Assign local copies of grid functions */ - gt11L = gt11[index]; - gt21L = gt21[index]; - gt22L = gt22[index]; - gt31L = gt31[index]; - gt32L = gt32[index]; - gt33L = gt33[index]; - - /* Assign local copies of subblock grid functions */ - - /* Include user supplied include files */ - - /* Precompute derivatives (new style) */ - PDstandard4th1gt11 = PDstandard4th1(gt11, i, j, k); - PDstandard4th2gt11 = PDstandard4th2(gt11, i, j, k); - PDstandard4th3gt11 = PDstandard4th3(gt11, i, j, k); - PDstandard4th1gt21 = PDstandard4th1(gt21, i, j, k); - PDstandard4th2gt21 = PDstandard4th2(gt21, i, j, k); - PDstandard4th3gt21 = PDstandard4th3(gt21, i, j, k); - PDstandard4th1gt22 = PDstandard4th1(gt22, i, j, k); - PDstandard4th2gt22 = PDstandard4th2(gt22, i, j, k); - PDstandard4th3gt22 = PDstandard4th3(gt22, i, j, k); - PDstandard4th1gt31 = PDstandard4th1(gt31, i, j, k); - PDstandard4th2gt31 = PDstandard4th2(gt31, i, j, k); - PDstandard4th3gt31 = PDstandard4th3(gt31, i, j, k); - PDstandard4th1gt32 = PDstandard4th1(gt32, i, j, k); - PDstandard4th2gt32 = PDstandard4th2(gt32, i, j, k); - PDstandard4th3gt32 = PDstandard4th3(gt32, i, j, k); - PDstandard4th1gt33 = PDstandard4th1(gt33, i, j, k); - PDstandard4th2gt33 = PDstandard4th2(gt33, i, j, k); - PDstandard4th3gt33 = PDstandard4th3(gt33, i, j, k); - - /* Precompute derivatives (old style) */ - - /* Calculate temporaries and grid functions */ - detgt = 2*gt21L*gt31L*gt32L + gt33L*(gt11L*gt22L - SQR(gt21L)) - gt22L*SQR(gt31L) - gt11L*SQR(gt32L); - - gtu11 = INV(detgt)*(gt22L*gt33L - SQR(gt32L)); - - gtu12 = (gt31L*gt32L - gt21L*gt33L)*INV(detgt); - - gtu13 = (-(gt22L*gt31L) + gt21L*gt32L)*INV(detgt); - - gtu22 = INV(detgt)*(gt11L*gt33L - SQR(gt31L)); - - gtu23 = (gt21L*gt31L - gt11L*gt32L)*INV(detgt); - - gtu33 = INV(detgt)*(gt11L*gt22L - SQR(gt21L)); - - Gt111 = khalf*(gtu11*PDstandard4th1gt11 + 2*(gtu12*PDstandard4th1gt21 + gtu13*PDstandard4th1gt31) - - gtu12*PDstandard4th2gt11 - gtu13*PDstandard4th3gt11); - - Gt211 = khalf*(gtu12*PDstandard4th1gt11 + 2*(gtu22*PDstandard4th1gt21 + gtu23*PDstandard4th1gt31) - - gtu22*PDstandard4th2gt11 - gtu23*PDstandard4th3gt11); - - Gt311 = khalf*(gtu13*PDstandard4th1gt11 + 2*(gtu23*PDstandard4th1gt21 + gtu33*PDstandard4th1gt31) - - gtu23*PDstandard4th2gt11 - gtu33*PDstandard4th3gt11); - - Gt121 = khalf*(gtu12*PDstandard4th1gt22 + gtu11*PDstandard4th2gt11 + - gtu13*(PDstandard4th1gt32 + PDstandard4th2gt31 - PDstandard4th3gt21)); - - Gt221 = khalf*(gtu22*PDstandard4th1gt22 + gtu12*PDstandard4th2gt11 + - gtu23*(PDstandard4th1gt32 + PDstandard4th2gt31 - PDstandard4th3gt21)); - - Gt321 = khalf*(gtu23*PDstandard4th1gt22 + gtu13*PDstandard4th2gt11 + - gtu33*(PDstandard4th1gt32 + PDstandard4th2gt31 - PDstandard4th3gt21)); - - Gt131 = khalf*(gtu13*PDstandard4th1gt33 + gtu11*PDstandard4th3gt11 + - gtu12*(PDstandard4th1gt32 - PDstandard4th2gt31 + PDstandard4th3gt21)); - - Gt231 = khalf*(gtu23*PDstandard4th1gt33 + gtu12*PDstandard4th3gt11 + - gtu22*(PDstandard4th1gt32 - PDstandard4th2gt31 + PDstandard4th3gt21)); - - Gt331 = khalf*(gtu33*PDstandard4th1gt33 + gtu13*PDstandard4th3gt11 + - gtu23*(PDstandard4th1gt32 - PDstandard4th2gt31 + PDstandard4th3gt21)); - - Gt122 = khalf*(gtu11*(-PDstandard4th1gt22 + 2*PDstandard4th2gt21) + gtu12*PDstandard4th2gt22 + - gtu13*(2*PDstandard4th2gt32 - PDstandard4th3gt22)); - - Gt222 = khalf*(gtu12*(-PDstandard4th1gt22 + 2*PDstandard4th2gt21) + gtu22*PDstandard4th2gt22 + - gtu23*(2*PDstandard4th2gt32 - PDstandard4th3gt22)); - - Gt322 = khalf*(gtu13*(-PDstandard4th1gt22 + 2*PDstandard4th2gt21) + gtu23*PDstandard4th2gt22 + - gtu33*(2*PDstandard4th2gt32 - PDstandard4th3gt22)); - - Gt132 = khalf*(gtu13*PDstandard4th2gt33 + gtu11*(-PDstandard4th1gt32 + PDstandard4th2gt31 + PDstandard4th3gt21) + - gtu12*PDstandard4th3gt22); - - Gt232 = khalf*(gtu23*PDstandard4th2gt33 + gtu12*(-PDstandard4th1gt32 + PDstandard4th2gt31 + PDstandard4th3gt21) + - gtu22*PDstandard4th3gt22); - - Gt332 = khalf*(gtu33*PDstandard4th2gt33 + gtu13*(-PDstandard4th1gt32 + PDstandard4th2gt31 + PDstandard4th3gt21) + - gtu23*PDstandard4th3gt22); - - Gt133 = khalf*(-(gtu11*PDstandard4th1gt33) - gtu12*PDstandard4th2gt33 + 2*gtu11*PDstandard4th3gt31 + - 2*gtu12*PDstandard4th3gt32 + gtu13*PDstandard4th3gt33); - - Gt233 = khalf*(-(gtu12*PDstandard4th1gt33) - gtu22*PDstandard4th2gt33 + 2*gtu12*PDstandard4th3gt31 + - 2*gtu22*PDstandard4th3gt32 + gtu23*PDstandard4th3gt33); - - Gt333 = khalf*(-(gtu13*PDstandard4th1gt33) - gtu23*PDstandard4th2gt33 + 2*gtu13*PDstandard4th3gt31 + - 2*gtu23*PDstandard4th3gt32 + gtu33*PDstandard4th3gt33); - - Xt1L = Gt111*gtu11 + Gt122*gtu22 + 2*(Gt121*gtu12 + Gt131*gtu13 + Gt132*gtu23) + Gt133*gtu33; - - Xt2L = Gt211*gtu11 + Gt222*gtu22 + 2*(Gt221*gtu12 + Gt231*gtu13 + Gt232*gtu23) + Gt233*gtu33; - - Xt3L = Gt311*gtu11 + Gt322*gtu22 + 2*(Gt321*gtu12 + Gt331*gtu13 + Gt332*gtu23) + Gt333*gtu33; - - - /* Copy local copies back to grid functions */ - Xt1[index] = Xt1L; - Xt2[index] = Xt2L; - Xt3[index] = Xt3L; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void BSSN_convertFromADMBaseGamma(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &BSSN_convertFromADMBaseGamma_Body); -} diff --git a/ML_BSSN/src/BSSN_convertToADMBase.c b/ML_BSSN/src/BSSN_convertToADMBase.c deleted file mode 100644 index 39177ba..0000000 --- a/ML_BSSN/src/BSSN_convertToADMBase.c +++ /dev/null @@ -1,267 +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 BSSN_convertToADMBase_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 BSSN_convertToADMBase_Body"); - } - - if (cctk_iteration % BSSN_convertToADMBase_calc_every != BSSN_convertToADMBase_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 e4phi = INITVALUE; - CCTK_REAL g11 = INITVALUE, g21 = INITVALUE, g22 = INITVALUE, g31 = INITVALUE, g32 = INITVALUE, g33 = INITVALUE; - CCTK_REAL K11 = INITVALUE, K21 = INITVALUE, K22 = INITVALUE, K31 = INITVALUE, K32 = INITVALUE, K33 = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL alpL = INITVALUE; - CCTK_REAL alphaL = INITVALUE; - CCTK_REAL At11L = INITVALUE, At21L = INITVALUE, At22L = INITVALUE, At31L = INITVALUE, At32L = INITVALUE, At33L = INITVALUE; - CCTK_REAL beta1L = INITVALUE, beta2L = INITVALUE, beta3L = INITVALUE; - CCTK_REAL betaxL = INITVALUE; - CCTK_REAL betayL = INITVALUE; - CCTK_REAL betazL = INITVALUE; - CCTK_REAL dtalpL = INITVALUE; - CCTK_REAL dtalphaL = INITVALUE; - CCTK_REAL dtbeta1L = INITVALUE, dtbeta2L = INITVALUE, dtbeta3L = INITVALUE; - CCTK_REAL dtbetaxL = INITVALUE; - CCTK_REAL dtbetayL = INITVALUE; - CCTK_REAL dtbetazL = INITVALUE; - CCTK_REAL gt11L = INITVALUE, gt21L = INITVALUE, gt22L = INITVALUE, gt31L = INITVALUE, gt32L = INITVALUE, gt33L = INITVALUE; - CCTK_REAL gxxL = INITVALUE; - CCTK_REAL gxyL = INITVALUE; - CCTK_REAL gxzL = INITVALUE; - CCTK_REAL gyyL = INITVALUE; - CCTK_REAL gyzL = INITVALUE; - CCTK_REAL gzzL = INITVALUE; - CCTK_REAL kxxL = INITVALUE; - CCTK_REAL kxyL = INITVALUE; - CCTK_REAL kxzL = INITVALUE; - CCTK_REAL kyyL = INITVALUE; - CCTK_REAL kyzL = INITVALUE; - CCTK_REAL kzzL = INITVALUE; - CCTK_REAL phiL = INITVALUE; - CCTK_REAL trKL = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - alphaL = alpha[index]; - At11L = At11[index]; - At21L = At21[index]; - At22L = At22[index]; - At31L = At31[index]; - At32L = At32[index]; - At33L = At33[index]; - beta1L = beta1[index]; - beta2L = beta2[index]; - beta3L = beta3[index]; - dtalphaL = dtalpha[index]; - dtbeta1L = dtbeta1[index]; - dtbeta2L = dtbeta2[index]; - dtbeta3L = dtbeta3[index]; - gt11L = gt11[index]; - gt21L = gt21[index]; - gt22L = gt22[index]; - gt31L = gt31[index]; - gt32L = gt32[index]; - gt33L = gt33[index]; - phiL = phi[index]; - trKL = trK[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 */ - e4phi = exp(4*phiL); - - g11 = e4phi*gt11L; - - g21 = e4phi*gt21L; - - g31 = e4phi*gt31L; - - g22 = e4phi*gt22L; - - g32 = e4phi*gt32L; - - g33 = e4phi*gt33L; - - gxxL = g11; - - gxyL = g21; - - gxzL = g31; - - gyyL = g22; - - gyzL = g32; - - gzzL = g33; - - K11 = At11L*e4phi + g11*kthird*trKL; - - K21 = At21L*e4phi + g21*kthird*trKL; - - K31 = At31L*e4phi + g31*kthird*trKL; - - K22 = At22L*e4phi + g22*kthird*trKL; - - K32 = At32L*e4phi + g32*kthird*trKL; - - K33 = At33L*e4phi + g33*kthird*trKL; - - kxxL = K11; - - kxyL = K21; - - kxzL = K31; - - kyyL = K22; - - kyzL = K32; - - kzzL = K33; - - alpL = alphaL; - - dtalpL = dtalphaL; - - betaxL = beta1L; - - betayL = beta2L; - - betazL = beta3L; - - dtbetaxL = dtbeta1L; - - dtbetayL = dtbeta2L; - - dtbetazL = dtbeta3L; - - - /* Copy local copies back to grid functions */ - alp[index] = alpL; - betax[index] = betaxL; - betay[index] = betayL; - betaz[index] = betazL; - dtalp[index] = dtalpL; - dtbetax[index] = dtbetaxL; - dtbetay[index] = dtbetayL; - dtbetaz[index] = dtbetazL; - gxx[index] = gxxL; - gxy[index] = gxyL; - gxz[index] = gxzL; - gyy[index] = gyyL; - gyz[index] = gyzL; - gzz[index] = gzzL; - kxx[index] = kxxL; - kxy[index] = kxyL; - kxz[index] = kxzL; - kyy[index] = kyyL; - kyz[index] = kyzL; - kzz[index] = kzzL; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void BSSN_convertToADMBase(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &BSSN_convertToADMBase_Body); -} diff --git a/ML_BSSN/src/BSSN_enforce.c b/ML_BSSN/src/BSSN_enforce.c deleted file mode 100644 index e6a22f1..0000000 --- a/ML_BSSN/src/BSSN_enforce.c +++ /dev/null @@ -1,179 +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 BSSN_enforce_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 BSSN_enforce_Body"); - } - - if (cctk_iteration % BSSN_enforce_calc_every != BSSN_enforce_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 */ - LC_LOOP3 (somename, - i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], - max[0]-min[0],max[1]-min[1],max[2]-min[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 gtu11 = INITVALUE, gtu12 = INITVALUE, gtu13 = INITVALUE, gtu22 = INITVALUE, gtu23 = INITVALUE, gtu33 = INITVALUE; - CCTK_REAL trA = INITVALUE; - - /* Declare local copies of grid functions */ - CCTK_REAL At11L = INITVALUE, At21L = INITVALUE, At22L = INITVALUE, At31L = INITVALUE, At32L = INITVALUE, At33L = INITVALUE; - CCTK_REAL gt11L = INITVALUE, gt21L = INITVALUE, gt22L = INITVALUE, gt31L = INITVALUE, gt32L = INITVALUE, gt33L = INITVALUE; - /* Declare precomputed derivatives*/ - - /* Declare derivatives */ - - /* Assign local copies of grid functions */ - At11L = At11[index]; - At21L = At21[index]; - At22L = At22[index]; - At31L = At31[index]; - At32L = At32[index]; - At33L = At33[index]; - gt11L = gt11[index]; - gt21L = gt21[index]; - gt22L = gt22[index]; - gt31L = gt31[index]; - gt32L = gt32[index]; - gt33L = gt33[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 */ - detgt = 2*gt21L*gt31L*gt32L + gt33L*(gt11L*gt22L - SQR(gt21L)) - gt22L*SQR(gt31L) - gt11L*SQR(gt32L); - - gtu11 = INV(detgt)*(gt22L*gt33L - SQR(gt32L)); - - gtu12 = (gt31L*gt32L - gt21L*gt33L)*INV(detgt); - - gtu13 = (-(gt22L*gt31L) + gt21L*gt32L)*INV(detgt); - - gtu22 = INV(detgt)*(gt11L*gt33L - SQR(gt31L)); - - gtu23 = (gt21L*gt31L - gt11L*gt32L)*INV(detgt); - - gtu33 = INV(detgt)*(gt11L*gt22L - SQR(gt21L)); - - trA = At11L*gtu11 + At22L*gtu22 + 2*(At21L*gtu12 + At31L*gtu13 + At32L*gtu23) + At33L*gtu33; - - At11L = At11L - gt11L*kthird*trA; - - At21L = At21L - gt21L*kthird*trA; - - At31L = At31L - gt31L*kthird*trA; - - At22L = At22L - gt22L*kthird*trA; - - At32L = At32L - gt32L*kthird*trA; - - At33L = At33L - gt33L*kthird*trA; - - - /* Copy local copies back to grid functions */ - At11[index] = At11L; - At21[index] = At21L; - At22[index] = At22L; - At31[index] = At31L; - At32[index] = At32L; - At33[index] = At33L; - - /* Copy local copies back to subblock grid functions */ - } - LC_ENDLOOP3 (somename); -} - -void BSSN_enforce(CCTK_ARGUMENTS) -{ - DECLARE_CCTK_ARGUMENTS - DECLARE_CCTK_PARAMETERS - - GenericFD_LoopOverEverything(cctkGH, &BSSN_enforce_Body); -} diff --git a/McLachlan.m b/McLachlan.m deleted file mode 100644 index 8762230..0000000 --- a/McLachlan.m +++ /dev/null @@ -1,544 +0,0 @@ -(* < McLachlan.m /Applications/Mathematica.app/Contents/MacOS/MathKernel | tee McLachlan.out *) - -$Path = Join[$Path, {"~/Calpha/Kranc-devel/Tools/CodeGen", - "~/Calpha/Kranc-devel/Tools/MathematicaMisc"}]; - -Get["KrancThorn`"]; - -SetEnhancedTimes[False]; -SetSourceLanguage["C"]; - -(******************************************************************************) -(* Derivatives *) -(******************************************************************************) - -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] - - (* - PDstandard6th[i_] -> StandardCenteredDifferenceOperator[1,3,i], - PDstandard6th[i_, i_] -> StandardCenteredDifferenceOperator[2,3,i], - PDstandard6th[i_, j_] -> StandardCenteredDifferenceOperator[1,3,i] - StandardCenteredDifferenceOperator[1,3,j] - *) -}; - -PD = PDstandard4th; - -KD = KroneckerDelta; - -(******************************************************************************) -(* Tensors *) -(******************************************************************************) - -(* Register the tensor quantities with the TensorTools package *) -Map [DefineTensor, - {g, K, alpha, beta, H, M, detg, gu, G, R, trR, Km, trK, - phi, gt, At, Xt, dtalpha, dtbeta, trA, cXt, cS, cA, - e4phi, em4phi, gtu, ddetg, ddetgt, Gt, Rt, Rphi, gK}]; -(* SetTensorAttribute[g, TensorWeight, 0]; *) - -Map [AssertSymmetricDecreasing, - {g[la,lb], K[la,lb], R[la,lb], - gt[la,lb], At[la,lb], Rt[la,lb], Rphi[la,lb]}]; -AssertSymmetricDecreasing [G[ua,lb,lc], lb, lc]; -AssertSymmetricDecreasing [Gt[ua,lb,lc], lb, lc]; -AssertSymmetricDecreasing [gK[la,lb,lc], la, lb]; -Map [AssertSymmetricIncreasing, {gu[ua,ub]}]; -Map [AssertSymmetricIncreasing, {gtu[ua,ub]}]; - -DefineConnection [CD, PD, G]; -DefineConnection [CDt, PD, Gt]; - -Map [DefineTensor, - {gxx, gxy, gxz, gyy, gyz, gzz, - kxx, kxy, kxz, kyy, kyz, kzz, - alp, - dtalp, - betax, betay, betaz, - dtbetax, dtbetay, dtbetaz}]; - -(******************************************************************************) -(* Expressions *) -(******************************************************************************) - -detgExpr = Det [MatrixOfComponents [g [la,lb]]]; -ddetgExpr[la_] = - Sum [D[Det[MatrixOfComponents[g[la, lb]]], X] PD[X, la], - {X, Union[Flatten[MatrixOfComponents[g[la, lb]]]]}]; - -detgtExpr = Det [MatrixOfComponents [gt[la,lb]]]; -ddetgtExpr[la_] = - Sum [D[Det[MatrixOfComponents[gt[la, lb]]], X] PD[X, la], - {X, Union[Flatten[MatrixOfComponents[gt[la, lb]]]]}]; - -(******************************************************************************) -(* Groups *) -(******************************************************************************) - -evolvedGroups = - {SetGroupName [CreateGroupFromTensor [g[la,lb]], "metric"], - SetGroupName [CreateGroupFromTensor [K[la,lb]], "curv" ], - SetGroupName [CreateGroupFromTensor [alpha ], "lapse" ], - SetGroupName [CreateGroupFromTensor [beta[ua]], "shift" ]}; -evaluatedGroups = - {SetGroupName [CreateGroupFromTensor [H ], "Ham"], - SetGroupName [CreateGroupFromTensor [M[la]], "mom"]}; - -declaredGroups = Join [evolvedGroups, evaluatedGroups]; -declaredGroupNames = Map [First, declaredGroups]; - -evolvedGroupsBSSN = - {SetGroupName [CreateGroupFromTensor [phi ], "log_confac"], - SetGroupName [CreateGroupFromTensor [gt[la,lb] ], "metric" ], - SetGroupName [CreateGroupFromTensor [Xt[ua] ], "Gamma" ], - SetGroupName [CreateGroupFromTensor [trK ], "trace_curv"], - SetGroupName [CreateGroupFromTensor [At[la,lb] ], "curv" ], - SetGroupName [CreateGroupFromTensor [alpha ], "lapse" ], - SetGroupName [CreateGroupFromTensor [dtalpha ], "dtlapse" ], - SetGroupName [CreateGroupFromTensor [beta[ua] ], "shift" ], - SetGroupName [CreateGroupFromTensor [dtbeta[ua]], "dtshift" ]}; -evaluatedGroupsBSSN = - {SetGroupName [CreateGroupFromTensor [H ], "Ham"], - SetGroupName [CreateGroupFromTensor [M[la] ], "mom"], - SetGroupName [CreateGroupFromTensor [cS ], "cons_detg"], - SetGroupName [CreateGroupFromTensor [cXt[ua]], "cons_Gamma"], - SetGroupName [CreateGroupFromTensor [cA ], "cons_traceA"]}; - -declaredGroupsBSSN = Join [evolvedGroupsBSSN, evaluatedGroupsBSSN]; -declaredGroupNamesBSSN = Map [First, declaredGroupsBSSN]; - - - -extraGroups = - {{"ADMBase::metric", {gxx, gxy, gxz, gyy, gyz, gzz}}, - {"ADMBase::curv", {kxx, kxy, kxz, kyy, kyz, kzz}}, - {"ADMBase::lapse", {alp}}, - {"ADMBase::dtlapse", {dtalp}}, - {"ADMBase::shift", {betax, betay, betaz}}, - {"ADMBase::dtshift", {dtbetax, dtbetay, dtbetaz}}}; - - - -groups = Join [declaredGroups, extraGroups]; -groupsBSSN = Join [declaredGroupsBSSN, extraGroups]; - -(******************************************************************************) -(* Initial data *) -(******************************************************************************) - -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], - K[la,lb] -> 0, - alpha -> 1, - beta[ua] -> 0 - } -} - -initialCalcBSSN = -{ - Name -> "ML_BSSN_Minkowski", - Schedule -> {"IN ADMBase_InitialData"}, - ConditionalOnKeyword -> {"my_initial_data", "Minkowski"}, - (* Where -> Boundary, *) - (* Where -> Interior, *) - Equations -> - { - phi -> 0, - gt[la,lb] -> KD[la,lb], - trK -> 0, - At[la,lb] -> 0, - Xt[ua] -> 0, - alpha -> 1, - dtalpha -> 0, - beta[ua] -> 0, - dtbeta[ua] -> 0 - } -} - -(******************************************************************************) -(* Convert from ADMBase *) -(******************************************************************************) - -convertFromADMBaseCalc = -{ - Name -> "ML_ADM_convertFromADMBase", - Schedule -> {"AT initial", "AFTER ADMBase_PostInitial"}, - ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, - Equations -> - { - g11 -> gxx, - g21 -> gxy, - g31 -> gxz, - g22 -> gyy, - g32 -> gyz, - g33 -> gzz, - K11 -> kxx, - K21 -> kxy, - K31 -> kxz, - K22 -> kyy, - K32 -> kyz, - K33 -> kzz, - (* TODO: this is incomplete; it ignores dtalp and dtbeta^i *) - alpha -> alp, - beta1 -> betax, - beta2 -> betay, - beta3 -> betaz - } -} - -convertFromADMBaseCalcBSSN = -{ - Name -> "ML_BSSN_convertFromADMBase", - Schedule -> {"AT initial", "AFTER ADMBase_PostInitial"}, - ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, - Shorthands -> {g[la,lb], detg, gu[ua,ub], em4phi, K[la,lb], Km[ua,lb]}, - Equations -> - { - g11 -> gxx, - g21 -> gxy, - g31 -> gxz, - g22 -> gyy, - g32 -> gyz, - g33 -> gzz, - - detg -> detgExpr, - gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], - - (* em4phi -> Exp [-4 phi], *) - em4phi -> 1 / detg^3, - phi -> Log [detg] / 12, - gt[la,lb] -> em4phi g[la,lb], - - K11 -> kxx, - K21 -> kxy, - K31 -> kxz, - K22 -> kyy, - K32 -> kyz, - K33 -> kzz, - - Km[ua,lb] -> gu[ua,uc] K[lc,lb], - trK -> Km[ua,la], - At[la,lb] -> em4phi (K[la,lb] - (1/3) g[la,lb] trK), - - alpha -> alp, - dtalpha -> dtalp, - - beta1 -> betax, - beta2 -> betay, - beta3 -> betaz, - (* TODO: this is wrong *) - dtbeta1 -> dtbetax, - dtbeta2 -> dtbetay, - dtbeta3 -> dtbetaz - } -} - -convertFromADMBaseCalcBSSNGamma = -{ - Name -> "ML_BSSN_convertFromADMBaseGamma", - Schedule -> {"AT initial", "AFTER ML_BSSN_convertFromADMBase"}, - ConditionalOnKeyword -> {"my_initial_data", "ADMBase"}, - Where -> Interior, - Shorthands -> {detgt, gtu[ua,ub], Gt[ua,lb,lc]}, - Equations -> - { - detgt -> detgtExpr, - gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], - Gt[ua,lb,lc] -> 1/2 gtu[ua,ud] - (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]), - Xt[ua] -> gtu[ub,uc] Gt[ua,lb,lc] - } -} - -(******************************************************************************) -(* Convert to ADMBase *) -(******************************************************************************) - -convertToADMBaseCalc = -{ - Name -> "ML_ADM_convertToADMBase", - Schedule -> {"IN MoL_PostStep", "AFTER ADM_ApplyBoundConds"}, - Equations -> - { - gxx -> g11, - gxy -> g21, - gxz -> g31, - gyy -> g22, - gyz -> g32, - gzz -> g33, - kxx -> K11, - kxy -> K21, - kxz -> K31, - kyy -> K22, - kyz -> K32, - kzz -> K33, - (* TODO: this is wrong; it sets dtalp and dtbeta^i incorrectly *) - alp -> alpha, - dtalp -> 0, - betax -> beta1, - betay -> beta2, - betaz -> beta3, - dtbetax -> 0, - dtbetay -> 0, - dtbetaz -> 0 - } -} - -convertToADMBaseCalcBSSN = -{ - Name -> "ML_BSSN_convertToADMBase", - Schedule -> {"IN MoL_PostStep", "AFTER ML_ADM_ApplyBoundConds"}, - Shorthands -> {e4phi, g[la,lb], K[la,lb]}, - Equations -> - { - e4phi -> Exp [4 phi], - g[la,lb] -> e4phi gt[la,lb], - gxx -> g11, - gxy -> g21, - gxz -> g31, - gyy -> g22, - gyz -> g32, - gzz -> g33, - K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, - kxx -> K11, - kxy -> K21, - kxz -> K31, - kyy -> K22, - kyz -> K32, - kzz -> K33, - alp -> alpha, - dtalp -> dtalpha, - betax -> beta1, - betay -> beta2, - betaz -> beta3, - dtbetax -> dtbeta1, - dtbetay -> dtbeta2, - dtbetaz -> dtbeta3 - } -} - -(******************************************************************************) -(* Evolution equations *) -(******************************************************************************) - -evolCalc = -{ - Name -> "ML_ADM_RHS", - Schedule -> {"IN MoL_CalcRHS"}, - Shorthands -> {detg, gu[ua,ub], G[ua,lb,lc], R[la,lb], Km[ua,lb], trK}, - Equations -> - { - detg -> detgExpr, - gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], - G[ua,lb,lc] -> 1/2 gu[ua,ud] - (PD[g[lb,ld],lc] + PD[g[lc,ld],lb] - PD[g[lb,lc],ld]), - R[la,lb] -> G[u1,l2,la] G[u2,l1,lb] - G[u1,la,lb] G[u2,l1,l2] - + 1/2 gu[u1,u2] (- PD[g[l1,l2],la,lb] + PD[g[l1,la],l2,lb] - - PD[g[la,lb],l1,l2] + PD[g[l2,lb],l1,la]), - Km[ua,lb] -> gu[ua,uc] K[lc,lb], - trK -> Km[ua,la], - - dot[g[la,lb]] -> -2 alpha K[la,lb] - + Lie[g[la,lb], beta], - dot[K[la,lb]] -> - CD[alpha,la,lb] - + alpha (+ R[la,lb] + K[la,lb] trK - 2 K[la,lc] Km[uc,lb]) - + Lie[K[la,lb], beta], - dot[alpha] -> 0, - dot[beta[ua]] -> 0 - } -} - -enforceCalcBSSN = -{ - Name -> "ML_BSSN_enforce", - Schedule -> {"IN MoL_PostStep"}, - Shorthands -> {detgt, gtu[ua,ub], trA}, - Equations -> - { - detgt -> detgtExpr, - gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], - - trA -> gtu[ua,ub] At[la,lb], - - At[la,lb] -> At[la,lb] - (1/3) gt[la,lb] trA - } -} - -(******************************************************************************) -(* Constraint equations *) -(******************************************************************************) - -constraintsCalc = -{ - Name -> "ML_ADM_constraints", - Schedule -> {"AT analysis"}, - Where -> Interior, - Shorthands -> {detg, gu[ua,ub], G[ua,lb,lc], R[la,lb], trR, Km[ua,lb], trK}, - Equations -> - { - detg -> detgExpr, - gu[ua,ub] -> 1/detg detgExpr MatrixInverse[g[ua,ub]], - G[ua,lb,lc] -> 1/2 gu[ua,ud] - (PD[g[lb,ld],lc] + PD[g[lc,ld],lb] - PD[g[lb,lc],ld]), - R[la,lb] -> G[u1,l2,la] G[l1,lb,u2] - G[u1,la,lb] G[l1,l2,u2] - + 1/2 gu[u1,u2] (- PD[g[l1,l2],la,lb] + PD[g[l1,la],l2,lb] - - PD[g[la,lb],l1,l2] + PD[g[l2,lb],l1,la]), - trR -> R[la,lb] gu[ua,ub], - Km[ua,lb] -> gu[ua,uc] K[lc,lb], - trK -> Km[ua,la], - - H -> trR - Km[ua,lb] Km[ub,la] + trK^2, - M[la] -> gu[ub,uc] (CD[K[lc,la], lb] - CD[K[lc,lb], la]) - } -} - -constraintsCalcBSSN = -{ - Name -> "ML_BSSN_constraints", - Schedule -> {"AT analysis"}, - Where -> Interior, - Shorthands -> {detgt, ddetgt[la], gtu[ua,ub], Gt[ua,lb,lc], e4phi, - g[la,lb], detg, gu[ua,ub], ddetg[la], G[ua,lb,lc], - Rt[la,lb], Rphi[la,lb], R[la,lb], trR, - K[la,lb], Km[la,lb], gK[la,lb,lc]}, - Equations -> - { - detgt -> detgtExpr, - ddetgt[la] -> ddetgtExpr[la], - gtu[ua,ub] -> 1/detgt detgtExpr MatrixInverse [gt[ua,ub]], - Gt[ua,lb,lc] -> 1/2 gtu[ua,ud] - (PD[gt[lb,ld],lc] + PD[gt[lc,ld],lb] - PD[gt[lb,lc],ld]), - - e4phi -> Exp [4 phi], - g[la,lb] -> e4phi gt[la,lb], - detg -> detgExpr, - gu[ua,ub] -> 1/detg detgExpr MatrixInverse [g[ua,ub]], - (* ddetg[la] -> PD[e4phi detg,la], *) - ddetg[la] -> e4phi ddetgt[la] + 4 detgt e4phi PD[phi,la], - G[ua,lb,lc] -> Gt[ua,lb,lc] - + 1/(2 detg) (+ KD[ua,lb] ddetg[lc] + KD[ua,lc] ddetg[lb] - - (1/3) g[lb,lc] gu[ua,ud] ddetg[ld]), - - (* PRD 62, 044034 (2000), eqn. (18) *) - 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] - + Xt[uk] gt[li,ln] Gt[un,lj,lk] + Xt[uk] gt[lj,ln] Gt[un,li,lk] - + gtu[ul,um] (+ 2 Gt[uk,ll,li] gt[lj,ln] Gt[un,lk,lm] - + 2 Gt[uk,ll,lj] gt[li,ln] Gt[un,lk,lm] - + Gt[uk,li,lm] gt[lk,ln] Gt[un,ll,lj]), - (* PRD 62, 044034 (2000), eqn. (15) *) - (* TODO: Check that CDt takes the tensor weight of phi into account *) - Rphi[li,lj] -> - 2 CDt[phi,lj,li] - - 2 gt[li,lj] gtu[ul,un] CDt[phi,ll,ln] - + 4 CDt[phi,li] CDt[phi,lj] - - 4 gt[li,lj] gtu[ul,un] CDt[phi,ln] CDt[phi,ll], - - R[la,lb] -> Rt[la,lb] + Rphi[la,lb], - trR -> gu[ua,ub] R[la,lb], - - K[la,lb] -> e4phi At[la,lb] + (1/3) g[la,lb] trK, - Km[ua,lb] -> gu[ua,uc] K[lc,lb], - - H -> trR - Km[ua,lb] Km[ub,la] + trK^2, - - (* 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]), - - (* det gamma-tilde *) - cS -> Log [detgt], - - (* Gamma constraint *) - cXt[ua] -> gtu[ub,uc] Gt[ua,lb,lc] - Xt[ua], - - (* trace A-tilde *) - cA -> gtu[ua,ub] At[la,lb] - } -} - -(******************************************************************************) -(* Implementations *) -(******************************************************************************) - -inheritedImplementations = {"ADMBase"}; - -(******************************************************************************) -(* Parameters *) -(******************************************************************************) - -inheritedKeywordParameters = {"ADMBase::initial_data"}; - -keywordParameters = -{{ - Name -> "my_initial_data", - (* Visibility -> "restricted", *) - (* Description -> "ddd", *) - AllowedValues -> {"ADMBase", "Minkowski"}, - Default -> "ADMBase" -}}; - -realParameters = {}; - -(******************************************************************************) -(* Construct the thorns *) -(******************************************************************************) - -calculations = -{ - initialCalc, - convertFromADMBaseCalc, - evolCalc, - convertToADMBaseCalc, - constraintsCalc -}; - -CreateKrancThornTT [groups, ".", "ML_ADM", - Calculations -> calculations, - DeclaredGroups -> declaredGroupNames, - PartialDerivatives -> derivatives, - InheritedImplementations -> inheritedImplementations, - InheritedKeywordParameters -> inheritedKeywordParameters, - KeywordParameters -> keywordParameters, - RealParameters -> realParameters -]; - -calculationsBSSN = -{ - initialCalcBSSN, - convertFromADMBaseCalcBSSN, - convertFromADMBaseCalcBSSNGamma, - enforceCalcBSSN, - convertToADMBaseCalcBSSN, - constraintsCalcBSSN -}; - -CreateKrancThornTT [groupsBSSN, ".", "ML_BSSN", - Calculations -> calculationsBSSN, - DeclaredGroups -> declaredGroupNamesBSSN, - PartialDerivatives -> derivatives, - InheritedImplementations -> inheritedImplementations, - InheritedKeywordParameters -> inheritedKeywordParameters, - KeywordParameters -> keywordParameters, - RealParameters -> realParameters -]; diff --git a/McLachlan.sh b/McLachlan.sh deleted file mode 100755 index fdf0810..0000000 --- a/McLachlan.sh +++ /dev/null @@ -1,6 +0,0 @@ -#! /bin/bash - -MATHEMATICA="/Applications/Mathematica.app/Contents/MacOS/MathKernel" - -# Run Mathematica to regenerate the code -< McLachlan.m "$MATHEMATICA" | tee McLachlan.out diff --git a/m/McLachlan.out b/m/McLachlan.out deleted file mode 100644 index 07b5f1e..0000000 --- a/m/McLachlan.out +++ /dev/null @@ -1,755 +0,0 @@ -Mathematica 6.0 for Mac OS X x86 (32-bit) -Copyright 1988-2007 Wolfram Research, Inc. - -In[1]:= -In[2]:= -In[2]:= -In[3]:= -In[3]:= -Set::wrsym: Symbol Quiet is Protected. - -General::obspkg: - LinearAlgebra`MatrixManipulation` - is now obsolete. The legacy version being loaded may conflict with - current Mathematica functionality. See the Compatibility Guide for - updating information. - -In[4]:= -In[4]:= -In[5]:= User set source language to C - -In[6]:= -In[6]:= -In[7]:= -In[8]:= -In[9]:= -In[9]:= -Syntax::newl: The newline character after - " PDstandard4th[i_, j_] -> StandardCe<<11>>renceOperator[1,2,i]" - is understood as a multiplication operator. - - -In[10]:= -In[10]:= -In[11]:= -In[11]:= -In[12]:= -In[12]:= -In[13]:= -In[14]:= -In[15]:= -In[15]:= -In[16]:= -In[17]:= -In[18]:= -In[18]:= -In[19]:= -In[20]:= -In[21]:= -In[22]:= -In[23]:= -In[24]:= -In[24]:= -In[25]:= -In[26]:= -In[26]:= -In[27]:= -In[27]:= -In[28]:= -In[29]:= -In[30]:= -In[30]:= -In[31]:= -In[32]:= -In[32]:= -In[33]:= -In[34]:= -In[34]:= -In[35]:= -In[36]:= -In[37]:= -In[37]:= -In[38]:= -In[39]:= -In[39]:= -In[40]:= -In[41]:= -In[41]:= -In[42]:= -In[43]:= -In[43]:= -In[44]:= -In[45]:= -In[45]:= -In[45]:= -In[45]:= -In[46]:= -In[46]:= -In[46]:= -In[46]:= -In[47]:= -In[48]:= -In[48]:= -In[49]:= -In[50]:= -In[51]:= -In[51]:= -Out[51]= {Name -> ML_ADM_Minkowski, Schedule -> {IN ADMBase_InitialData}, - -> ConditionalOnKeyword -> {my_initial_data, Minkowski}, - -> Equations -> - -> {g -> KroneckerDelta[ , ], K -> 0, alpha -> 1, - a b a b a b - - a -> beta -> 0}} - -In[52]:= -In[52]:= -Out[52]= {Name -> ML_BSSN_Minkowski, Schedule -> {IN ADMBase_InitialData}, - -> ConditionalOnKeyword -> {my_initial_data, Minkowski}, - -> Equations -> - -> {phi -> 0, gt -> KroneckerDelta[ , ], trK -> 0, - a b a b - - a a -> At -> 0, Xt -> 0, alpha -> 1, dtalpha -> 0, beta -> 0, - a b - - a -> dtbeta -> 0}} - -In[53]:= -In[53]:= -In[54]:= -In[55]:= -In[56]:= -In[56]:= -Out[56]= {Name -> ML_ADM_convertFromADMBase, - -> Schedule -> {AT initial, AFTER ADMBase_PostInitial}, - -> ConditionalOnKeyword -> {my_initial_data, ADMBase}, - -> Equations -> - -> {g11 -> gxx, g21 -> gxy, g31 -> gxz, g22 -> gyy, g32 -> gyz, - -> g33 -> gzz, K11 -> kxx, K21 -> kxy, K31 -> kxz, K22 -> kyy, - -> K32 -> kyz, K33 -> kzz, alpha -> alp, beta1 -> betax, beta2 -> betay, - -> beta3 -> betaz}} - -In[57]:= -In[57]:= -Syntax::bktwrn: - "em4phi (K[la,lb] - (1/3) g[la,lb] trK)" represents multiplication; use - "em4phi [K[la,lb] - (1/3) g[la,lb] trK]" to represent a function. - - -Out[57]= {Name -> ML_BSSN_convertFromADMBase, - -> Schedule -> {AT initial, AFTER ADMBase_PostInitial}, - -> ConditionalOnKeyword -> {my_initial_data, ADMBase}, - -> Shorthands -> - - a b a -> {g , detg, gu , em4phi, K , Km }, - a b a b b - -> Equations -> - -> {g11 -> gxx, g21 -> gxy, g31 -> gxz, g22 -> gyy, g32 -> gyz, - -> g33 -> gzz, detg -> - - 2 2 2 -> -(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + g11 g22 g33, - - a b -> gu -> - - 2 2 2 -> ((-(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + g11 g22 g33) - - a b -3 -> MatrixInverse[g ]) / detg, em4phi -> detg , - - Log[detg] -> phi -> ---------, gt -> em4phi g , K11 -> kxx, - 12 a b a b - -> K21 -> kxy, K31 -> kxz, K22 -> kyy, K32 -> kyz, K33 -> kzz, - - a a c a -> Km -> K gu , trK -> Km , - b c b a - - -trK -> At -> em4phi (---- g + K ), alpha -> alp, - a b 3 a b a b - -> dtalpha -> dtalp, beta1 -> betax, beta2 -> betay, beta3 -> betaz, - -> dtbeta1 -> dtbetax, dtbeta2 -> dtbetay, dtbeta3 -> dtbetaz}} - -In[58]:= -In[58]:= -Out[58]= {Name -> ML_BSSN_convertFromADMBaseGamma, - -> Schedule -> {AT initial, AFTER ML_BSSN_convertFromADMBase}, - -> ConditionalOnKeyword -> {my_initial_data, ADMBase}, Where -> Interior, - - a b a -> Shorthands -> {detgt, gtu , Gt }, - b c - -> Equations -> - -> {detgt -> - - 2 2 2 -> -(gt22 gt31 ) + 2 gt21 gt31 gt32 - gt11 gt32 - gt21 gt33 + - - a b -> gt11 gt22 gt33, gtu -> - - 2 2 2 -> ((-(gt22 gt31 ) + 2 gt21 gt31 gt32 - gt11 gt32 - gt21 gt33 + - - a b -> gt11 gt22 gt33) MatrixInverse[gt ]) / detgt, - - a -> Gt -> - b c - -> (-PDstandard4th[gt , ] + PDstandard4th[gt , ] + - b c d b d c - - a d -> PDstandard4th[gt , ]) / 2 gtu , - c d b - - a a b c -> Xt -> Gt gtu }} - b c - -In[59]:= -In[59]:= -In[60]:= -In[61]:= -In[62]:= -In[62]:= -Out[62]= {Name -> ML_ADM_convertToADMBase, - -> Schedule -> {IN MoL_PostStep, AFTER ADM_ApplyBoundConds}, - -> Equations -> - -> {gxx -> g11, gxy -> g21, gxz -> g31, gyy -> g22, gyz -> g32, - -> gzz -> g33, kxx -> K11, kxy -> K21, kxz -> K31, kyy -> K22, - -> kyz -> K32, kzz -> K33, alp -> alpha, dtalp -> 0, betax -> beta1, - -> betay -> beta2, betaz -> beta3, dtbetax -> 0, dtbetay -> 0, - -> dtbetaz -> 0}} - -In[63]:= -In[63]:= -Out[63]= {Name -> ML_BSSN_convertToADMBase, - -> Schedule -> {IN MoL_PostStep, AFTER ML_ADM_ApplyBoundConds}, - -> Shorthands -> {e4phi, g , K }, - a b a b - -> Equations -> - - 4 phi -> {e4phi -> E , g -> e4phi gt , gxx -> g11, - a b a b - -> gxy -> g21, gxz -> g31, gyy -> g22, gyz -> g32, gzz -> g33, - - trK -> K -> e4phi At + --- g , kxx -> K11, kxy -> K21, - a b a b 3 a b - -> kxz -> K31, kyy -> K22, kyz -> K32, kzz -> K33, alp -> alpha, - -> dtalp -> dtalpha, betax -> beta1, betay -> beta2, betaz -> beta3, - -> dtbetax -> dtbeta1, dtbetay -> dtbeta2, dtbetaz -> dtbeta3}} - -In[64]:= -In[64]:= -In[65]:= -In[66]:= -In[67]:= -In[67]:= -Syntax::bktwrn: - "alpha (+ R[la,lb] + K[la,lb] trK - 2 K[la,lc] Km[uc,lb])" - represents multiplication; use - "alpha [+ R[la,lb] + K[la,lb] trK - 2 K[la,lc] Km[uc,lb]]" - to represent a function. - - -Out[67]= {Name -> ML_ADM_RHS, Schedule -> {IN MoL_CalcRHS}, - -> Shorthands -> - - a b a a -> {detg, gu , G , R , Km , trK}, - b c a b b - -> Equations -> - - 2 2 2 -> {detg -> -(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + - - a b -> g11 g22 g33, gu -> - - 2 2 2 -> ((-(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + g11 g22 g33) - - a b -> MatrixInverse[g ]) / detg, - - a -> G -> - b c - -> (-PDstandard4th[g , ] + PDstandard4th[g , ] + - b c d b d c - - a d -> PDstandard4th[g , ]) / 2 gu , - c d b - -> R -> - a b - -> (-PDstandard4th[g , , ] + - 1 2 a b - -> PDstandard4th[g , , ] + - 1 a 2 b - -> PDstandard4th[g , , ] - - 2 b 1 a - - 1 2 -> PDstandard4th[g , , ]) / 2 gu - - a b 1 2 - - 2 1 2 1 -> G G + G G , - 1 2 a b 1 b 2 a - - a a c a -> Km -> K gu , trK -> Km , - b c b a - -> dot[g ] -> Lie[g , beta] - 2 alpha K , - a b a b a b - -> dot[K ] -> - a b - -> -CD[alpha, , ] + Lie[K , beta] + - a b a b - - c -> alpha (trK K - 2 Km K + R ), - a b b a c a b - - a -> dot[alpha] -> 0, dot[beta ] -> 0}} - -In[68]:= -In[68]:= -Out[68]= {Name -> ML_BSSN_enforce, Schedule -> {IN MoL_PostStep}, - - a b -> Shorthands -> {detgt, gtu , trA}, - -> Equations -> - -> {detgt -> - - 2 2 2 -> -(gt22 gt31 ) + 2 gt21 gt31 gt32 - gt11 gt32 - gt21 gt33 + - - a b -> gt11 gt22 gt33, gtu -> - - 2 2 2 -> ((-(gt22 gt31 ) + 2 gt21 gt31 gt32 - gt11 gt32 - gt21 gt33 + - - a b -> gt11 gt22 gt33) MatrixInverse[gt ]) / detgt, - - a b -> trA -> gtu At , - a b - - trA -> At -> At - --- gt }} - a b a b 3 a b - -In[69]:= -In[69]:= -In[70]:= -In[71]:= -In[72]:= -In[72]:= -Out[72]= {Name -> ML_ADM_constraints, Schedule -> {AT analysis}, - -> Where -> Interior, Shorthands -> - - a b a a -> {detg, gu , G , R , trR, Km , trK}, - b c a b b - -> Equations -> - - 2 2 2 -> {detg -> -(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + - - a b -> g11 g22 g33, gu -> - - 2 2 2 -> ((-(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + g11 g22 g33) - - a b -> MatrixInverse[g ]) / detg, - - a -> G -> - b c - -> (-PDstandard4th[g , ] + PDstandard4th[g , ] + - b c d b d c - - a d -> PDstandard4th[g , ]) / 2 gu , - c d b - -> R -> - a b - -> (-PDstandard4th[g , , ] + - 1 2 a b - -> PDstandard4th[g , , ] + - 1 a 2 b - -> PDstandard4th[g , , ] - - 2 b 1 a - - 1 2 -> PDstandard4th[g , , ]) / 2 gu + - a b 1 2 - - 1 2 1 2 -> G G - G G , - 2 a 1 b a b 1 2 - - a b a a c -> trR -> R gu , Km -> K gu , - a b b c b - - a 2 b a -> trK -> Km , H -> trK + trR - Km Km , - a a b - - b c -> M -> CD[K , ] - CD[K , ] gu }} - a c a b c b a - -In[73]:= -In[73]:= -Out[73]= {Name -> ML_BSSN_constraints, Schedule -> {AT analysis}, - -> Where -> Interior, Shorthands -> - - a b a -> {detgt, ddetgt , gtu , Gt , e4phi, g , - a b c a b - - a b a -> detg, gu , ddetg , G , Rt , Rphi , - a b c a b a b - -> R , trR, K , Km , gK }, - a b a b a b a b c - -> Equations -> - -> {detgt -> - - 2 2 2 -> -(gt22 gt31 ) + 2 gt21 gt31 gt32 - gt11 gt32 - gt21 gt33 + - -> gt11 gt22 gt33, ddetgt -> - a - - 2 -> (-gt32 + gt22 gt33) PDstandard4th[gt11, ] + - a - -> (2 gt31 gt32 - 2 gt21 gt33) PDstandard4th[gt21, ] + - a - - 2 -> (-gt31 + gt11 gt33) PDstandard4th[gt22, ] + - a - -> (-2 gt22 gt31 + 2 gt21 gt32) PDstandard4th[gt31, ] + - a - -> (2 gt21 gt31 - 2 gt11 gt32) PDstandard4th[gt32, ] + - a - - 2 -> (-gt21 + gt11 gt22) PDstandard4th[gt33, ], - a - - a b -> gtu -> - - 2 2 2 -> ((-(gt22 gt31 ) + 2 gt21 gt31 gt32 - gt11 gt32 - gt21 gt33 + - - a b -> gt11 gt22 gt33) MatrixInverse[gt ]) / detgt, - - a -> Gt -> - b c - -> (-PDstandard4th[gt , ] + PDstandard4th[gt , ] + - b c d b d c - - a d -> PDstandard4th[gt , ]) / 2 gtu , - c d b - - 4 phi -> e4phi -> E , g -> e4phi gt , - a b a b - - 2 2 2 -> detg -> -(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + - - a b -> g11 g22 g33, gu -> - - 2 2 2 -> ((-(g22 g31 ) + 2 g21 g31 g32 - g11 g32 - g21 g33 + g11 g22 g33) - - a b -> MatrixInverse[g ]) / detg, - -> ddetg -> - a - -> 4 detgt e4phi PDstandard4th[phi, ] + e4phi ddetgt , - a a - - a -> G -> - b c - - a -> (KroneckerDelta[ , ] ddetg + - c b - - a -> KroneckerDelta[ , ] ddetg - - b c - - 1 a d a -> - gu g ddetg ) / (2 detg) + Gt , - 3 b c d b c - -> Rt -> - i j - - k -> PDstandard4th[Xt , ] gt + - j k i - - k -> PDstandard4th[Xt , ] gt - - i k j - - PDstandard4th[gt , , ] - i j l m l m -> ------------------------------------- gtu + - 2 - - n k -> Gt gt Xt + - i k j n - - n k -> Gt gt Xt + - j k i n - - n k -> 2 Gt Gt gt + - k m l i j n - - n k -> 2 Gt Gt gt + - k m l j i n - - n k l m -> Gt Gt gt gtu , - l j i m k n - -> Rphi -> - i j - -> 4 CDt[phi, ] CDt[phi, ] - 2 CDt[phi, , ] - - i j j i - - l n -> 4 CDt[phi, ] CDt[phi, ] gtu gt - - l n i j - - l n -> 2 CDt[phi, , ] gtu gt , - l n i j - - a b -> R -> Rphi + Rt , trR -> R gu , - a b a b a b a b - - trK -> K -> e4phi At + --- g , - a b a b 3 a b - - a a c -> Km -> K gu , - b c b - - 2 b a -> H -> trK + trR - Km Km , - a b - -> gK -> - a b c - -> e4phi CD[At , ] + - a b c - -> 4 e4phi PDstandard4th[phi, ] At + - c a b - - PDstandard4th[trK, ] - c -> ------------------------ g , - 3 a b - - b c -> M -> gK - gK gu , cS -> Log[detgt], - a c a b c b a - - a a a b c -> cXt -> -1 Xt + Gt gtu , - b c - - a b -> cA -> gtu At }} - a b - -In[74]:= -In[74]:= -In[75]:= -In[76]:= -In[77]:= -In[77]:= -In[78]:= -In[78]:= -In[79]:= -In[80]:= -In[81]:= -In[81]:= -In[82]:= -In[82]:= -In[83]:= -In[83]:= -In[84]:= -In[84]:= -In[85]:= -In[86]:= -In[87]:= -In[87]:= -In[88]:= -In[88]:= Processing tensorial arguments -Computing reflection symmetries -Creating (component-based) Kranc thorn -Processing arguments to CreateKrancThorn -Verifying arguments -Creating startup file -Creating configuration file -Creating interface file -Creating param file -Creating schedule file -Creating MoL registration file -Creating symmetry registration file -Creating differencing header file -Creating calculation source files -Creating calculation function: ML_ADM_Minkowski -Creating calculation function: ML_ADM_convertFromADMBase -Creating calculation function: ML_ADM_RHS -Creating calculation function: ML_ADM_convertToADMBase -Creating calculation function: ML_ADM_constraints -Creating make file -Creating thorn -Creating thorns in directory ./ML_ADM - -In[89]:= -In[89]:= -In[90]:= -In[90]:= Processing tensorial arguments -Computing reflection symmetries -Creating (component-based) Kranc thorn -Processing arguments to CreateKrancThorn -Verifying arguments -Creating startup file -Creating configuration file -Creating interface file -Creating param file -Creating schedule file -Creating MoL registration file -Creating symmetry registration file -Creating differencing header file -Creating calculation source files -Creating calculation function: ML_BSSN_Minkowski -Creating calculation function: ML_BSSN_convertFromADMBase -Creating calculation function: ML_BSSN_convertFromADMBaseGamma -Creating calculation function: ML_BSSN_enforce -Creating calculation function: ML_BSSN_convertToADMBase -Creating calculation function: ML_BSSN_constraints -Creating make file -Creating thorn -Creating thorns in directory ./ML_BSSN - -In[91]:= |