From 70fd92da159408e939381abf14f37c01d523760d Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Wed, 2 Apr 2008 18:42:43 -0500 Subject: Create several specialised thorns: ML_ADM and ML_BSSN as standard cases, ML_BSSN_MP for multi-patch simulations, and ML_BSSN_M for AMR simulations with matter. Add helper thorns for all BSSN thorns. Helper thorns are required at run time, e.g. to register gauge conditions. Split ADM and BSSN equations into their own source files. --- ML_ADM/interface.ccl | 42 +++--- ML_ADM/param.ccl | 40 +++--- ML_ADM/schedule.ccl | 62 ++++----- ML_ADM/src/Boundaries.c | 216 +++++++++++++++---------------- ML_ADM/src/ML_ADM_Minkowski.c | 5 +- ML_ADM/src/ML_ADM_constraints_boundary.c | 5 +- ML_ADM/src/ML_ADM_convertFromADMBase.c | 5 +- ML_ADM/src/ML_ADM_convertToADMBase.c | 5 +- 8 files changed, 196 insertions(+), 184 deletions(-) (limited to 'ML_ADM') diff --git a/ML_ADM/interface.ccl b/ML_ADM/interface.ccl index 79752b5..4b49535 100644 --- a/ML_ADM/interface.ccl +++ b/ML_ADM/interface.ccl @@ -5,7 +5,7 @@ implements: ML_ADM -inherits: ADMBase TmunuBase Grid GenericFD Boundary +inherits: ADMBase Grid GenericFD Boundary @@ -31,21 +31,21 @@ CCTK_INT FUNCTION Boundary_SelectVarForBC(CCTK_POINTER_TO_CONST IN GH, CCTK_INT USES FUNCTION Boundary_SelectVarForBC public: -CCTK_REAL Ham type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' +CCTK_REAL ML_Ham type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' { H -} "Ham" +} "ML_Ham" public: -CCTK_REAL mom type=GF timelevels=1 tags='tensortypealias="D" tensorweight=1.0000000000000000000' +CCTK_REAL ML_mom type=GF timelevels=1 tags='tensortypealias="D" tensorweight=1.0000000000000000000' { M1, M2, M3 -} "mom" +} "ML_mom" public: -CCTK_REAL ml_curvrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' +CCTK_REAL ML_curvrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' { K11rhs, K12rhs, @@ -53,16 +53,16 @@ CCTK_REAL ml_curvrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorw K22rhs, K23rhs, K33rhs -} "ml_curvrhs" +} "ML_curvrhs" public: -CCTK_REAL ml_lapserhs type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' +CCTK_REAL ML_lapserhs type=GF timelevels=1 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' { alpharhs -} "ml_lapserhs" +} "ML_lapserhs" public: -CCTK_REAL ml_metricrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' +CCTK_REAL ML_metricrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' { g11rhs, g12rhs, @@ -70,18 +70,18 @@ CCTK_REAL ml_metricrhs type=GF timelevels=1 tags='tensortypealias="DD_sym" tenso g22rhs, g23rhs, g33rhs -} "ml_metricrhs" +} "ML_metricrhs" public: -CCTK_REAL ml_shiftrhs type=GF timelevels=1 tags='tensortypealias="U" tensorweight=1.0000000000000000000' +CCTK_REAL ML_shiftrhs type=GF timelevels=1 tags='tensortypealias="U" tensorweight=1.0000000000000000000' { beta1rhs, beta2rhs, beta3rhs -} "ml_shiftrhs" +} "ML_shiftrhs" public: -CCTK_REAL ml_curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' +CCTK_REAL ML_curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' { K11, K12, @@ -89,16 +89,16 @@ CCTK_REAL ml_curv type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweig K22, K23, K33 -} "ml_curv" +} "ML_curv" public: -CCTK_REAL ml_lapse type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' +CCTK_REAL ML_lapse type=GF timelevels=3 tags='tensortypealias="Scalar" tensorweight=1.0000000000000000000' { alpha -} "ml_lapse" +} "ML_lapse" public: -CCTK_REAL ml_metric type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' +CCTK_REAL ML_metric type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorweight=1.0000000000000000000' { g11, g12, @@ -106,12 +106,12 @@ CCTK_REAL ml_metric type=GF timelevels=3 tags='tensortypealias="DD_sym" tensorwe g22, g23, g33 -} "ml_metric" +} "ML_metric" public: -CCTK_REAL ml_shift type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000' +CCTK_REAL ML_shift type=GF timelevels=3 tags='tensortypealias="U" tensorweight=1.0000000000000000000' { beta1, beta2, beta3 -} "ml_shift" +} "ML_shift" diff --git a/ML_ADM/param.ccl b/ML_ADM/param.ccl index 52675c8..660f718 100644 --- a/ML_ADM/param.ccl +++ b/ML_ADM/param.ccl @@ -19,9 +19,9 @@ USES CCTK_INT MoL_Num_Evolved_Vars USES CCTK_INT MoL_Num_Constrained_Vars restricted: -CCTK_REAL verbose "verbose" +CCTK_INT verbose "verbose" { - "*:*" :: "" + *:* :: "" } 0 private: @@ -47,8 +47,8 @@ CCTK_INT ML_ADM_MaxNumEvolvedVars "Number of evolved variables used by this thor restricted: CCTK_INT ML_ADM_MaxNumConstrainedVars "Number of constrained variables used by this thorn" ACCUMULATOR-BASE=MethodofLines::MoL_Num_Constrained_Vars { - 38:38 :: "Number of constrained variables used by this thorn" -} 38 + 65:65 :: "Number of constrained variables used by this thorn" +} 65 restricted: CCTK_INT timelevels "Number of active timelevels" @@ -333,7 +333,7 @@ KEYWORD beta3_bound "Boundary condition to implement" } "skip" private: -KEYWORD ml_curv_bound "Boundary condition to implement" +KEYWORD ML_curv_bound "Boundary condition to implement" { "flat" :: "Flat boundary condition" "none" :: "No boundary condition" @@ -345,7 +345,7 @@ KEYWORD ml_curv_bound "Boundary condition to implement" } "skip" private: -KEYWORD ml_lapse_bound "Boundary condition to implement" +KEYWORD ML_lapse_bound "Boundary condition to implement" { "flat" :: "Flat boundary condition" "none" :: "No boundary condition" @@ -357,7 +357,7 @@ KEYWORD ml_lapse_bound "Boundary condition to implement" } "skip" private: -KEYWORD ml_metric_bound "Boundary condition to implement" +KEYWORD ML_metric_bound "Boundary condition to implement" { "flat" :: "Flat boundary condition" "none" :: "No boundary condition" @@ -369,7 +369,7 @@ KEYWORD ml_metric_bound "Boundary condition to implement" } "skip" private: -KEYWORD ml_shift_bound "Boundary condition to implement" +KEYWORD ML_shift_bound "Boundary condition to implement" { "flat" :: "Flat boundary condition" "none" :: "No boundary condition" @@ -477,25 +477,25 @@ CCTK_REAL beta3_bound_speed "characteristic speed at boundary" } 1. private: -CCTK_REAL ml_curv_bound_speed "characteristic speed at boundary" +CCTK_REAL ML_curv_bound_speed "characteristic speed at boundary" { "0:*" :: "outgoing characteristic speed > 0" } 1. private: -CCTK_REAL ml_lapse_bound_speed "characteristic speed at boundary" +CCTK_REAL ML_lapse_bound_speed "characteristic speed at boundary" { "0:*" :: "outgoing characteristic speed > 0" } 1. private: -CCTK_REAL ml_metric_bound_speed "characteristic speed at boundary" +CCTK_REAL ML_metric_bound_speed "characteristic speed at boundary" { "0:*" :: "outgoing characteristic speed > 0" } 1. private: -CCTK_REAL ml_shift_bound_speed "characteristic speed at boundary" +CCTK_REAL ML_shift_bound_speed "characteristic speed at boundary" { "0:*" :: "outgoing characteristic speed > 0" } 1. @@ -597,25 +597,25 @@ CCTK_REAL beta3_bound_limit "limit value for r -> infinity" } 0. private: -CCTK_REAL ml_curv_bound_limit "limit value for r -> infinity" +CCTK_REAL ML_curv_bound_limit "limit value for r -> infinity" { "*:*" :: "value of limit value is unrestricted" } 0. private: -CCTK_REAL ml_lapse_bound_limit "limit value for r -> infinity" +CCTK_REAL ML_lapse_bound_limit "limit value for r -> infinity" { "*:*" :: "value of limit value is unrestricted" } 0. private: -CCTK_REAL ml_metric_bound_limit "limit value for r -> infinity" +CCTK_REAL ML_metric_bound_limit "limit value for r -> infinity" { "*:*" :: "value of limit value is unrestricted" } 0. private: -CCTK_REAL ml_shift_bound_limit "limit value for r -> infinity" +CCTK_REAL ML_shift_bound_limit "limit value for r -> infinity" { "*:*" :: "value of limit value is unrestricted" } 0. @@ -717,25 +717,25 @@ CCTK_REAL beta3_bound_scalar "Dirichlet boundary value" } 0. private: -CCTK_REAL ml_curv_bound_scalar "Dirichlet boundary value" +CCTK_REAL ML_curv_bound_scalar "Dirichlet boundary value" { "*:*" :: "unrestricted" } 0. private: -CCTK_REAL ml_lapse_bound_scalar "Dirichlet boundary value" +CCTK_REAL ML_lapse_bound_scalar "Dirichlet boundary value" { "*:*" :: "unrestricted" } 0. private: -CCTK_REAL ml_metric_bound_scalar "Dirichlet boundary value" +CCTK_REAL ML_metric_bound_scalar "Dirichlet boundary value" { "*:*" :: "unrestricted" } 0. private: -CCTK_REAL ml_shift_bound_scalar "Dirichlet boundary value" +CCTK_REAL ML_shift_bound_scalar "Dirichlet boundary value" { "*:*" :: "unrestricted" } 0. diff --git a/ML_ADM/schedule.ccl b/ML_ADM/schedule.ccl index 5dbca42..e5e3c84 100644 --- a/ML_ADM/schedule.ccl +++ b/ML_ADM/schedule.ccl @@ -4,68 +4,68 @@ # Mathematica script written by Ian Hinder and Sascha Husa -STORAGE: Ham[1] +STORAGE: ML_Ham[1] -STORAGE: mom[1] +STORAGE: ML_mom[1] -STORAGE: ml_curvrhs[1] +STORAGE: ML_curvrhs[1] -STORAGE: ml_lapserhs[1] +STORAGE: ML_lapserhs[1] -STORAGE: ml_metricrhs[1] +STORAGE: ML_metricrhs[1] -STORAGE: ml_shiftrhs[1] +STORAGE: ML_shiftrhs[1] if (timelevels == 1) { - STORAGE: ml_curv[1] + STORAGE: ML_curv[1] } if (timelevels == 2) { - STORAGE: ml_curv[2] + STORAGE: ML_curv[2] } if (timelevels == 3) { - STORAGE: ml_curv[3] + STORAGE: ML_curv[3] } if (timelevels == 1) { - STORAGE: ml_lapse[1] + STORAGE: ML_lapse[1] } if (timelevels == 2) { - STORAGE: ml_lapse[2] + STORAGE: ML_lapse[2] } if (timelevels == 3) { - STORAGE: ml_lapse[3] + STORAGE: ML_lapse[3] } if (timelevels == 1) { - STORAGE: ml_metric[1] + STORAGE: ML_metric[1] } if (timelevels == 2) { - STORAGE: ml_metric[2] + STORAGE: ML_metric[2] } if (timelevels == 3) { - STORAGE: ml_metric[3] + STORAGE: ML_metric[3] } if (timelevels == 1) { - STORAGE: ml_shift[1] + STORAGE: ML_shift[1] } if (timelevels == 2) { - STORAGE: ml_shift[2] + STORAGE: ML_shift[2] } if (timelevels == 3) { - STORAGE: ml_shift[3] + STORAGE: ML_shift[3] } schedule ML_ADM_Startup at STARTUP @@ -112,10 +112,10 @@ schedule ML_ADM_RHS IN MoL_CalcRHS schedule ML_ADM_RHS AT analysis { LANG: C - SYNC: ml_curvrhs - SYNC: ml_lapserhs - SYNC: ml_metricrhs - SYNC: ml_shiftrhs + SYNC: ML_curvrhs + SYNC: ML_lapserhs + SYNC: ML_metricrhs + SYNC: ML_shiftrhs } "ML_ADM_RHS" @@ -127,7 +127,7 @@ if (CCTK_EQUALS(my_boundary_condition, "Minkowski")) } "ML_ADM_boundary" } -schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER (ML_ADM_ApplyBCs ML_ADM_boundary) +schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER ML_ADM_ApplyBCs { LANG: C } "ML_ADM_convertToADMBase" @@ -135,10 +135,10 @@ schedule ML_ADM_convertToADMBase IN MoL_PostStep AFTER (ML_ADM_ApplyBCs ML_ADM_b schedule ML_ADM_constraints AT analysis { LANG: C - SYNC: Ham - SYNC: mom - TRIGGERS: Ham - TRIGGERS: mom + SYNC: ML_Ham + SYNC: ML_mom + TRIGGERS: ML_Ham + TRIGGERS: ML_mom } "ML_ADM_constraints" schedule ML_ADM_constraints_boundary AT analysis AFTER ML_ADM_constraints @@ -150,10 +150,10 @@ schedule ML_ADM_ApplyBoundConds in MoL_PostStep { LANG: C OPTIONS: level - SYNC: ml_curv - SYNC: ml_lapse - SYNC: ml_metric - SYNC: ml_shift + SYNC: ML_curv + SYNC: ML_lapse + SYNC: ML_metric + SYNC: ML_shift } "apply boundary conditions" schedule ML_ADM_CheckBoundaries at BASEGRID diff --git a/ML_ADM/src/Boundaries.c b/ML_ADM/src/Boundaries.c index 956080a..c2971da 100644 --- a/ML_ADM/src/Boundaries.c +++ b/ML_ADM/src/Boundaries.c @@ -35,48 +35,48 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS) CCTK_INT ierr = 0; - if (CCTK_EQUALS(ml_curv_bound, "none" ) || - CCTK_EQUALS(ml_curv_bound, "static") || - CCTK_EQUALS(ml_curv_bound, "flat" ) || - CCTK_EQUALS(ml_curv_bound, "zero" ) ) + if (CCTK_EQUALS(ML_curv_bound, "none" ) || + CCTK_EQUALS(ML_curv_bound, "static") || + CCTK_EQUALS(ML_curv_bound, "flat" ) || + CCTK_EQUALS(ML_curv_bound, "zero" ) ) { ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, - "ML_ADM::ml_curv", ml_curv_bound); + "ML_ADM::ML_curv", ML_curv_bound); if (ierr < 0) - CCTK_WARN(0, "Failed to register ml_curv_bound BC for ML_ADM::ml_curv!"); + CCTK_WARN(0, "Failed to register ML_curv_bound BC for ML_ADM::ML_curv!"); } - if (CCTK_EQUALS(ml_lapse_bound, "none" ) || - CCTK_EQUALS(ml_lapse_bound, "static") || - CCTK_EQUALS(ml_lapse_bound, "flat" ) || - CCTK_EQUALS(ml_lapse_bound, "zero" ) ) + if (CCTK_EQUALS(ML_lapse_bound, "none" ) || + CCTK_EQUALS(ML_lapse_bound, "static") || + CCTK_EQUALS(ML_lapse_bound, "flat" ) || + CCTK_EQUALS(ML_lapse_bound, "zero" ) ) { ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, - "ML_ADM::ml_lapse", ml_lapse_bound); + "ML_ADM::ML_lapse", ML_lapse_bound); if (ierr < 0) - CCTK_WARN(0, "Failed to register ml_lapse_bound BC for ML_ADM::ml_lapse!"); + CCTK_WARN(0, "Failed to register ML_lapse_bound BC for ML_ADM::ML_lapse!"); } - if (CCTK_EQUALS(ml_metric_bound, "none" ) || - CCTK_EQUALS(ml_metric_bound, "static") || - CCTK_EQUALS(ml_metric_bound, "flat" ) || - CCTK_EQUALS(ml_metric_bound, "zero" ) ) + if (CCTK_EQUALS(ML_metric_bound, "none" ) || + CCTK_EQUALS(ML_metric_bound, "static") || + CCTK_EQUALS(ML_metric_bound, "flat" ) || + CCTK_EQUALS(ML_metric_bound, "zero" ) ) { ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, - "ML_ADM::ml_metric", ml_metric_bound); + "ML_ADM::ML_metric", ML_metric_bound); if (ierr < 0) - CCTK_WARN(0, "Failed to register ml_metric_bound BC for ML_ADM::ml_metric!"); + CCTK_WARN(0, "Failed to register ML_metric_bound BC for ML_ADM::ML_metric!"); } - if (CCTK_EQUALS(ml_shift_bound, "none" ) || - CCTK_EQUALS(ml_shift_bound, "static") || - CCTK_EQUALS(ml_shift_bound, "flat" ) || - CCTK_EQUALS(ml_shift_bound, "zero" ) ) + if (CCTK_EQUALS(ML_shift_bound, "none" ) || + CCTK_EQUALS(ML_shift_bound, "static") || + CCTK_EQUALS(ML_shift_bound, "flat" ) || + CCTK_EQUALS(ML_shift_bound, "zero" ) ) { ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, -1, - "ML_ADM::ml_shift", ml_shift_bound); + "ML_ADM::ML_shift", ML_shift_bound); if (ierr < 0) - CCTK_WARN(0, "Failed to register ml_shift_bound BC for ML_ADM::ml_shift!"); + CCTK_WARN(0, "Failed to register ML_shift_bound BC for ML_ADM::ML_shift!"); } if (CCTK_EQUALS(K11_bound, "none" ) || @@ -255,79 +255,79 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS) CCTK_WARN(0, "Failed to register beta3_bound BC for ML_ADM::beta3!"); } - if (CCTK_EQUALS(ml_curv_bound, "radiative")) + if (CCTK_EQUALS(ML_curv_bound, "radiative")) { /* apply radiation boundary condition */ - static CCTK_INT handle_ml_curv_bound = -1; - if (handle_ml_curv_bound < 0) handle_ml_curv_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_curv_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_curv_bound , ml_curv_bound_limit, "LIMIT") < 0) + static CCTK_INT handle_ML_curv_bound = -1; + if (handle_ML_curv_bound < 0) handle_ML_curv_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_curv_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_curv_bound , ML_curv_bound_limit, "LIMIT") < 0) CCTK_WARN(0, "could not set LIMIT value in table!"); - if (Util_TableSetReal(handle_ml_curv_bound ,ml_curv_bound_speed, "SPEED") < 0) + if (Util_TableSetReal(handle_ML_curv_bound ,ML_curv_bound_speed, "SPEED") < 0) CCTK_WARN(0, "could not set SPEED value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_curv_bound, - "ML_ADM::ml_curv", "Radiation"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_curv_bound, + "ML_ADM::ML_curv", "Radiation"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ml_curv!"); + CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_curv!"); } - if (CCTK_EQUALS(ml_lapse_bound, "radiative")) + if (CCTK_EQUALS(ML_lapse_bound, "radiative")) { /* apply radiation boundary condition */ - static CCTK_INT handle_ml_lapse_bound = -1; - if (handle_ml_lapse_bound < 0) handle_ml_lapse_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_lapse_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_lapse_bound , ml_lapse_bound_limit, "LIMIT") < 0) + static CCTK_INT handle_ML_lapse_bound = -1; + if (handle_ML_lapse_bound < 0) handle_ML_lapse_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_lapse_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_lapse_bound , ML_lapse_bound_limit, "LIMIT") < 0) CCTK_WARN(0, "could not set LIMIT value in table!"); - if (Util_TableSetReal(handle_ml_lapse_bound ,ml_lapse_bound_speed, "SPEED") < 0) + if (Util_TableSetReal(handle_ML_lapse_bound ,ML_lapse_bound_speed, "SPEED") < 0) CCTK_WARN(0, "could not set SPEED value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_lapse_bound, - "ML_ADM::ml_lapse", "Radiation"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_lapse_bound, + "ML_ADM::ML_lapse", "Radiation"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ml_lapse!"); + CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_lapse!"); } - if (CCTK_EQUALS(ml_metric_bound, "radiative")) + if (CCTK_EQUALS(ML_metric_bound, "radiative")) { /* apply radiation boundary condition */ - static CCTK_INT handle_ml_metric_bound = -1; - if (handle_ml_metric_bound < 0) handle_ml_metric_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_metric_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_metric_bound , ml_metric_bound_limit, "LIMIT") < 0) + static CCTK_INT handle_ML_metric_bound = -1; + if (handle_ML_metric_bound < 0) handle_ML_metric_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_metric_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_metric_bound , ML_metric_bound_limit, "LIMIT") < 0) CCTK_WARN(0, "could not set LIMIT value in table!"); - if (Util_TableSetReal(handle_ml_metric_bound ,ml_metric_bound_speed, "SPEED") < 0) + if (Util_TableSetReal(handle_ML_metric_bound ,ML_metric_bound_speed, "SPEED") < 0) CCTK_WARN(0, "could not set SPEED value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_metric_bound, - "ML_ADM::ml_metric", "Radiation"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_metric_bound, + "ML_ADM::ML_metric", "Radiation"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ml_metric!"); + CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_metric!"); } - if (CCTK_EQUALS(ml_shift_bound, "radiative")) + if (CCTK_EQUALS(ML_shift_bound, "radiative")) { /* apply radiation boundary condition */ - static CCTK_INT handle_ml_shift_bound = -1; - if (handle_ml_shift_bound < 0) handle_ml_shift_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_shift_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_shift_bound , ml_shift_bound_limit, "LIMIT") < 0) + static CCTK_INT handle_ML_shift_bound = -1; + if (handle_ML_shift_bound < 0) handle_ML_shift_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_shift_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_shift_bound , ML_shift_bound_limit, "LIMIT") < 0) CCTK_WARN(0, "could not set LIMIT value in table!"); - if (Util_TableSetReal(handle_ml_shift_bound ,ml_shift_bound_speed, "SPEED") < 0) + if (Util_TableSetReal(handle_ML_shift_bound ,ML_shift_bound_speed, "SPEED") < 0) CCTK_WARN(0, "could not set SPEED value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_shift_bound, - "ML_ADM::ml_shift", "Radiation"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_shift_bound, + "ML_ADM::ML_shift", "Radiation"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ml_shift!"); + CCTK_WARN(0, "Failed to register Radiation BC for ML_ADM::ML_shift!"); } @@ -635,71 +635,71 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS) } - if (CCTK_EQUALS(ml_curv_bound, "scalar")) + if (CCTK_EQUALS(ML_curv_bound, "scalar")) { /* apply scalar boundary condition */ - static CCTK_INT handle_ml_curv_bound = -1; - if (handle_ml_curv_bound < 0) handle_ml_curv_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_curv_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_curv_bound ,ml_curv_bound_scalar, "SCALAR") < 0) + static CCTK_INT handle_ML_curv_bound = -1; + if (handle_ML_curv_bound < 0) handle_ML_curv_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_curv_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_curv_bound ,ML_curv_bound_scalar, "SCALAR") < 0) CCTK_WARN(0, "could not set SCALAR value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_curv_bound, - "ML_ADM::ml_curv", "scalar"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_curv_bound, + "ML_ADM::ML_curv", "scalar"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ml_curv!"); + CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_curv!"); } - if (CCTK_EQUALS(ml_lapse_bound, "scalar")) + if (CCTK_EQUALS(ML_lapse_bound, "scalar")) { /* apply scalar boundary condition */ - static CCTK_INT handle_ml_lapse_bound = -1; - if (handle_ml_lapse_bound < 0) handle_ml_lapse_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_lapse_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_lapse_bound ,ml_lapse_bound_scalar, "SCALAR") < 0) + static CCTK_INT handle_ML_lapse_bound = -1; + if (handle_ML_lapse_bound < 0) handle_ML_lapse_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_lapse_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_lapse_bound ,ML_lapse_bound_scalar, "SCALAR") < 0) CCTK_WARN(0, "could not set SCALAR value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_lapse_bound, - "ML_ADM::ml_lapse", "scalar"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_lapse_bound, + "ML_ADM::ML_lapse", "scalar"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ml_lapse!"); + CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_lapse!"); } - if (CCTK_EQUALS(ml_metric_bound, "scalar")) + if (CCTK_EQUALS(ML_metric_bound, "scalar")) { /* apply scalar boundary condition */ - static CCTK_INT handle_ml_metric_bound = -1; - if (handle_ml_metric_bound < 0) handle_ml_metric_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_metric_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_metric_bound ,ml_metric_bound_scalar, "SCALAR") < 0) + static CCTK_INT handle_ML_metric_bound = -1; + if (handle_ML_metric_bound < 0) handle_ML_metric_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_metric_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_metric_bound ,ML_metric_bound_scalar, "SCALAR") < 0) CCTK_WARN(0, "could not set SCALAR value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_metric_bound, - "ML_ADM::ml_metric", "scalar"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_metric_bound, + "ML_ADM::ML_metric", "scalar"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ml_metric!"); + CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_metric!"); } - if (CCTK_EQUALS(ml_shift_bound, "scalar")) + if (CCTK_EQUALS(ML_shift_bound, "scalar")) { /* apply scalar boundary condition */ - static CCTK_INT handle_ml_shift_bound = -1; - if (handle_ml_shift_bound < 0) handle_ml_shift_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); - if (handle_ml_shift_bound < 0) CCTK_WARN(0, "could not create table!"); - if (Util_TableSetReal(handle_ml_shift_bound ,ml_shift_bound_scalar, "SCALAR") < 0) + static CCTK_INT handle_ML_shift_bound = -1; + if (handle_ML_shift_bound < 0) handle_ML_shift_bound = Util_TableCreate(UTIL_TABLE_FLAGS_CASE_INSENSITIVE); + if (handle_ML_shift_bound < 0) CCTK_WARN(0, "could not create table!"); + if (Util_TableSetReal(handle_ML_shift_bound ,ML_shift_bound_scalar, "SCALAR") < 0) CCTK_WARN(0, "could not set SCALAR value in table!"); - ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ml_shift_bound, - "ML_ADM::ml_shift", "scalar"); + ierr = Boundary_SelectGroupForBC(cctkGH, CCTK_ALL_FACES, 1, handle_ML_shift_bound, + "ML_ADM::ML_shift", "scalar"); if (ierr < 0) - CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ml_shift!"); + CCTK_WARN(0, "Failed to register Scalar BC for ML_ADM::ML_shift!"); } @@ -980,25 +980,25 @@ void ML_ADM_ApplyBoundConds(CCTK_ARGUMENTS) /* template for entries in parameter file: -#$bound$#ML_ADM::ml_curv_bound = "skip" -#$bound$#ML_ADM::ml_curv_bound_speed = 1.0 -#$bound$#ML_ADM::ml_curv_bound_limit = 0.0 -#$bound$#ML_ADM::ml_curv_bound_scalar = 0.0 +#$bound$#ML_ADM::ML_curv_bound = "skip" +#$bound$#ML_ADM::ML_curv_bound_speed = 1.0 +#$bound$#ML_ADM::ML_curv_bound_limit = 0.0 +#$bound$#ML_ADM::ML_curv_bound_scalar = 0.0 -#$bound$#ML_ADM::ml_lapse_bound = "skip" -#$bound$#ML_ADM::ml_lapse_bound_speed = 1.0 -#$bound$#ML_ADM::ml_lapse_bound_limit = 0.0 -#$bound$#ML_ADM::ml_lapse_bound_scalar = 0.0 +#$bound$#ML_ADM::ML_lapse_bound = "skip" +#$bound$#ML_ADM::ML_lapse_bound_speed = 1.0 +#$bound$#ML_ADM::ML_lapse_bound_limit = 0.0 +#$bound$#ML_ADM::ML_lapse_bound_scalar = 0.0 -#$bound$#ML_ADM::ml_metric_bound = "skip" -#$bound$#ML_ADM::ml_metric_bound_speed = 1.0 -#$bound$#ML_ADM::ml_metric_bound_limit = 0.0 -#$bound$#ML_ADM::ml_metric_bound_scalar = 0.0 +#$bound$#ML_ADM::ML_metric_bound = "skip" +#$bound$#ML_ADM::ML_metric_bound_speed = 1.0 +#$bound$#ML_ADM::ML_metric_bound_limit = 0.0 +#$bound$#ML_ADM::ML_metric_bound_scalar = 0.0 -#$bound$#ML_ADM::ml_shift_bound = "skip" -#$bound$#ML_ADM::ml_shift_bound_speed = 1.0 -#$bound$#ML_ADM::ml_shift_bound_limit = 0.0 -#$bound$#ML_ADM::ml_shift_bound_scalar = 0.0 +#$bound$#ML_ADM::ML_shift_bound = "skip" +#$bound$#ML_ADM::ML_shift_bound_speed = 1.0 +#$bound$#ML_ADM::ML_shift_bound_limit = 0.0 +#$bound$#ML_ADM::ML_shift_bound_scalar = 0.0 #$bound$#ML_ADM::K11_bound = "skip" #$bound$#ML_ADM::K11_bound_speed = 1.0 diff --git a/ML_ADM/src/ML_ADM_Minkowski.c b/ML_ADM/src/ML_ADM_Minkowski.c index 9af8b2e..62e14e2 100644 --- a/ML_ADM/src/ML_ADM_Minkowski.c +++ b/ML_ADM/src/ML_ADM_Minkowski.c @@ -5,7 +5,10 @@ #define KRANC_C +#include #include +#include +#include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -84,7 +87,7 @@ void ML_ADM_Minkowski_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL n pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ - _Pragma ("omp parallel") + #pragma omp parallel LC_LOOP3 (ML_ADM_Minkowski, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) diff --git a/ML_ADM/src/ML_ADM_constraints_boundary.c b/ML_ADM/src/ML_ADM_constraints_boundary.c index ae07d3b..58854a8 100644 --- a/ML_ADM/src/ML_ADM_constraints_boundary.c +++ b/ML_ADM/src/ML_ADM_constraints_boundary.c @@ -5,7 +5,10 @@ #define KRANC_C +#include #include +#include +#include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -84,7 +87,7 @@ void ML_ADM_constraints_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ - _Pragma ("omp parallel") + #pragma omp parallel LC_LOOP3 (ML_ADM_constraints_boundary, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) diff --git a/ML_ADM/src/ML_ADM_convertFromADMBase.c b/ML_ADM/src/ML_ADM_convertFromADMBase.c index f564d16..82493f6 100644 --- a/ML_ADM/src/ML_ADM_convertFromADMBase.c +++ b/ML_ADM/src/ML_ADM_convertFromADMBase.c @@ -5,7 +5,10 @@ #define KRANC_C +#include #include +#include +#include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -84,7 +87,7 @@ void ML_ADM_convertFromADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CC pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ - _Pragma ("omp parallel") + #pragma omp parallel LC_LOOP3 (ML_ADM_convertFromADMBase, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) diff --git a/ML_ADM/src/ML_ADM_convertToADMBase.c b/ML_ADM/src/ML_ADM_convertToADMBase.c index 277235b..4b015f1 100644 --- a/ML_ADM/src/ML_ADM_convertToADMBase.c +++ b/ML_ADM/src/ML_ADM_convertToADMBase.c @@ -5,7 +5,10 @@ #define KRANC_C +#include #include +#include +#include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -84,7 +87,7 @@ void ML_ADM_convertToADMBase_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ - _Pragma ("omp parallel") + #pragma omp parallel LC_LOOP3 (ML_ADM_convertToADMBase, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) -- cgit v1.2.3 From 5ac1c6cdfbb3df3885bc9c1de00b5726c13fbd3b Mon Sep 17 00:00:00 2001 From: Erik Schnetter Date: Mon, 5 May 2008 17:04:39 -0500 Subject: Update generated ADM thorn --- ML_ADM/src/ML_ADM_RHS.c | 408 +++++++++++++++++++++++++--------------- ML_ADM/src/ML_ADM_boundary.c | 5 +- ML_ADM/src/ML_ADM_constraints.c | 352 ++++++++++++++++++++++------------ 3 files changed, 491 insertions(+), 274 deletions(-) (limited to 'ML_ADM') diff --git a/ML_ADM/src/ML_ADM_RHS.c b/ML_ADM/src/ML_ADM_RHS.c index dac53f7..cfe4ff1 100644 --- a/ML_ADM/src/ML_ADM_RHS.c +++ b/ML_ADM/src/ML_ADM_RHS.c @@ -5,7 +5,10 @@ #define KRANC_C +#include #include +#include +#include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -84,7 +87,7 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ - _Pragma ("omp parallel") + #pragma omp parallel LC_LOOP3 (ML_ADM_RHS, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) @@ -99,7 +102,8 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE; CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE; CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE; + CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE; + CCTK_REAL gu31 = INITVALUE, gu32 = 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, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE; @@ -136,31 +140,30 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ CCTK_REAL PDstandardNth1g11 = INITVALUE; CCTK_REAL PDstandardNth2g11 = INITVALUE; CCTK_REAL PDstandardNth3g11 = INITVALUE; + CCTK_REAL PDstandardNth11g11 = INITVALUE; CCTK_REAL PDstandardNth22g11 = INITVALUE; CCTK_REAL PDstandardNth33g11 = INITVALUE; CCTK_REAL PDstandardNth12g11 = INITVALUE; CCTK_REAL PDstandardNth13g11 = INITVALUE; - CCTK_REAL PDstandardNth21g11 = INITVALUE; CCTK_REAL PDstandardNth23g11 = INITVALUE; - CCTK_REAL PDstandardNth31g11 = INITVALUE; - CCTK_REAL PDstandardNth32g11 = INITVALUE; CCTK_REAL PDstandardNth1g12 = INITVALUE; CCTK_REAL PDstandardNth2g12 = INITVALUE; CCTK_REAL PDstandardNth3g12 = INITVALUE; + CCTK_REAL PDstandardNth11g12 = INITVALUE; + CCTK_REAL PDstandardNth22g12 = INITVALUE; CCTK_REAL PDstandardNth33g12 = INITVALUE; CCTK_REAL PDstandardNth12g12 = INITVALUE; CCTK_REAL PDstandardNth13g12 = INITVALUE; CCTK_REAL PDstandardNth21g12 = INITVALUE; CCTK_REAL PDstandardNth23g12 = INITVALUE; - CCTK_REAL PDstandardNth31g12 = INITVALUE; - CCTK_REAL PDstandardNth32g12 = INITVALUE; CCTK_REAL PDstandardNth1g13 = INITVALUE; CCTK_REAL PDstandardNth2g13 = INITVALUE; CCTK_REAL PDstandardNth3g13 = INITVALUE; + CCTK_REAL PDstandardNth11g13 = INITVALUE; CCTK_REAL PDstandardNth22g13 = INITVALUE; + CCTK_REAL PDstandardNth33g13 = INITVALUE; CCTK_REAL PDstandardNth12g13 = INITVALUE; CCTK_REAL PDstandardNth13g13 = INITVALUE; - CCTK_REAL PDstandardNth21g13 = INITVALUE; CCTK_REAL PDstandardNth23g13 = INITVALUE; CCTK_REAL PDstandardNth31g13 = INITVALUE; CCTK_REAL PDstandardNth32g13 = INITVALUE; @@ -168,17 +171,18 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ CCTK_REAL PDstandardNth2g22 = INITVALUE; CCTK_REAL PDstandardNth3g22 = INITVALUE; CCTK_REAL PDstandardNth11g22 = INITVALUE; + CCTK_REAL PDstandardNth22g22 = INITVALUE; CCTK_REAL PDstandardNth33g22 = INITVALUE; CCTK_REAL PDstandardNth12g22 = INITVALUE; CCTK_REAL PDstandardNth13g22 = INITVALUE; CCTK_REAL PDstandardNth21g22 = INITVALUE; CCTK_REAL PDstandardNth23g22 = INITVALUE; - CCTK_REAL PDstandardNth31g22 = INITVALUE; - CCTK_REAL PDstandardNth32g22 = INITVALUE; CCTK_REAL PDstandardNth1g23 = INITVALUE; CCTK_REAL PDstandardNth2g23 = INITVALUE; CCTK_REAL PDstandardNth3g23 = INITVALUE; CCTK_REAL PDstandardNth11g23 = INITVALUE; + CCTK_REAL PDstandardNth22g23 = INITVALUE; + CCTK_REAL PDstandardNth33g23 = INITVALUE; CCTK_REAL PDstandardNth12g23 = INITVALUE; CCTK_REAL PDstandardNth13g23 = INITVALUE; CCTK_REAL PDstandardNth21g23 = INITVALUE; @@ -190,29 +194,20 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ CCTK_REAL PDstandardNth3g33 = INITVALUE; CCTK_REAL PDstandardNth11g33 = INITVALUE; CCTK_REAL PDstandardNth22g33 = INITVALUE; + CCTK_REAL PDstandardNth33g33 = INITVALUE; CCTK_REAL PDstandardNth12g33 = INITVALUE; CCTK_REAL PDstandardNth13g33 = INITVALUE; - CCTK_REAL PDstandardNth21g33 = INITVALUE; CCTK_REAL PDstandardNth23g33 = INITVALUE; CCTK_REAL PDstandardNth31g33 = INITVALUE; CCTK_REAL PDstandardNth32g33 = INITVALUE; CCTK_REAL PDstandardNth1K11 = INITVALUE; - CCTK_REAL PDstandardNth2K11 = INITVALUE; - CCTK_REAL PDstandardNth3K11 = INITVALUE; CCTK_REAL PDstandardNth1K12 = INITVALUE; CCTK_REAL PDstandardNth2K12 = INITVALUE; - CCTK_REAL PDstandardNth3K12 = INITVALUE; CCTK_REAL PDstandardNth1K13 = INITVALUE; - CCTK_REAL PDstandardNth2K13 = INITVALUE; CCTK_REAL PDstandardNth3K13 = INITVALUE; - CCTK_REAL PDstandardNth1K22 = INITVALUE; CCTK_REAL PDstandardNth2K22 = INITVALUE; - CCTK_REAL PDstandardNth3K22 = INITVALUE; - CCTK_REAL PDstandardNth1K23 = INITVALUE; CCTK_REAL PDstandardNth2K23 = INITVALUE; CCTK_REAL PDstandardNth3K23 = INITVALUE; - CCTK_REAL PDstandardNth1K33 = INITVALUE; - CCTK_REAL PDstandardNth2K33 = INITVALUE; CCTK_REAL PDstandardNth3K33 = INITVALUE; /* Assign local copies of grid functions */ @@ -259,12 +254,17 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ PDstandardNth1g11 = PDstandardNth1(g11, i, j, k); PDstandardNth2g11 = PDstandardNth2(g11, i, j, k); PDstandardNth3g11 = PDstandardNth3(g11, i, j, k); + PDstandardNth11g11 = PDstandardNth11(g11, i, j, k); PDstandardNth22g11 = PDstandardNth22(g11, i, j, k); PDstandardNth33g11 = PDstandardNth33(g11, i, j, k); + PDstandardNth12g11 = PDstandardNth12(g11, i, j, k); + PDstandardNth13g11 = PDstandardNth13(g11, i, j, k); PDstandardNth23g11 = PDstandardNth23(g11, i, j, k); PDstandardNth1g12 = PDstandardNth1(g12, i, j, k); PDstandardNth2g12 = PDstandardNth2(g12, i, j, k); PDstandardNth3g12 = PDstandardNth3(g12, i, j, k); + PDstandardNth11g12 = PDstandardNth11(g12, i, j, k); + PDstandardNth22g12 = PDstandardNth22(g12, i, j, k); PDstandardNth33g12 = PDstandardNth33(g12, i, j, k); PDstandardNth12g12 = PDstandardNth12(g12, i, j, k); PDstandardNth13g12 = PDstandardNth13(g12, i, j, k); @@ -272,7 +272,9 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ PDstandardNth1g13 = PDstandardNth1(g13, i, j, k); PDstandardNth2g13 = PDstandardNth2(g13, i, j, k); PDstandardNth3g13 = PDstandardNth3(g13, i, j, k); + PDstandardNth11g13 = PDstandardNth11(g13, i, j, k); PDstandardNth22g13 = PDstandardNth22(g13, i, j, k); + PDstandardNth33g13 = PDstandardNth33(g13, i, j, k); PDstandardNth12g13 = PDstandardNth12(g13, i, j, k); PDstandardNth13g13 = PDstandardNth13(g13, i, j, k); PDstandardNth23g13 = PDstandardNth23(g13, i, j, k); @@ -280,12 +282,17 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ PDstandardNth2g22 = PDstandardNth2(g22, i, j, k); PDstandardNth3g22 = PDstandardNth3(g22, i, j, k); PDstandardNth11g22 = PDstandardNth11(g22, i, j, k); + PDstandardNth22g22 = PDstandardNth22(g22, i, j, k); PDstandardNth33g22 = PDstandardNth33(g22, i, j, k); + PDstandardNth12g22 = PDstandardNth12(g22, i, j, k); PDstandardNth13g22 = PDstandardNth13(g22, i, j, k); + PDstandardNth23g22 = PDstandardNth23(g22, i, j, k); PDstandardNth1g23 = PDstandardNth1(g23, i, j, k); PDstandardNth2g23 = PDstandardNth2(g23, i, j, k); PDstandardNth3g23 = PDstandardNth3(g23, i, j, k); PDstandardNth11g23 = PDstandardNth11(g23, i, j, k); + PDstandardNth22g23 = PDstandardNth22(g23, i, j, k); + PDstandardNth33g23 = PDstandardNth33(g23, i, j, k); PDstandardNth12g23 = PDstandardNth12(g23, i, j, k); PDstandardNth13g23 = PDstandardNth13(g23, i, j, k); PDstandardNth23g23 = PDstandardNth23(g23, i, j, k); @@ -294,206 +301,297 @@ void ML_ADM_RHS_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL normal[ PDstandardNth3g33 = PDstandardNth3(g33, i, j, k); PDstandardNth11g33 = PDstandardNth11(g33, i, j, k); PDstandardNth22g33 = PDstandardNth22(g33, i, j, k); + PDstandardNth33g33 = PDstandardNth33(g33, i, j, k); PDstandardNth12g33 = PDstandardNth12(g33, i, j, k); + PDstandardNth13g33 = PDstandardNth13(g33, i, j, k); + PDstandardNth23g33 = PDstandardNth23(g33, i, j, k); PDstandardNth1K11 = PDstandardNth1(K11, i, j, k); - PDstandardNth2K11 = PDstandardNth2(K11, i, j, k); - PDstandardNth3K11 = PDstandardNth3(K11, i, j, k); PDstandardNth1K12 = PDstandardNth1(K12, i, j, k); PDstandardNth2K12 = PDstandardNth2(K12, i, j, k); - PDstandardNth3K12 = PDstandardNth3(K12, i, j, k); PDstandardNth1K13 = PDstandardNth1(K13, i, j, k); - PDstandardNth2K13 = PDstandardNth2(K13, i, j, k); PDstandardNth3K13 = PDstandardNth3(K13, i, j, k); - PDstandardNth1K22 = PDstandardNth1(K22, i, j, k); PDstandardNth2K22 = PDstandardNth2(K22, i, j, k); - PDstandardNth3K22 = PDstandardNth3(K22, i, j, k); - PDstandardNth1K23 = PDstandardNth1(K23, i, j, k); PDstandardNth2K23 = PDstandardNth2(K23, i, j, k); PDstandardNth3K23 = PDstandardNth3(K23, i, j, k); - PDstandardNth1K33 = PDstandardNth1(K33, i, j, k); - PDstandardNth2K33 = PDstandardNth2(K33, i, j, k); PDstandardNth3K33 = PDstandardNth3(K33, i, j, k); /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L); + detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - + g11L*SQR(g23L); gu11 = INV(detg)*(g22L*g33L - SQR(g23L)); - gu21 = (g13L*g23L - g12L*g33L)*INV(detg); + gu12 = (g13L*g23L - g12L*g33L)*INV(detg); - gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + + gu21 = (g13L*g23L - g12L*g33L)*INV(detg); gu22 = INV(detg)*(g11L*g33L - SQR(g13L)); + gu23 = (g12L*g13L - g11L*g23L)*INV(detg); + + gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + gu32 = (g12L*g13L - g11L*g23L)*INV(detg); gu33 = INV(detg)*(g11L*g22L - SQR(g12L)); - G111 = khalf*(gu11*PDstandardNth1g11 + 2*(gu21*PDstandardNth1g12 + gu31*PDstandardNth1g13) - gu21*PDstandardNth2g11 - - gu31*PDstandardNth3g11); + G111 = khalf*(gu11*PDstandardNth1g11 + + 2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - + gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11); - G211 = khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu32*PDstandardNth1g13) - gu22*PDstandardNth2g11 - - gu32*PDstandardNth3g11); + G211 = khalf*(gu21*PDstandardNth1g11 + + 2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - + gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11); - G311 = khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - - gu33*PDstandardNth3g11); + G311 = khalf*(gu31*PDstandardNth1g11 + + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - + gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11); - G112 = khalf*(gu21*PDstandardNth1g22 + gu11*PDstandardNth2g11 + - gu31*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); + G112 = khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + + gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); G212 = khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + - gu32*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); + gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); G312 = khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G113 = khalf*(gu31*PDstandardNth1g33 + gu11*PDstandardNth3g11 + - gu21*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); + G113 = khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + + gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G213 = khalf*(gu32*PDstandardNth1g33 + gu21*PDstandardNth3g11 + + G213 = khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); G313 = khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu21*PDstandardNth2g22 + - gu31*(2*PDstandardNth2g23 - PDstandardNth3g22)); + G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + + gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + - gu32*(2*PDstandardNth2g23 - PDstandardNth3g22)); + G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + + gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + - gu33*(2*PDstandardNth2g23 - PDstandardNth3g22)); + G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + + gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G123 = khalf*(gu31*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + - gu21*PDstandardNth3g22); + G123 = khalf*(gu13*PDstandardNth2g33 + + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + gu12*PDstandardNth3g22); - G223 = khalf*(gu32*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + G223 = khalf*(gu23*PDstandardNth2g33 + + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu22*PDstandardNth3g22); - G323 = khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + G323 = khalf*(gu33*PDstandardNth2g33 + + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu32*PDstandardNth3g22); - G133 = khalf*(-(gu11*PDstandardNth1g33) - gu21*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + - 2*gu21*PDstandardNth3g23 + gu31*PDstandardNth3g33); - - G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + - 2*gu22*PDstandardNth3g23 + gu32*PDstandardNth3g33); - - G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + - 2*gu32*PDstandardNth3g23 + gu33*PDstandardNth3g33); - - R11 = khalf*(4*G213*G312 + G211*(2*G112 - 2*G222 - 2*G323) + G311*(2*G113 - 2*G333) - gu22*PDstandardNth11g22 - - 2*(G111*G212 + G223*G311 + G111*G313 + gu32*PDstandardNth11g23) - gu33*PDstandardNth11g33 + - 2*gu22*PDstandardNth12g12 + 2*gu32*PDstandardNth12g13 + 2*gu32*PDstandardNth13g12 + 2*gu33*PDstandardNth13g13 - - gu22*PDstandardNth22g11 - 2*gu32*PDstandardNth23g11 - gu33*PDstandardNth33g11 + 2*SQR(G212) + 2*SQR(G313)); - - R12 = khalf*(2*(G122*G211 + G123*G311 + G213*G322 + G313*G323) - - 2*(G112*G212 + G112*G313 + G212*G323 + G312*G333 + gu21*PDstandardNth12g12) - gu32*PDstandardNth12g23 - - gu33*PDstandardNth12g33 + gu31*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12) + - gu32*PDstandardNth13g22 + gu33*PDstandardNth13g23 + gu21*(PDstandardNth11g22 + PDstandardNth22g11) + - gu32*PDstandardNth22g13 + gu31*PDstandardNth23g11 - gu32*PDstandardNth23g12 + gu33*PDstandardNth23g13 - - gu33*PDstandardNth33g12); - - R13 = khalf*(2*(G123*G211 + G212*G223 + G133*G311 + G233*G312) - - 2*(G213*G222 + G223*G313 + G113*(G212 + G313) + gu31*PDstandardNth13g13) + - gu21*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12 + PDstandardNth23g11) + - gu22*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + - gu31*(PDstandardNth11g33 + PDstandardNth33g11) + - gu32*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12)); - - R22 = khalf*(4*G123*G312 + G122*(2*G212 - 2*(G111 + G313)) + G322*(2*G223 - 2*G333) - - 2*(G113*G322 + G222*(G112 + G323) + gu31*PDstandardNth13g22) + - gu11*(-PDstandardNth11g22 + 2*PDstandardNth12g12 - PDstandardNth22g11) + - gu31*(-2*PDstandardNth22g13 + 2*(PDstandardNth12g23 + PDstandardNth23g12)) + - gu33*(-PDstandardNth22g33 + 2*PDstandardNth23g23 - PDstandardNth33g22) + 2*(SQR(G112) + SQR(G323))); - - R23 = khalf*(2*(G112*G113 + G122*G213 + G133*G312 + G233*G322) + - gu11*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + - gu21*(-PDstandardNth12g23 + PDstandardNth13g22 + PDstandardNth22g13 - PDstandardNth23g12) - - 2*(G111*G123 + G113*G323 + G223*(G112 + G323) + gu32*PDstandardNth23g23) + - gu31*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12) + - gu32*(PDstandardNth22g33 + PDstandardNth33g22)); - - R33 = khalf*(4*G123*G213 - gu11*PDstandardNth11g33 - - 2*(G111*G133 + G133*G212 + G112*G233 + G222*G233 + G113*G333 + G223*G333 + gu21*PDstandardNth12g33) + - 2*(G133*G313 + G233*G323 + gu11*PDstandardNth13g13) + 2*gu21*PDstandardNth13g23 - gu22*PDstandardNth22g33 + - 2*gu21*PDstandardNth23g13 + 2*gu22*PDstandardNth23g23 - gu11*PDstandardNth33g11 - 2*gu21*PDstandardNth33g12 - - gu22*PDstandardNth33g22 + 2*SQR(G113) + 2*SQR(G223)); - - Km11 = gu11*K11L + gu21*K12L + gu31*K13L; - - Km21 = gu21*K11L + gu22*K12L + gu32*K13L; + G133 = khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + + 2*gu11*PDstandardNth3g13 + 2*gu12*PDstandardNth3g23 + + gu13*PDstandardNth3g33); + + G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + + 2*gu21*PDstandardNth3g13 + 2*gu22*PDstandardNth3g23 + + gu23*PDstandardNth3g33); + + G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + + 2*gu31*PDstandardNth3g13 + 2*gu32*PDstandardNth3g23 + + gu33*PDstandardNth3g33); + + R11 = 2*(G112*G211 + G113*G311 + G213*G312) - G111*(G111 + G212 + G313) - + G211*(G112 + G222 + G323) - G311*(G113 + G223 + G333) + + khalf*(-6*gu11*PDstandardNth11g11 + + gu12*(-8*PDstandardNth11g12 + 2*PDstandardNth12g12) + + gu22*(-9*PDstandardNth11g22 + 3*PDstandardNth12g22) + + gu13*(-8*PDstandardNth11g13 + 2*PDstandardNth13g13) + + gu33*(-9*PDstandardNth11g33 + 3*PDstandardNth13g33) + + gu21*(-6*PDstandardNth11g12 - PDstandardNth12g12 + PDstandardNth22g12) + + gu23*(-9*PDstandardNth11g23 + 3*PDstandardNth13g23 + + PDstandardNth22g23 - PDstandardNth23g23) + + gu31*(-6*PDstandardNth11g13 - PDstandardNth13g13 + PDstandardNth33g13) + + gu32*(-9*PDstandardNth11g23 + 3*PDstandardNth12g23 - + PDstandardNth23g23 + PDstandardNth33g23)) + SQR(G111) + SQR(G212) + + SQR(G313); + + R12 = khalf*(-2*(G112*G212 + G112*G313 + G212*G323 + G312*G333) - + 6*(gu11*PDstandardNth12g11 + gu31*PDstandardNth12g13) - + 9*(gu22*PDstandardNth12g22 + gu23*PDstandardNth12g23 + + gu32*PDstandardNth12g23 + gu33*PDstandardNth12g33) - + gu31*PDstandardNth13g13 + gu21* + (-7*PDstandardNth12g12 + PDstandardNth22g12) + + gu12*(PDstandardNth11g12 - 10*PDstandardNth12g12 + + 3*PDstandardNth22g12) + 3*gu22*PDstandardNth22g22 + + gu23*PDstandardNth22g23 + 3*gu32*PDstandardNth22g23 + + gu13*(PDstandardNth11g13 - 9*PDstandardNth12g13 - PDstandardNth13g13 + + 3*PDstandardNth23g13) - gu32*PDstandardNth23g23 + + 2*(G122*G211 + G123*G311 + G213*G322 + G313*G323 + + gu23*PDstandardNth23g23) + 3*gu33*PDstandardNth23g33 + + gu31*PDstandardNth33g13 + gu32*PDstandardNth33g23); + + R13 = khalf*(-2*(G213*G222 + G223*G313 + G113*(G212 + G313)) - + 6*gu11*PDstandardNth13g11 - + 9*(gu22*PDstandardNth13g22 + (gu23 + gu32)*PDstandardNth13g23 + + gu33*PDstandardNth13g33) + + gu21*(-PDstandardNth12g12 - 6*PDstandardNth13g12 + PDstandardNth22g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 9*PDstandardNth13g12 + + 3*PDstandardNth23g12) + + gu23*(PDstandardNth22g23 - PDstandardNth23g23) + + 2*(G123*G211 + G212*G223 + G133*G311 + G233*G312 + + gu32*PDstandardNth23g23) + + gu31*(-7*PDstandardNth13g13 + PDstandardNth33g13) + + gu13*(PDstandardNth11g13 - 10*PDstandardNth13g13 + + 3*PDstandardNth33g13) + gu32*PDstandardNth33g23 + + 3*(gu22*PDstandardNth23g22 + gu23*PDstandardNth33g23 + + gu33*PDstandardNth33g33)); + + R22 = -(G122*(G111 + G212 + G313)) + 2*(G122*G212 + G123*G312 + G223*G322) - + G222*(G112 + G222 + G323) - G322*(G113 + G223 + G333) + + khalf*(3*gu11*(PDstandardNth12g11 - 3*PDstandardNth22g11) + + gu21*(2*PDstandardNth12g12 - 8*PDstandardNth22g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth22g12) - + 6*gu22*PDstandardNth22g22 + + gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + + 3*PDstandardNth23g13) - + gu23*(8*PDstandardNth22g23 - 2*PDstandardNth23g23) + + gu33*(-9*PDstandardNth22g33 + 3*PDstandardNth23g33) + + gu31*(3*PDstandardNth12g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + + PDstandardNth33g13) + gu32* + (-6*PDstandardNth22g23 - PDstandardNth23g23 + PDstandardNth33g23)) + + SQR(G112) + SQR(G222) + SQR(G323); + + R23 = khalf*(-2*(G111*G123 + G113*G323 + G223*(G112 + G323)) + + 2*(G112*G113 + G122*G213 + G133*G312 + G233*G322 + + gu31*PDstandardNth13g13) + + gu11*(3*PDstandardNth13g11 - 9*PDstandardNth23g11) + + gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - + 9*PDstandardNth23g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth23g12) - + 6*gu22*PDstandardNth23g22 - + 9*(gu31*PDstandardNth23g13 + gu33*PDstandardNth23g33) + + gu31*PDstandardNth33g13 + gu13* + (PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth23g13 + + 3*PDstandardNth33g13) + + gu32*(-7*PDstandardNth23g23 + PDstandardNth33g23) + + gu23*(PDstandardNth22g23 - 10*PDstandardNth23g23 + + 3*PDstandardNth33g23) + 3*gu33*PDstandardNth33g33); + + R33 = -(G133*(G111 + G212 + G313)) + 2*(G123*G213 + G133*G313) + 2*G233*G323 - + G233*(G112 + G222 + G323) - G333*(G113 + G223 + G333) + + khalf*(3*gu11*(PDstandardNth13g11 - 3*PDstandardNth33g11) + + gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - + 9*PDstandardNth33g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 + 3*PDstandardNth23g12 - + 9*PDstandardNth33g12) + + gu31*(2*PDstandardNth13g13 - 8*PDstandardNth33g13) + + gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 6*PDstandardNth33g13) + + 3*gu22*(PDstandardNth23g22 - 3*PDstandardNth33g22) + + gu32*(2*PDstandardNth23g23 - 8*PDstandardNth33g23) + + gu23*(PDstandardNth22g23 - PDstandardNth23g23 - 6*PDstandardNth33g23) - + 6*gu33*PDstandardNth33g33) + SQR(G113) + SQR(G223) + SQR(G333); + + Km11 = gu11*K11L + gu12*K12L + gu13*K13L; + + Km21 = gu21*K11L + gu22*K12L + gu23*K13L; Km31 = gu31*K11L + gu32*K12L + gu33*K13L; - Km12 = gu11*K12L + gu21*K22L + gu31*K23L; + Km12 = gu11*K12L + gu12*K22L + gu13*K23L; - Km22 = gu21*K12L + gu22*K22L + gu32*K23L; + Km22 = gu21*K12L + gu22*K22L + gu23*K23L; Km32 = gu31*K12L + gu32*K22L + gu33*K23L; - Km13 = gu11*K13L + gu21*K23L + gu31*K33L; + Km13 = gu11*K13L + gu12*K23L + gu13*K33L; - Km23 = gu21*K13L + gu22*K23L + gu32*K33L; + Km23 = gu21*K13L + gu22*K23L + gu23*K33L; Km33 = gu31*K13L + gu32*K23L + gu33*K33L; trK = Km11 + Km22 + Km33; - g11rhsL = -2*alphaL*K11L + 2*(g11L*PDstandardNth1beta1 + g12L*PDstandardNth1beta2 + g13L*PDstandardNth1beta3) + - beta1L*PDstandardNth1g11 + beta2L*PDstandardNth2g11 + beta3L*PDstandardNth3g11; - - g12rhsL = -2*alphaL*K12L + g22L*PDstandardNth1beta2 + g23L*PDstandardNth1beta3 + beta1L*PDstandardNth1g12 + - g11L*PDstandardNth2beta1 + g12L*(PDstandardNth1beta1 + PDstandardNth2beta2) + g13L*PDstandardNth2beta3 + - beta2L*PDstandardNth2g12 + beta3L*PDstandardNth3g12; - - g13rhsL = -2*alphaL*K13L + g23L*PDstandardNth1beta2 + g33L*PDstandardNth1beta3 + beta1L*PDstandardNth1g13 + - beta2L*PDstandardNth2g13 + g11L*PDstandardNth3beta1 + g12L*PDstandardNth3beta2 + - g13L*(PDstandardNth1beta1 + PDstandardNth3beta3) + beta3L*PDstandardNth3g13; - - g22rhsL = -2*alphaL*K22L + beta1L*PDstandardNth1g22 + - 2*(g12L*PDstandardNth2beta1 + g22L*PDstandardNth2beta2 + g23L*PDstandardNth2beta3) + beta2L*PDstandardNth2g22 + - beta3L*PDstandardNth3g22; - - g23rhsL = -2*alphaL*K23L + beta1L*PDstandardNth1g23 + g13L*PDstandardNth2beta1 + g33L*PDstandardNth2beta3 + - beta2L*PDstandardNth2g23 + g12L*PDstandardNth3beta1 + g22L*PDstandardNth3beta2 + - g23L*(PDstandardNth2beta2 + PDstandardNth3beta3) + beta3L*PDstandardNth3g23; - - g33rhsL = -2*alphaL*K33L + beta1L*PDstandardNth1g33 + beta2L*PDstandardNth2g33 + - 2*(g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + g33L*PDstandardNth3beta3) + beta3L*PDstandardNth3g33; - - K11rhsL = -PDstandardNth11alpha + G111*PDstandardNth1alpha + - 2*(K11L*PDstandardNth1beta1 + K12L*PDstandardNth1beta2 + K13L*PDstandardNth1beta3) + beta1L*PDstandardNth1K11 + - G211*PDstandardNth2alpha + beta2L*PDstandardNth2K11 + G311*PDstandardNth3alpha + beta3L*PDstandardNth3K11 + + g11rhsL = -2*alphaL*K11L + 4*(g11L*PDstandardNth1beta1 + + g12L*PDstandardNth1beta2 + g13L*PDstandardNth1beta3) + + beta1L*PDstandardNth1g11 + beta2L*PDstandardNth2g12 + + beta3L*PDstandardNth3g13; + + g12rhsL = -2*alphaL*K12L + 3*(g12L*PDstandardNth1beta1 + + g22L*PDstandardNth1beta2 + g23L*PDstandardNth1beta3) + + beta1L*PDstandardNth1g12 + g12L*PDstandardNth2beta1 + + g22L*PDstandardNth2beta2 + g23L*PDstandardNth2beta3 + + beta2L*PDstandardNth2g22 + beta3L*PDstandardNth3g23; + + g13rhsL = -2*alphaL*K13L + 3*(g13L*PDstandardNth1beta1 + + g23L*PDstandardNth1beta2 + g33L*PDstandardNth1beta3) + + beta1L*PDstandardNth1g13 + beta2L*PDstandardNth2g23 + + g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + + g33L*PDstandardNth3beta3 + beta3L*PDstandardNth3g33; + + g22rhsL = -2*alphaL*K22L + beta1L*PDstandardNth1g12 + + 4*(g12L*PDstandardNth2beta1 + g22L*PDstandardNth2beta2 + + g23L*PDstandardNth2beta3) + beta2L*PDstandardNth2g22 + + beta3L*PDstandardNth3g23; + + g23rhsL = -2*alphaL*K23L + beta1L*PDstandardNth1g13 + + 3*(g13L*PDstandardNth2beta1 + g23L*PDstandardNth2beta2 + + g33L*PDstandardNth2beta3) + beta2L*PDstandardNth2g23 + + g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + + g33L*PDstandardNth3beta3 + beta3L*PDstandardNth3g33; + + g33rhsL = -2*alphaL*K33L + beta1L*PDstandardNth1g13 + + beta2L*PDstandardNth2g23 + 4* + (g13L*PDstandardNth3beta1 + g23L*PDstandardNth3beta2 + + g33L*PDstandardNth3beta3) + beta3L*PDstandardNth3g33; + + K11rhsL = -9*PDstandardNth11alpha + G111*PDstandardNth1alpha + + 4*(K11L*PDstandardNth1beta1 + K12L*PDstandardNth1beta2 + + K13L*PDstandardNth1beta3) + beta1L*PDstandardNth1K11 + + G212*PDstandardNth2alpha + beta2L*PDstandardNth2K12 + + G313*PDstandardNth3alpha + beta3L*PDstandardNth3K13 + alphaL*(-2*(K11L*Km11 + K12L*Km21 + K13L*Km31) + R11 + K11L*trK); - K12rhsL = -PDstandardNth12alpha + G112*PDstandardNth1alpha + K22L*PDstandardNth1beta2 + K23L*PDstandardNth1beta3 + - beta1L*PDstandardNth1K12 + G212*PDstandardNth2alpha + K11L*PDstandardNth2beta1 + - K12L*(PDstandardNth1beta1 + PDstandardNth2beta2) + K13L*PDstandardNth2beta3 + beta2L*PDstandardNth2K12 + - G312*PDstandardNth3alpha + beta3L*PDstandardNth3K12 + - alphaL*(-2*(K11L*Km12 + K12L*Km22 + K13L*Km32) + R12 + K12L*trK); - - K13rhsL = -PDstandardNth13alpha + G113*PDstandardNth1alpha + K23L*PDstandardNth1beta2 + K33L*PDstandardNth1beta3 + - beta1L*PDstandardNth1K13 + G213*PDstandardNth2alpha + beta2L*PDstandardNth2K13 + G313*PDstandardNth3alpha + - K11L*PDstandardNth3beta1 + K12L*PDstandardNth3beta2 + K13L*(PDstandardNth1beta1 + PDstandardNth3beta3) + - beta3L*PDstandardNth3K13 + alphaL*(-2*(K11L*Km13 + K12L*Km23 + K13L*Km33) + R13 + K13L*trK); - - K22rhsL = G122*PDstandardNth1alpha + beta1L*PDstandardNth1K22 - PDstandardNth22alpha + G222*PDstandardNth2alpha + - 2*(K12L*PDstandardNth2beta1 + K22L*PDstandardNth2beta2 + K23L*PDstandardNth2beta3) + beta2L*PDstandardNth2K22 + - G322*PDstandardNth3alpha + beta3L*PDstandardNth3K22 + + K12rhsL = -9*PDstandardNth12alpha + G112*PDstandardNth1alpha + + 3*(K12L*PDstandardNth1beta1 + K22L*PDstandardNth1beta2 + + K23L*PDstandardNth1beta3) + beta1L*PDstandardNth1K12 + + G222*PDstandardNth2alpha + K12L*PDstandardNth2beta1 + + K22L*PDstandardNth2beta2 + K23L*PDstandardNth2beta3 + + beta2L*PDstandardNth2K22 + G323*PDstandardNth3alpha + + beta3L*PDstandardNth3K23 + alphaL* + (-2*(K11L*Km12 + K12L*Km22 + K13L*Km32) + R12 + K12L*trK); + + K13rhsL = -9*PDstandardNth13alpha + G113*PDstandardNth1alpha + + 3*(K13L*PDstandardNth1beta1 + K23L*PDstandardNth1beta2 + + K33L*PDstandardNth1beta3) + beta1L*PDstandardNth1K13 + + G223*PDstandardNth2alpha + beta2L*PDstandardNth2K23 + + G333*PDstandardNth3alpha + K13L*PDstandardNth3beta1 + + K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3 + + beta3L*PDstandardNth3K33 + alphaL* + (-2*(K11L*Km13 + K12L*Km23 + K13L*Km33) + R13 + K13L*trK); + + K22rhsL = G112*PDstandardNth1alpha + beta1L*PDstandardNth1K12 - + 9*PDstandardNth22alpha + G222*PDstandardNth2alpha + + 4*(K12L*PDstandardNth2beta1 + K22L*PDstandardNth2beta2 + + K23L*PDstandardNth2beta3) + beta2L*PDstandardNth2K22 + + G323*PDstandardNth3alpha + beta3L*PDstandardNth3K23 + alphaL*(-2*(K12L*Km12 + K22L*Km22 + K23L*Km32) + R22 + K22L*trK); - K23rhsL = G123*PDstandardNth1alpha + beta1L*PDstandardNth1K23 - PDstandardNth23alpha + G223*PDstandardNth2alpha + - K13L*PDstandardNth2beta1 + K33L*PDstandardNth2beta3 + beta2L*PDstandardNth2K23 + G323*PDstandardNth3alpha + - K12L*PDstandardNth3beta1 + K22L*PDstandardNth3beta2 + K23L*(PDstandardNth2beta2 + PDstandardNth3beta3) + - beta3L*PDstandardNth3K23 + alphaL*(-2*(K12L*Km13 + K22L*Km23 + K23L*Km33) + R23 + K23L*trK); - - K33rhsL = G133*PDstandardNth1alpha + beta1L*PDstandardNth1K33 + G233*PDstandardNth2alpha + beta2L*PDstandardNth2K33 - - PDstandardNth33alpha + G333*PDstandardNth3alpha + - 2*(K13L*PDstandardNth3beta1 + K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3) + beta3L*PDstandardNth3K33 + + K23rhsL = G113*PDstandardNth1alpha + beta1L*PDstandardNth1K13 - + 9*PDstandardNth23alpha + G223*PDstandardNth2alpha + + 3*(K13L*PDstandardNth2beta1 + K23L*PDstandardNth2beta2 + + K33L*PDstandardNth2beta3) + beta2L*PDstandardNth2K23 + + G333*PDstandardNth3alpha + K13L*PDstandardNth3beta1 + + K23L*PDstandardNth3beta2 + K33L*PDstandardNth3beta3 + + beta3L*PDstandardNth3K33 + alphaL* + (-2*(K12L*Km13 + K22L*Km23 + K23L*Km33) + R23 + K23L*trK); + + K33rhsL = G113*PDstandardNth1alpha + beta1L*PDstandardNth1K13 + + G223*PDstandardNth2alpha + beta2L*PDstandardNth2K23 - + 9*PDstandardNth33alpha + G333*PDstandardNth3alpha + + 4*(K13L*PDstandardNth3beta1 + K23L*PDstandardNth3beta2 + + K33L*PDstandardNth3beta3) + beta3L*PDstandardNth3K33 + alphaL*(-2*(K13L*Km13 + K23L*Km23 + K33L*Km33) + R33 + K33L*trK); alpharhsL = 0; diff --git a/ML_ADM/src/ML_ADM_boundary.c b/ML_ADM/src/ML_ADM_boundary.c index f4f0289..f3fd84f 100644 --- a/ML_ADM/src/ML_ADM_boundary.c +++ b/ML_ADM/src/ML_ADM_boundary.c @@ -5,7 +5,10 @@ #define KRANC_C +#include #include +#include +#include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -84,7 +87,7 @@ void ML_ADM_boundary_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL no pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ - _Pragma ("omp parallel") + #pragma omp parallel LC_LOOP3 (ML_ADM_boundary, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) diff --git a/ML_ADM/src/ML_ADM_constraints.c b/ML_ADM/src/ML_ADM_constraints.c index 7ee74e2..a2e08a2 100644 --- a/ML_ADM/src/ML_ADM_constraints.c +++ b/ML_ADM/src/ML_ADM_constraints.c @@ -5,7 +5,10 @@ #define KRANC_C +#include #include +#include +#include #include "cctk.h" #include "cctk_Arguments.h" #include "cctk_Parameters.h" @@ -84,7 +87,7 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL pm1o12dz2 = -pow(dz,-2)/12.; /* Loop over the grid points */ - _Pragma ("omp parallel") + #pragma omp parallel LC_LOOP3 (ML_ADM_constraints, i,j,k, min[0],min[1],min[2], max[0],max[1],max[2], cctk_lsh[0],cctk_lsh[1],cctk_lsh[2]) @@ -99,7 +102,8 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL CCTK_REAL G111 = INITVALUE, G112 = INITVALUE, G113 = INITVALUE, G122 = INITVALUE, G123 = INITVALUE, G133 = INITVALUE; CCTK_REAL G211 = INITVALUE, G212 = INITVALUE, G213 = INITVALUE, G222 = INITVALUE, G223 = INITVALUE, G233 = INITVALUE; CCTK_REAL G311 = INITVALUE, G312 = INITVALUE, G313 = INITVALUE, G322 = INITVALUE, G323 = INITVALUE, G333 = INITVALUE; - CCTK_REAL gu11 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu31 = INITVALUE, gu32 = INITVALUE, gu33 = INITVALUE; + CCTK_REAL gu11 = INITVALUE, gu12 = INITVALUE, gu13 = INITVALUE, gu21 = INITVALUE, gu22 = INITVALUE, gu23 = INITVALUE; + CCTK_REAL gu31 = INITVALUE, gu32 = 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, R12 = INITVALUE, R13 = INITVALUE, R22 = INITVALUE, R23 = INITVALUE, R33 = INITVALUE; @@ -117,31 +121,30 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL CCTK_REAL PDstandardNth1g11 = INITVALUE; CCTK_REAL PDstandardNth2g11 = INITVALUE; CCTK_REAL PDstandardNth3g11 = INITVALUE; + CCTK_REAL PDstandardNth11g11 = INITVALUE; CCTK_REAL PDstandardNth22g11 = INITVALUE; CCTK_REAL PDstandardNth33g11 = INITVALUE; CCTK_REAL PDstandardNth12g11 = INITVALUE; CCTK_REAL PDstandardNth13g11 = INITVALUE; - CCTK_REAL PDstandardNth21g11 = INITVALUE; CCTK_REAL PDstandardNth23g11 = INITVALUE; - CCTK_REAL PDstandardNth31g11 = INITVALUE; - CCTK_REAL PDstandardNth32g11 = INITVALUE; CCTK_REAL PDstandardNth1g12 = INITVALUE; CCTK_REAL PDstandardNth2g12 = INITVALUE; CCTK_REAL PDstandardNth3g12 = INITVALUE; + CCTK_REAL PDstandardNth11g12 = INITVALUE; + CCTK_REAL PDstandardNth22g12 = INITVALUE; CCTK_REAL PDstandardNth33g12 = INITVALUE; CCTK_REAL PDstandardNth12g12 = INITVALUE; CCTK_REAL PDstandardNth13g12 = INITVALUE; CCTK_REAL PDstandardNth21g12 = INITVALUE; CCTK_REAL PDstandardNth23g12 = INITVALUE; - CCTK_REAL PDstandardNth31g12 = INITVALUE; - CCTK_REAL PDstandardNth32g12 = INITVALUE; CCTK_REAL PDstandardNth1g13 = INITVALUE; CCTK_REAL PDstandardNth2g13 = INITVALUE; CCTK_REAL PDstandardNth3g13 = INITVALUE; + CCTK_REAL PDstandardNth11g13 = INITVALUE; CCTK_REAL PDstandardNth22g13 = INITVALUE; + CCTK_REAL PDstandardNth33g13 = INITVALUE; CCTK_REAL PDstandardNth12g13 = INITVALUE; CCTK_REAL PDstandardNth13g13 = INITVALUE; - CCTK_REAL PDstandardNth21g13 = INITVALUE; CCTK_REAL PDstandardNth23g13 = INITVALUE; CCTK_REAL PDstandardNth31g13 = INITVALUE; CCTK_REAL PDstandardNth32g13 = INITVALUE; @@ -149,17 +152,18 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL CCTK_REAL PDstandardNth2g22 = INITVALUE; CCTK_REAL PDstandardNth3g22 = INITVALUE; CCTK_REAL PDstandardNth11g22 = INITVALUE; + CCTK_REAL PDstandardNth22g22 = INITVALUE; CCTK_REAL PDstandardNth33g22 = INITVALUE; CCTK_REAL PDstandardNth12g22 = INITVALUE; CCTK_REAL PDstandardNth13g22 = INITVALUE; CCTK_REAL PDstandardNth21g22 = INITVALUE; CCTK_REAL PDstandardNth23g22 = INITVALUE; - CCTK_REAL PDstandardNth31g22 = INITVALUE; - CCTK_REAL PDstandardNth32g22 = INITVALUE; CCTK_REAL PDstandardNth1g23 = INITVALUE; CCTK_REAL PDstandardNth2g23 = INITVALUE; CCTK_REAL PDstandardNth3g23 = INITVALUE; CCTK_REAL PDstandardNth11g23 = INITVALUE; + CCTK_REAL PDstandardNth22g23 = INITVALUE; + CCTK_REAL PDstandardNth33g23 = INITVALUE; CCTK_REAL PDstandardNth12g23 = INITVALUE; CCTK_REAL PDstandardNth13g23 = INITVALUE; CCTK_REAL PDstandardNth21g23 = INITVALUE; @@ -171,12 +175,13 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL CCTK_REAL PDstandardNth3g33 = INITVALUE; CCTK_REAL PDstandardNth11g33 = INITVALUE; CCTK_REAL PDstandardNth22g33 = INITVALUE; + CCTK_REAL PDstandardNth33g33 = INITVALUE; CCTK_REAL PDstandardNth12g33 = INITVALUE; CCTK_REAL PDstandardNth13g33 = INITVALUE; - CCTK_REAL PDstandardNth21g33 = INITVALUE; CCTK_REAL PDstandardNth23g33 = INITVALUE; CCTK_REAL PDstandardNth31g33 = INITVALUE; CCTK_REAL PDstandardNth32g33 = INITVALUE; + CCTK_REAL PDstandardNth1K11 = INITVALUE; CCTK_REAL PDstandardNth2K11 = INITVALUE; CCTK_REAL PDstandardNth3K11 = INITVALUE; CCTK_REAL PDstandardNth1K12 = INITVALUE; @@ -186,12 +191,14 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL CCTK_REAL PDstandardNth2K13 = INITVALUE; CCTK_REAL PDstandardNth3K13 = INITVALUE; CCTK_REAL PDstandardNth1K22 = INITVALUE; + CCTK_REAL PDstandardNth2K22 = INITVALUE; CCTK_REAL PDstandardNth3K22 = INITVALUE; CCTK_REAL PDstandardNth1K23 = INITVALUE; CCTK_REAL PDstandardNth2K23 = INITVALUE; CCTK_REAL PDstandardNth3K23 = INITVALUE; CCTK_REAL PDstandardNth1K33 = INITVALUE; CCTK_REAL PDstandardNth2K33 = INITVALUE; + CCTK_REAL PDstandardNth3K33 = INITVALUE; /* Assign local copies of grid functions */ g11L = g11[index]; @@ -215,12 +222,17 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL PDstandardNth1g11 = PDstandardNth1(g11, i, j, k); PDstandardNth2g11 = PDstandardNth2(g11, i, j, k); PDstandardNth3g11 = PDstandardNth3(g11, i, j, k); + PDstandardNth11g11 = PDstandardNth11(g11, i, j, k); PDstandardNth22g11 = PDstandardNth22(g11, i, j, k); PDstandardNth33g11 = PDstandardNth33(g11, i, j, k); + PDstandardNth12g11 = PDstandardNth12(g11, i, j, k); + PDstandardNth13g11 = PDstandardNth13(g11, i, j, k); PDstandardNth23g11 = PDstandardNth23(g11, i, j, k); PDstandardNth1g12 = PDstandardNth1(g12, i, j, k); PDstandardNth2g12 = PDstandardNth2(g12, i, j, k); PDstandardNth3g12 = PDstandardNth3(g12, i, j, k); + PDstandardNth11g12 = PDstandardNth11(g12, i, j, k); + PDstandardNth22g12 = PDstandardNth22(g12, i, j, k); PDstandardNth33g12 = PDstandardNth33(g12, i, j, k); PDstandardNth12g12 = PDstandardNth12(g12, i, j, k); PDstandardNth13g12 = PDstandardNth13(g12, i, j, k); @@ -228,7 +240,9 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL PDstandardNth1g13 = PDstandardNth1(g13, i, j, k); PDstandardNth2g13 = PDstandardNth2(g13, i, j, k); PDstandardNth3g13 = PDstandardNth3(g13, i, j, k); + PDstandardNth11g13 = PDstandardNth11(g13, i, j, k); PDstandardNth22g13 = PDstandardNth22(g13, i, j, k); + PDstandardNth33g13 = PDstandardNth33(g13, i, j, k); PDstandardNth12g13 = PDstandardNth12(g13, i, j, k); PDstandardNth13g13 = PDstandardNth13(g13, i, j, k); PDstandardNth23g13 = PDstandardNth23(g13, i, j, k); @@ -236,12 +250,17 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL PDstandardNth2g22 = PDstandardNth2(g22, i, j, k); PDstandardNth3g22 = PDstandardNth3(g22, i, j, k); PDstandardNth11g22 = PDstandardNth11(g22, i, j, k); + PDstandardNth22g22 = PDstandardNth22(g22, i, j, k); PDstandardNth33g22 = PDstandardNth33(g22, i, j, k); + PDstandardNth12g22 = PDstandardNth12(g22, i, j, k); PDstandardNth13g22 = PDstandardNth13(g22, i, j, k); + PDstandardNth23g22 = PDstandardNth23(g22, i, j, k); PDstandardNth1g23 = PDstandardNth1(g23, i, j, k); PDstandardNth2g23 = PDstandardNth2(g23, i, j, k); PDstandardNth3g23 = PDstandardNth3(g23, i, j, k); PDstandardNth11g23 = PDstandardNth11(g23, i, j, k); + PDstandardNth22g23 = PDstandardNth22(g23, i, j, k); + PDstandardNth33g23 = PDstandardNth33(g23, i, j, k); PDstandardNth12g23 = PDstandardNth12(g23, i, j, k); PDstandardNth13g23 = PDstandardNth13(g23, i, j, k); PDstandardNth23g23 = PDstandardNth23(g23, i, j, k); @@ -250,7 +269,11 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL PDstandardNth3g33 = PDstandardNth3(g33, i, j, k); PDstandardNth11g33 = PDstandardNth11(g33, i, j, k); PDstandardNth22g33 = PDstandardNth22(g33, i, j, k); + PDstandardNth33g33 = PDstandardNth33(g33, i, j, k); PDstandardNth12g33 = PDstandardNth12(g33, i, j, k); + PDstandardNth13g33 = PDstandardNth13(g33, i, j, k); + PDstandardNth23g33 = PDstandardNth23(g33, i, j, k); + PDstandardNth1K11 = PDstandardNth1(K11, i, j, k); PDstandardNth2K11 = PDstandardNth2(K11, i, j, k); PDstandardNth3K11 = PDstandardNth3(K11, i, j, k); PDstandardNth1K12 = PDstandardNth1(K12, i, j, k); @@ -260,170 +283,263 @@ void ML_ADM_constraints_Body(cGH *cctkGH, CCTK_INT dir, CCTK_INT face, CCTK_REAL PDstandardNth2K13 = PDstandardNth2(K13, i, j, k); PDstandardNth3K13 = PDstandardNth3(K13, i, j, k); PDstandardNth1K22 = PDstandardNth1(K22, i, j, k); + PDstandardNth2K22 = PDstandardNth2(K22, i, j, k); PDstandardNth3K22 = PDstandardNth3(K22, i, j, k); PDstandardNth1K23 = PDstandardNth1(K23, i, j, k); PDstandardNth2K23 = PDstandardNth2(K23, i, j, k); PDstandardNth3K23 = PDstandardNth3(K23, i, j, k); PDstandardNth1K33 = PDstandardNth1(K33, i, j, k); PDstandardNth2K33 = PDstandardNth2(K33, i, j, k); + PDstandardNth3K33 = PDstandardNth3(K33, i, j, k); /* Precompute derivatives (old style) */ /* Calculate temporaries and grid functions */ - detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - g11L*SQR(g23L); + detg = 2*g12L*g13L*g23L + g33L*(g11L*g22L - SQR(g12L)) - g22L*SQR(g13L) - + g11L*SQR(g23L); gu11 = INV(detg)*(g22L*g33L - SQR(g23L)); - gu21 = (g13L*g23L - g12L*g33L)*INV(detg); + gu12 = (g13L*g23L - g12L*g33L)*INV(detg); - gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + gu13 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + + gu21 = (g13L*g23L - g12L*g33L)*INV(detg); gu22 = INV(detg)*(g11L*g33L - SQR(g13L)); + gu23 = (g12L*g13L - g11L*g23L)*INV(detg); + + gu31 = (-(g13L*g22L) + g12L*g23L)*INV(detg); + gu32 = (g12L*g13L - g11L*g23L)*INV(detg); gu33 = INV(detg)*(g11L*g22L - SQR(g12L)); - G111 = khalf*(gu11*PDstandardNth1g11 + 2*(gu21*PDstandardNth1g12 + gu31*PDstandardNth1g13) - gu21*PDstandardNth2g11 - - gu31*PDstandardNth3g11); + G111 = khalf*(gu11*PDstandardNth1g11 + + 2*(gu12*PDstandardNth1g12 + gu13*PDstandardNth1g13) - + gu12*PDstandardNth2g11 - gu13*PDstandardNth3g11); - G211 = khalf*(gu21*PDstandardNth1g11 + 2*(gu22*PDstandardNth1g12 + gu32*PDstandardNth1g13) - gu22*PDstandardNth2g11 - - gu32*PDstandardNth3g11); + G211 = khalf*(gu21*PDstandardNth1g11 + + 2*(gu22*PDstandardNth1g12 + gu23*PDstandardNth1g13) - + gu22*PDstandardNth2g11 - gu23*PDstandardNth3g11); - G311 = khalf*(gu31*PDstandardNth1g11 + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - gu32*PDstandardNth2g11 - - gu33*PDstandardNth3g11); + G311 = khalf*(gu31*PDstandardNth1g11 + + 2*(gu32*PDstandardNth1g12 + gu33*PDstandardNth1g13) - + gu32*PDstandardNth2g11 - gu33*PDstandardNth3g11); - G112 = khalf*(gu21*PDstandardNth1g22 + gu11*PDstandardNth2g11 + - gu31*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); + G112 = khalf*(gu12*PDstandardNth1g22 + gu11*PDstandardNth2g11 + + gu13*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); G212 = khalf*(gu22*PDstandardNth1g22 + gu21*PDstandardNth2g11 + - gu32*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); + gu23*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); G312 = khalf*(gu32*PDstandardNth1g22 + gu31*PDstandardNth2g11 + gu33*(PDstandardNth1g23 + PDstandardNth2g13 - PDstandardNth3g12)); - G113 = khalf*(gu31*PDstandardNth1g33 + gu11*PDstandardNth3g11 + - gu21*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); + G113 = khalf*(gu13*PDstandardNth1g33 + gu11*PDstandardNth3g11 + + gu12*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G213 = khalf*(gu32*PDstandardNth1g33 + gu21*PDstandardNth3g11 + + G213 = khalf*(gu23*PDstandardNth1g33 + gu21*PDstandardNth3g11 + gu22*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); G313 = khalf*(gu33*PDstandardNth1g33 + gu31*PDstandardNth3g11 + gu32*(PDstandardNth1g23 - PDstandardNth2g13 + PDstandardNth3g12)); - G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu21*PDstandardNth2g22 + - gu31*(2*PDstandardNth2g23 - PDstandardNth3g22)); + G122 = khalf*(gu11*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + + gu12*PDstandardNth2g22 + gu13*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu22*PDstandardNth2g22 + - gu32*(2*PDstandardNth2g23 - PDstandardNth3g22)); + G222 = khalf*(gu21*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + + gu22*PDstandardNth2g22 + gu23*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + gu32*PDstandardNth2g22 + - gu33*(2*PDstandardNth2g23 - PDstandardNth3g22)); + G322 = khalf*(gu31*(-PDstandardNth1g22 + 2*PDstandardNth2g12) + + gu32*PDstandardNth2g22 + gu33*(2*PDstandardNth2g23 - PDstandardNth3g22)); - G123 = khalf*(gu31*PDstandardNth2g33 + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + - gu21*PDstandardNth3g22); + G123 = khalf*(gu13*PDstandardNth2g33 + + gu11*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + gu12*PDstandardNth3g22); - G223 = khalf*(gu32*PDstandardNth2g33 + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + G223 = khalf*(gu23*PDstandardNth2g33 + + gu21*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu22*PDstandardNth3g22); - G323 = khalf*(gu33*PDstandardNth2g33 + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + + G323 = khalf*(gu33*PDstandardNth2g33 + + gu31*(-PDstandardNth1g23 + PDstandardNth2g13 + PDstandardNth3g12) + gu32*PDstandardNth3g22); - G133 = khalf*(-(gu11*PDstandardNth1g33) - gu21*PDstandardNth2g33 + 2*gu11*PDstandardNth3g13 + - 2*gu21*PDstandardNth3g23 + gu31*PDstandardNth3g33); - - G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + 2*gu21*PDstandardNth3g13 + - 2*gu22*PDstandardNth3g23 + gu32*PDstandardNth3g33); - - G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + 2*gu31*PDstandardNth3g13 + - 2*gu32*PDstandardNth3g23 + gu33*PDstandardNth3g33); - - R11 = khalf*(-(gu22*PDstandardNth11g22) - 2*(G111*G122 + G111*G133 + G211*G222 + G211*G233 + G311*G322 + G311*G333 + - gu32*PDstandardNth11g23) - gu33*PDstandardNth11g33 + 2*gu22*PDstandardNth12g12 + 2*gu32*PDstandardNth12g13 + - 2*gu32*PDstandardNth13g12 + 2*gu33*PDstandardNth13g13 - gu22*PDstandardNth22g11 - 2*gu32*PDstandardNth23g11 - - gu33*PDstandardNth33g11 + 2*SQR(G112) + 2*SQR(G113) + 2*SQR(G212) + 2*SQR(G213) + 2*SQR(G312) + 2*SQR(G313)); - - R12 = khalf*(2*(G113*G123 + G213*G223 + G313*G323) - 2*(G112*G133 + G212*G233 + G312*G333 + gu21*PDstandardNth12g12) - - gu32*PDstandardNth12g23 - gu33*PDstandardNth12g33 + - gu31*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12) + gu32*PDstandardNth13g22 + - gu33*PDstandardNth13g23 + gu21*(PDstandardNth11g22 + PDstandardNth22g11) + gu32*PDstandardNth22g13 + - gu31*PDstandardNth23g11 - gu32*PDstandardNth23g12 + gu33*PDstandardNth23g13 - gu33*PDstandardNth33g12); - - R13 = khalf*(2*(G112*G123 + G212*G223 + G312*G323) - 2*(G113*G122 + G213*G222 + G313*G322 + gu31*PDstandardNth13g13) + - gu21*(PDstandardNth11g23 - PDstandardNth12g13 - PDstandardNth13g12 + PDstandardNth23g11) + - gu22*(PDstandardNth12g23 - PDstandardNth13g22 - PDstandardNth22g13 + PDstandardNth23g12) + - gu31*(PDstandardNth11g33 + PDstandardNth33g11) + - gu32*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12)); - - R22 = khalf*(-2*(G122*(G111 + G133) + G222*(G211 + G233) + G322*(G311 + G333) + gu31*PDstandardNth13g22) + - gu11*(-PDstandardNth11g22 + 2*PDstandardNth12g12 - PDstandardNth22g11) + - gu31*(-2*PDstandardNth22g13 + 2*(PDstandardNth12g23 + PDstandardNth23g12)) + - gu33*(-PDstandardNth22g33 + 2*PDstandardNth23g23 - PDstandardNth33g22) + - 2*(SQR(G112) + SQR(G123) + SQR(G212) + SQR(G223) + SQR(G312) + SQR(G323))); - - R23 = khalf*(2*(G112*G113 + G212*G213 + G312*G313) + - gu11*(-PDstandardNth11g23 + PDstandardNth12g13 + PDstandardNth13g12 - PDstandardNth23g11) + - gu21*(-PDstandardNth12g23 + PDstandardNth13g22 + PDstandardNth22g13 - PDstandardNth23g12) - - 2*(G111*G123 + G211*G223 + G311*G323 + gu32*PDstandardNth23g23) + - gu31*(PDstandardNth12g33 - PDstandardNth13g23 - PDstandardNth23g13 + PDstandardNth33g12) + - gu32*(PDstandardNth22g33 + PDstandardNth33g22)); - - R33 = khalf*(gu11*(-PDstandardNth11g33 + 2*PDstandardNth13g13 - PDstandardNth33g11) - - 2*((G111 + G122)*G133 + (G211 + G222)*G233 + (G311 + G322)*G333 + - gu21*(PDstandardNth12g33 + PDstandardNth33g12)) + - gu22*(-PDstandardNth22g33 + 2*PDstandardNth23g23 - PDstandardNth33g22) + - 2*(gu21*(PDstandardNth13g23 + PDstandardNth23g13) + SQR(G113) + SQR(G123) + SQR(G213) + SQR(G223) + SQR(G313) + - SQR(G323))); - - trR = gu11*R11 + gu22*R22 + 2*(gu21*R12 + gu31*R13 + gu32*R23) + gu33*R33; - - Km11 = gu11*K11L + gu21*K12L + gu31*K13L; - - Km21 = gu21*K11L + gu22*K12L + gu32*K13L; + G133 = khalf*(-(gu11*PDstandardNth1g33) - gu12*PDstandardNth2g33 + + 2*gu11*PDstandardNth3g13 + 2*gu12*PDstandardNth3g23 + + gu13*PDstandardNth3g33); + + G233 = khalf*(-(gu21*PDstandardNth1g33) - gu22*PDstandardNth2g33 + + 2*gu21*PDstandardNth3g13 + 2*gu22*PDstandardNth3g23 + + gu23*PDstandardNth3g33); + + G333 = khalf*(-(gu31*PDstandardNth1g33) - gu32*PDstandardNth2g33 + + 2*gu31*PDstandardNth3g13 + 2*gu32*PDstandardNth3g23 + + gu33*PDstandardNth3g33); + + R11 = -(G111*(G111 + G122 + G133)) - G211*(G211 + G222 + G233) - + G311*(G311 + G322 + G333) + khalf* + (-6*gu11*PDstandardNth11g11 + + gu12*(-8*PDstandardNth11g12 + 2*PDstandardNth12g12) + + gu22*(-9*PDstandardNth11g22 + 3*PDstandardNth12g22) + + gu13*(-8*PDstandardNth11g13 + 2*PDstandardNth13g13) + + gu33*(-9*PDstandardNth11g33 + 3*PDstandardNth13g33) + + gu21*(-6*PDstandardNth11g12 - PDstandardNth12g12 + PDstandardNth22g12) + + gu23*(-9*PDstandardNth11g23 + 3*PDstandardNth13g23 + + PDstandardNth22g23 - PDstandardNth23g23) + + gu31*(-6*PDstandardNth11g13 - PDstandardNth13g13 + PDstandardNth33g13) + + gu32*(-9*PDstandardNth11g23 + 3*PDstandardNth12g23 - + PDstandardNth23g23 + PDstandardNth33g23)) + SQR(G111) + SQR(G112) + + SQR(G113) + SQR(G211) + SQR(G212) + SQR(G213) + SQR(G311) + SQR(G312) + + SQR(G313); + + R12 = khalf*(-2*(G112*G133 + G212*G233 + G312*G333) - + 6*(gu11*PDstandardNth12g11 + gu31*PDstandardNth12g13) - + 9*(gu22*PDstandardNth12g22 + gu23*PDstandardNth12g23 + + gu32*PDstandardNth12g23 + gu33*PDstandardNth12g33) - + gu31*PDstandardNth13g13 + gu21* + (-7*PDstandardNth12g12 + PDstandardNth22g12) + + gu12*(PDstandardNth11g12 - 10*PDstandardNth12g12 + + 3*PDstandardNth22g12) + 3*gu22*PDstandardNth22g22 + + gu23*PDstandardNth22g23 + 3*gu32*PDstandardNth22g23 + + gu13*(PDstandardNth11g13 - 9*PDstandardNth12g13 - PDstandardNth13g13 + + 3*PDstandardNth23g13) - gu32*PDstandardNth23g23 + + 2*(G113*G123 + G213*G223 + G313*G323 + gu23*PDstandardNth23g23) + + 3*gu33*PDstandardNth23g33 + gu31*PDstandardNth33g13 + + gu32*PDstandardNth33g23); + + R13 = khalf*(-2*(G113*G122 + G213*G222 + G313*G322) - + 6*gu11*PDstandardNth13g11 - + 9*(gu22*PDstandardNth13g22 + (gu23 + gu32)*PDstandardNth13g23 + + gu33*PDstandardNth13g33) + + gu21*(-PDstandardNth12g12 - 6*PDstandardNth13g12 + PDstandardNth22g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 9*PDstandardNth13g12 + + 3*PDstandardNth23g12) + + gu23*(PDstandardNth22g23 - PDstandardNth23g23) + + 2*(G112*G123 + G212*G223 + G312*G323 + gu32*PDstandardNth23g23) + + gu31*(-7*PDstandardNth13g13 + PDstandardNth33g13) + + gu13*(PDstandardNth11g13 - 10*PDstandardNth13g13 + + 3*PDstandardNth33g13) + gu32*PDstandardNth33g23 + + 3*(gu22*PDstandardNth23g22 + gu23*PDstandardNth33g23 + + gu33*PDstandardNth33g33)); + + R22 = -(G122*(G111 + G122 + G133)) - G222*(G211 + G222 + G233) - + G322*(G311 + G322 + G333) + khalf* + (3*gu11*(PDstandardNth12g11 - 3*PDstandardNth22g11) + + gu21*(2*PDstandardNth12g12 - 8*PDstandardNth22g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth22g12) - + 6*gu22*PDstandardNth22g22 + + gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + + 3*PDstandardNth23g13) - + gu23*(8*PDstandardNth22g23 - 2*PDstandardNth23g23) + + gu33*(-9*PDstandardNth22g33 + 3*PDstandardNth23g33) + + gu31*(3*PDstandardNth12g13 - PDstandardNth13g13 - 9*PDstandardNth22g13 + + PDstandardNth33g13) + gu32* + (-6*PDstandardNth22g23 - PDstandardNth23g23 + PDstandardNth33g23)) + + SQR(G112) + SQR(G122) + SQR(G123) + SQR(G212) + SQR(G222) + SQR(G223) + + SQR(G312) + SQR(G322) + SQR(G323); + + R23 = khalf*(-2*(G111*G123 + G211*G223 + G311*G323) + + 2*(G112*G113 + G212*G213 + G312*G313 + gu31*PDstandardNth13g13) + + gu11*(3*PDstandardNth13g11 - 9*PDstandardNth23g11) + + gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - + 9*PDstandardNth23g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 - 6*PDstandardNth23g12) - + 6*gu22*PDstandardNth23g22 - + 9*(gu31*PDstandardNth23g13 + gu33*PDstandardNth23g33) + + gu31*PDstandardNth33g13 + gu13* + (PDstandardNth11g13 - PDstandardNth13g13 - 9*PDstandardNth23g13 + + 3*PDstandardNth33g13) + + gu32*(-7*PDstandardNth23g23 + PDstandardNth33g23) + + gu23*(PDstandardNth22g23 - 10*PDstandardNth23g23 + + 3*PDstandardNth33g23) + 3*gu33*PDstandardNth33g33); + + R33 = -(G133*(G111 + G122 + G133)) - G233*(G211 + G222 + G233) - + G333*(G311 + G322 + G333) + khalf* + (3*gu11*(PDstandardNth13g11 - 3*PDstandardNth33g11) + + gu21*(-PDstandardNth12g12 + 3*PDstandardNth13g12 + PDstandardNth22g12 - + 9*PDstandardNth33g12) + + gu12*(PDstandardNth11g12 - PDstandardNth12g12 + 3*PDstandardNth23g12 - + 9*PDstandardNth33g12) + + gu31*(2*PDstandardNth13g13 - 8*PDstandardNth33g13) + + gu13*(PDstandardNth11g13 - PDstandardNth13g13 - 6*PDstandardNth33g13) + + 3*gu22*(PDstandardNth23g22 - 3*PDstandardNth33g22) + + gu32*(2*PDstandardNth23g23 - 8*PDstandardNth33g23) + + gu23*(PDstandardNth22g23 - PDstandardNth23g23 - 6*PDstandardNth33g23) - + 6*gu33*PDstandardNth33g33) + SQR(G113) + SQR(G123) + SQR(G133) + + SQR(G213) + SQR(G223) + SQR(G233) + SQR(G313) + SQR(G323) + SQR(G333); + + trR = gu11*R11 + (gu12 + gu21)*R12 + (gu13 + gu31)*R13 + gu22*R22 + + (gu23 + gu32)*R23 + gu33*R33; + + Km11 = gu11*K11L + gu12*K12L + gu13*K13L; + + Km21 = gu21*K11L + gu22*K12L + gu23*K13L; Km31 = gu31*K11L + gu32*K12L + gu33*K13L; - Km12 = gu11*K12L + gu21*K22L + gu31*K23L; + Km12 = gu11*K12L + gu12*K22L + gu13*K23L; - Km22 = gu21*K12L + gu22*K22L + gu32*K23L; + Km22 = gu21*K12L + gu22*K22L + gu23*K23L; Km32 = gu31*K12L + gu32*K22L + gu33*K23L; - Km13 = gu11*K13L + gu21*K23L + gu31*K33L; + Km13 = gu11*K13L + gu12*K23L + gu13*K33L; - Km23 = gu21*K13L + gu22*K23L + gu32*K33L; + Km23 = gu21*K13L + gu22*K23L + gu23*K33L; Km33 = gu31*K13L + gu32*K23L + gu33*K33L; trK = Km11 + Km22 + Km33; - HL = -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - SQR(Km33) + SQR(trK); - - M1L = gu21*(-(G112*K11L) + G111*K12L - G212*K12L - G312*K13L + G211*K22L + G311*K23L - PDstandardNth1K12 + - PDstandardNth2K11) + gu22*(-(G122*K11L) + G112*K12L - G222*K12L - G322*K13L + G212*K22L + G312*K23L - - PDstandardNth1K22 + PDstandardNth2K12) + gu31* - (-(G113*K11L) - G213*K12L + G111*K13L - G313*K13L + G211*K23L + G311*K33L - PDstandardNth1K13 + PDstandardNth3K11)\ - + gu32*(G113*K12L + G112*K13L + G213*K22L + (G212 + G313)*K23L + G312*K33L - - 2*(G123*K11L + G223*K12L + G323*K13L + PDstandardNth1K23) + PDstandardNth2K13 + PDstandardNth3K12) + - gu33*(-(G133*K11L) - G233*K12L + G113*K13L - G333*K13L + G213*K23L + G313*K33L - PDstandardNth1K33 + - PDstandardNth3K13); - - M2L = gu11*(G112*K11L + (-G111 + G212)*K12L + G312*K13L - G211*K22L - G311*K23L + PDstandardNth1K12 - - PDstandardNth2K11) + gu21*(G122*K11L + (-G112 + G222)*K12L + G322*K13L - G212*K22L - G312*K23L + - PDstandardNth1K22 - PDstandardNth2K12) + gu31* - (G123*K11L + (-2*G113 + G223)*K12L + (G112 + G323)*K13L + G212*K23L + G312*K33L + PDstandardNth1K23 - - 2*(G213*K22L + G313*K23L + PDstandardNth2K13) + PDstandardNth3K12) + - gu32*(-(G123*K12L) + G122*K13L - G223*K22L + G222*K23L - G323*K23L + G322*K33L - PDstandardNth2K23 + - PDstandardNth3K22) + gu33*(-(G133*K12L) + G123*K13L - G233*K22L + G223*K23L - G333*K23L + G323*K33L - - PDstandardNth2K33 + PDstandardNth3K23); - - M3L = gu11*(G113*K11L + G213*K12L + (-G111 + G313)*K13L - G211*K23L - G311*K33L + PDstandardNth1K13 - - PDstandardNth3K11) + gu21*(G123*K11L + (G113 + G223)*K12L + (-2*G112 + G323)*K13L + G213*K22L + - (-2*G212 + G313)*K23L + PDstandardNth1K23 + PDstandardNth2K13 - 2*(G312*K33L + PDstandardNth3K12)) + - gu31*(G133*K11L + G233*K12L + (-G113 + G333)*K13L - G213*K23L - G313*K33L + PDstandardNth1K33 - - PDstandardNth3K13) + gu22*(G123*K12L - G122*K13L + G223*K22L - G222*K23L + G323*K23L - G322*K33L + - PDstandardNth2K23 - PDstandardNth3K22) + gu32* - (G133*K12L - G123*K13L + G233*K22L - G223*K23L + G333*K23L - G323*K33L + PDstandardNth2K33 - PDstandardNth3K23); + HL = -2*(Km12*Km21 + Km13*Km31 + Km23*Km32) + trR - SQR(Km11) - SQR(Km22) - + SQR(Km33) + SQR(trK); + + M1L = -2*(gu11*PDstandardNth1K11 + gu12*PDstandardNth1K12 + + gu13*PDstandardNth1K13) + + gu21*(-(G112*K11L) + G111*K12L - G212*K12L - G312*K13L + G211*K22L + + G311*K23L - 3*PDstandardNth1K12 + PDstandardNth2K12) + + gu22*(-(G122*K11L) + G112*K12L - G222*K12L - G322*K13L + G212*K22L + + G312*K23L - 3*PDstandardNth1K22 + PDstandardNth2K22) + + gu23*(-(G123*K11L) + G113*K12L - G223*K12L - G323*K13L + G213*K22L + + G313*K23L - 3*PDstandardNth1K23 + PDstandardNth2K23) + + gu31*(-(G113*K11L) - G213*K12L + G111*K13L - G313*K13L + G211*K23L + + G311*K33L - 3*PDstandardNth1K13 + PDstandardNth3K13) + + gu32*(-(G123*K11L) - G223*K12L + G112*K13L - G323*K13L + G212*K23L + + G312*K33L - 3*PDstandardNth1K23 + PDstandardNth3K23) + + gu33*(-(G133*K11L) - G233*K12L + G113*K13L - G333*K13L + G213*K23L + + G313*K33L - 3*PDstandardNth1K33 + PDstandardNth3K33); + + M2L = gu11*(G112*K11L + (-G111 + G212)*K12L + G312*K13L - G211*K22L - + G311*K23L + PDstandardNth1K11 - 3*PDstandardNth2K11) + + gu12*(G122*K11L + (-G112 + G222)*K12L + G322*K13L - G212*K22L - G312*K23L + + PDstandardNth1K12 - 3*PDstandardNth2K12) + + gu13*(G123*K11L + (-G113 + G223)*K12L + G323*K13L - G213*K22L - G313*K23L + + PDstandardNth1K13 - 3*PDstandardNth2K13) - + 2*(gu21*PDstandardNth2K12 + gu22*PDstandardNth2K22 + + gu23*PDstandardNth2K23) + + gu31*(-(G113*K12L) + G112*K13L - G213*K22L + G212*K23L - G313*K23L + + G312*K33L - 3*PDstandardNth2K13 + PDstandardNth3K13) + + gu32*(-(G123*K12L) + G122*K13L - G223*K22L + G222*K23L - G323*K23L + + G322*K33L - 3*PDstandardNth2K23 + PDstandardNth3K23) + + gu33*(-(G133*K12L) + G123*K13L - G233*K22L + G223*K23L - G333*K23L + + G323*K33L - 3*PDstandardNth2K33 + PDstandardNth3K33); + + M3L = gu11*(G113*K11L + G213*K12L + (-G111 + G313)*K13L - G211*K23L - + G311*K33L + PDstandardNth1K11 - 3*PDstandardNth3K11) + + gu12*(G123*K11L + G223*K12L + (-G112 + G323)*K13L - G212*K23L - G312*K33L + + PDstandardNth1K12 - 3*PDstandardNth3K12) + + gu21*(G113*K12L - G112*K13L + G213*K22L - G212*K23L + G313*K23L - + G312*K33L + PDstandardNth2K12 - 3*PDstandardNth3K12) + + gu13*(G133*K11L + G233*K12L + (-G113 + G333)*K13L - G213*K23L - G313*K33L + + PDstandardNth1K13 - 3*PDstandardNth3K13) + + gu22*(G123*K12L - G122*K13L + G223*K22L - G222*K23L + G323*K23L - + G322*K33L + PDstandardNth2K22 - 3*PDstandardNth3K22) + + gu23*(G133*K12L - G123*K13L + G233*K22L - G223*K23L + G333*K23L - + G323*K33L + PDstandardNth2K23 - 3*PDstandardNth3K23) - + 2*(gu31*PDstandardNth3K13 + gu32*PDstandardNth3K23 + gu33*PDstandardNth3K33); /* Copy local copies back to grid functions */ -- cgit v1.2.3