From c7f22fc40fcd6d905975f4588c31ce6d5cf0aba5 Mon Sep 17 00:00:00 2001 From: Liwei Ji Date: Wed, 24 Jul 2024 16:53:00 -0500 Subject: [PATCH] Z4coo: remove thorn --- Z4coo/README.md | 12 - Z4coo/configuration.ccl | 3 - Z4coo/interface.ccl | 52 - Z4coo/param.ccl | 103 -- Z4coo/schedule.ccl | 222 ---- Z4coo/src/adm.cxx | 90 -- Z4coo/src/constraint.cxx | 192 --- Z4coo/src/derivs.hxx | 492 -------- Z4coo/src/enforce.cxx | 128 -- Z4coo/src/initial1.cxx | 131 -- Z4coo/src/initial2.cxx | 108 -- Z4coo/src/make.code.defn | 14 - Z4coo/src/newradx.cxx | 36 - Z4coo/src/rhs.cxx | 253 ---- Z4coo/wolfram/Z4coo_set_ADM.hxx | 126 -- Z4coo/wolfram/Z4coo_set_constraint.hxx | 1034 ---------------- Z4coo/wolfram/Z4coo_set_rhs.hxx | 1523 ------------------------ 17 files changed, 4519 deletions(-) delete mode 100644 Z4coo/README.md delete mode 100644 Z4coo/configuration.ccl delete mode 100644 Z4coo/interface.ccl delete mode 100644 Z4coo/param.ccl delete mode 100644 Z4coo/schedule.ccl delete mode 100644 Z4coo/src/adm.cxx delete mode 100644 Z4coo/src/constraint.cxx delete mode 100644 Z4coo/src/derivs.hxx delete mode 100644 Z4coo/src/enforce.cxx delete mode 100644 Z4coo/src/initial1.cxx delete mode 100644 Z4coo/src/initial2.cxx delete mode 100644 Z4coo/src/make.code.defn delete mode 100644 Z4coo/src/newradx.cxx delete mode 100644 Z4coo/src/rhs.cxx delete mode 100644 Z4coo/wolfram/Z4coo_set_ADM.hxx delete mode 100644 Z4coo/wolfram/Z4coo_set_constraint.hxx delete mode 100644 Z4coo/wolfram/Z4coo_set_rhs.hxx diff --git a/Z4coo/README.md b/Z4coo/README.md deleted file mode 100644 index 26d5433c..00000000 --- a/Z4coo/README.md +++ /dev/null @@ -1,12 +0,0 @@ -# Z4coo - -|Authors |Liwei Ji | -|:----------------|:---------------------------| -|Maintainer |Liwei Ji -|Licence |GNU GPL version 2 - -## How to generate header files - -* Install [Generato](https://github.com/lwJi/Generato) - -* Run `wolframscript -f Z4coo_set_rhs.wl` diff --git a/Z4coo/configuration.ccl b/Z4coo/configuration.ccl deleted file mode 100644 index c68469d1..00000000 --- a/Z4coo/configuration.ccl +++ /dev/null @@ -1,3 +0,0 @@ -# Configuration definitions for thorn Z4coo - -REQUIRES Arith Loop NewRadX diff --git a/Z4coo/interface.ccl b/Z4coo/interface.ccl deleted file mode 100644 index 47ce8f41..00000000 --- a/Z4coo/interface.ccl +++ /dev/null @@ -1,52 +0,0 @@ -# Interface definition for thorn Z4coo - -IMPLEMENTS: Z4coo - -INHERITS: ADMBaseX TmunuBaseX - -USES INCLUDE HEADER: defs.hxx -USES INCLUDE HEADER: div.hxx -USES INCLUDE HEADER: dual.hxx -USES INCLUDE HEADER: loop_device.hxx -USES INCLUDE HEADER: mat.hxx -USES INCLUDE HEADER: simd.hxx -USES INCLUDE HEADER: sum.hxx -USES INCLUDE HEADER: vec.hxx -USES INCLUDE HEADER: vect.hxx - -USES INCLUDE HEADER: newradx.hxx - - -CCTK_INT FUNCTION GetCallFunctionCount() -REQUIRES FUNCTION GetCallFunctionCount - - - -# All variables have been shifted so that they tend to zero in flat space - -CCTK_REAL chi TYPE=gf TAGS='rhs="chi_rhs" dependents="ADMBaseX::metric"' "chi" -CCTK_REAL gamma_tilde TYPE=gf TAGS='parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} rhs="gamma_tilde_rhs" dependents="ADMBaseX::metric"' { gammatxx gammatxy gammatxz gammatyy gammatyz gammatzz } "gamma-tilde" -CCTK_REAL K_hat TYPE=gf TAGS='rhs="K_hat_rhs" dependents="ADMBaseX::curv"' { Kh } "K-hat" -CCTK_REAL A_tilde TYPE=gf TAGS='parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} rhs="A_tilde_rhs" dependents="ADMBaseX::curv"' { Atxx Atxy Atxz Atyy Atyz Atzz } "A-tilde" -CCTK_REAL Gam_tilde TYPE=gf TAGS='parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} rhs="Gam_tilde_rhs"' { Gamtx Gamty Gamtz } "Gamma-tilde" -CCTK_REAL Theta TYPE=gf TAGS='rhs="Theta_rhs" dependents="ADMBaseX::curv"' "Theta" -CCTK_REAL alphaG TYPE=gf TAGS='rhs="alphaG_rhs" dependents="ADMBaseX::lapse ADMBaseX::dtlapse"' "alpha" -CCTK_REAL betaG TYPE=gf TAGS='parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} rhs="betaG_rhs" dependents="ADMBaseX::shift ADMBaseX::dtshift"' { betaGx betaGy betaGz } "beta" - - - -CCTK_REAL ZtC TYPE=gf TAGS='parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} checkpoint="no"' { ZtCx ZtCy ZtCz } "Z-tilde" -CCTK_REAL HC TYPE=gf TAGS='checkpoint="no"' "H" -CCTK_REAL MtC TYPE=gf TAGS='parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} checkpoint="no"' { MtCx MtCy MtCz } "M-tilde" -CCTK_REAL allC TYPE=gf TAGS='checkpoint="no"' "constraint monitor" - - - -CCTK_REAL chi_rhs TYPE=gf TAGS='checkpoint="no"' "chi" -CCTK_REAL gamma_tilde_rhs TYPE=gf TAGS='parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} checkpoint="no"' { gammatxx_rhs gammatxy_rhs gammatxz_rhs gammatyy_rhs gammatyz_rhs gammatzz_rhs } "gamma-tilde" -CCTK_REAL K_hat_rhs TYPE=gf TAGS='checkpoint="no"' { Kh_rhs } "K-hat" -CCTK_REAL A_tilde_rhs TYPE=gf TAGS='parities={+1 +1 +1 -1 -1 +1 -1 +1 -1 +1 +1 +1 +1 -1 -1 +1 +1 +1} checkpoint="no"' { Atxx_rhs Atxy_rhs Atxz_rhs Atyy_rhs Atyz_rhs Atzz_rhs } "A-tilde" -CCTK_REAL Gam_tilde_rhs TYPE=gf TAGS='parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} checkpoint="no"' { Gamtx_rhs Gamty_rhs Gamtz_rhs } "Gamma-tilde" -CCTK_REAL Theta_rhs TYPE=gf TAGS='checkpoint="no"' "Theta" -CCTK_REAL alphaG_rhs TYPE=gf TAGS='checkpoint="no"' "alpha" -CCTK_REAL betaG_rhs TYPE=gf TAGS='parities={-1 +1 +1 +1 -1 +1 +1 +1 -1} checkpoint="no"' { betaGx_rhs betaGy_rhs betaGz_rhs } "beta" diff --git a/Z4coo/param.ccl b/Z4coo/param.ccl deleted file mode 100644 index 22139125..00000000 --- a/Z4coo/param.ccl +++ /dev/null @@ -1,103 +0,0 @@ -# Parameter definitions for thorn Z4coo - -BOOLEAN calc_ADM_vars "Calculate ADM variables" STEERABLE=recover -{ -} yes - -BOOLEAN calc_ADMRHS_vars "Calculate RHS of ADM variables" STEERABLE=recover -{ -} yes - -BOOLEAN calc_constraints "Calculate constraints" STEERABLE=recover -{ -} yes - - -BOOLEAN set_Theta_zero "set Theta to zero, which converts Z4c to BSSN" -{ -} no - -CCTK_REAL kappa1 "kappa1" STEERABLE=always -{ - *:* :: "" -} 0.02 - -CCTK_REAL kappa2 "kappa2" STEERABLE=always -{ - *:* :: "" -} 0.0 - -CCTK_REAL f_mu_L "mu_L = f_mu_L / alpha" STEERABLE=always -{ - *:* :: "" -} 2.0 - -CCTK_REAL f_mu_S "mu_S = f_mu_S / alpha^2" STEERABLE=always -{ - *:* :: "" -} 1.0 - -CCTK_REAL eta "eta" STEERABLE=always -{ - *:* :: "" -} 2.0 - -CCTK_REAL chi_floor "Floor for chi" STEERABLE=always -{ - (0:* :: "" -} 1.0e-10 - -CCTK_REAL alphaG_floor "Floor for alphaG" STEERABLE=always -{ - (0:* :: "" -} 1.0e-10 - -CCTK_REAL epsdiss "Dissipation coefficient " STEERABLE=always -{ - 0.0:* :: "" -} 0.32 - -KEYWORD boundary_conditions "boundary conditions"{ - "NewRadX" :: "radiative boundary conditions using NewRadX" - "CarpetX" :: "use CarpetX default boundary conditions" -} "CarpetX" - -CCTK_INT n_chi "n power of outgoing boundary r^n fall off rate for chi" -{ - 0:2 :: "1 is reasonable" -} 1 - -CCTK_INT n_gammat "n power of outgoing boundary r^n fall off rate for gammat_ij" -{ - 0:2 :: "1 is reasonable" -} 1 - -CCTK_INT n_Kh "n power of outgoing boundary r^n fall off rate for Kh" -{ - 0:2 :: "Maybe 1?" -} 1 - -CCTK_INT n_At "n power of outgoing boundary r^n fall off rate for At_ij" -{ - 0:2 :: "Maybe 1?" -} 1 - -CCTK_INT n_Gamt "n power of outgoing boundary r^n fall off rate for Gamt^i" -{ - 0:2 :: "Maybe 1?" -} 1 - -CCTK_INT n_Theta "n power of outgoing boundary r^n fall off rate for Theta" -{ - 0:2 :: "Maybe 1?" -} 1 - -CCTK_INT n_alphaG "n power of outgoing boundary r^n fall off rate for alpha" -{ - 0:2 :: "1 is my guess" -} 1 - -CCTK_INT n_betaG "n power of outgoing boundary r^n fall off rate for beta" -{ - 0:2 :: "1 is my guess" -} 1 diff --git a/Z4coo/schedule.ccl b/Z4coo/schedule.ccl deleted file mode 100644 index a7d03811..00000000 --- a/Z4coo/schedule.ccl +++ /dev/null @@ -1,222 +0,0 @@ -# Schedule definitions for thorn Z4coo - -STORAGE: chi -STORAGE: gamma_tilde -STORAGE: K_hat -STORAGE: A_tilde -STORAGE: Gam_tilde -STORAGE: Theta -STORAGE: alphaG -STORAGE: betaG - -STORAGE: ZtC -STORAGE: HC -STORAGE: MtC -STORAGE: allC - -STORAGE: chi_rhs -STORAGE: gamma_tilde_rhs -STORAGE: K_hat_rhs -STORAGE: A_tilde_rhs -STORAGE: Gam_tilde_rhs -STORAGE: Theta_rhs -STORAGE: alphaG_rhs -STORAGE: betaG_rhs - - -################################################################################ -# We have 4 schedule groups: -# 1. initial: set up core Z4c variables from ADM variables -# 2. poststep: post-process core Z4c variables and calculate other variables -# 3. analysis: calculate constraints etc. -# 4. rhs: calculate RHS of Z4c variables - -SCHEDULE GROUP Z4coo_InitialGroup AT initial AFTER ADMBaseX_PostInitial -{ -} "Convert ADM to Z4c variables" - -SCHEDULE GROUP Z4coo_PostStepGroup AT initial AFTER Z4coo_InitialGroup BEFORE ADMBaseX_SetADMVars -{ -} "Post-process Z4c variables" - -SCHEDULE GROUP Z4coo_PostStepGroup2 AT initial AFTER (TmunuBaseX_SetTmunuVars, Z4coo_PostStepGroup) BEFORE ADMBaseX_SetADMRHS -{ -} "Post-process Z4c variables, part 2" - -SCHEDULE GROUP Z4coo_PostStepGroup AT postregrid BEFORE ADMBaseX_SetADMVars -{ -} "Post-process Z4c variables" - -SCHEDULE GROUP Z4coo_PostStepGroup2 AT postregrid AFTER (TmunuBaseX_SetTmunuVars, Z4coo_PostStepGroup) BEFORE ADMBaseX_SetADMRHS -{ -} "Post-process Z4c variables, part 2" - -SCHEDULE GROUP Z4coo_PostStepGroup IN ODESolvers_PostStep BEFORE ADMBaseX_SetADMVars -{ -} "Post-process Z4c variables" - -SCHEDULE GROUP Z4coo_PostStepGroup2 IN ODESolvers_PostStep AFTER (TmunuBaseX_SetTmunuVars, Z4coo_PostStepGroup) BEFORE ADMBaseX_SetADMRHS -{ -} "Post-process Z4c variables, part 2" - -SCHEDULE GROUP Z4coo_AnalysisGroup AT analysis -{ -} "Analyse Z4c variables" - -SCHEDULE GROUP Z4coo_RHSGroup IN ODESolvers_RHS -{ -} "Calculate Z4c RHS" - -################################################################################ - - -SCHEDULE Z4coo_Initial1 IN Z4coo_InitialGroup -{ - LANG: C - READS: ADMBaseX::metric(interior) - READS: ADMBaseX::curv(interior) - READS: ADMBaseX::lapse(interior) - READS: ADMBaseX::shift(interior) - WRITES: chi(interior) - WRITES: gamma_tilde(interior) - WRITES: K_hat(interior) - WRITES: A_tilde(interior) - WRITES: Theta(interior) - WRITES: alphaG(interior) - WRITES: betaG(interior) - SYNC: gamma_tilde -} "Convert ADM to Z4c variables, part 1" - -SCHEDULE Z4coo_Initial2 IN Z4coo_InitialGroup AFTER Z4coo_Initial1 -{ - LANG: C - READS: gamma_tilde(everywhere) - WRITES: Gam_tilde(interior) -} "Convert ADM to Z4c variables, part 2" - -SCHEDULE Z4coo_Enforce IN Z4coo_PostStepGroup -{ - LANG: C - READS: chi(interior) - READS: gamma_tilde(interior) - READS: A_tilde(interior) - READS: alphaG(interior) - WRITES: chi(interior) - WRITES: gamma_tilde(interior) - WRITES: A_tilde(interior) - WRITES: alphaG(interior) - SYNC: chi - SYNC: gamma_tilde - SYNC: K_hat - SYNC: A_tilde - SYNC: Gam_tilde - SYNC: Theta - SYNC: alphaG - SYNC: betaG -} "Enforce algebraic Z4c constraints" - -if (calc_ADM_vars) { - SCHEDULE Z4coo_ADM IN Z4coo_PostStepGroup AFTER Z4coo_Enforce - { - LANG: C - READS: chi(everywhere) - READS: gamma_tilde(everywhere) - READS: K_hat(everywhere) - READS: A_tilde(everywhere) - #READS: Gam_tilde(everywhere) - READS: Theta(everywhere) - READS: alphaG(everywhere) - READS: betaG(everywhere) - WRITES: ADMBaseX::metric(everywhere) - WRITES: ADMBaseX::curv(everywhere) - WRITES: ADMBaseX::lapse(everywhere) - WRITES: ADMBaseX::shift(everywhere) - #WRITES: ADMBaseX::dtlapse(everywhere) - #WRITES: ADMBaseX::dtshift(everywhere) - } "Convert Z4c to ADM variables" -} - -if (calc_constraints) { - SCHEDULE Z4coo_Constraints IN Z4coo_AnalysisGroup - { - LANG: C - READS: chi(everywhere) - READS: gamma_tilde(everywhere) - READS: K_hat(everywhere) - READS: A_tilde(everywhere) - READS: Gam_tilde(everywhere) - READS: Theta(everywhere) - READS: alphaG(everywhere) - READS: betaG(everywhere) - READS: TmunuBaseX::eTtt(interior) - READS: TmunuBaseX::eTti(interior) - READS: TmunuBaseX::eTij(interior) - WRITES: ZtC(interior) - WRITES: HC(interior) - WRITES: MtC(interior) - #WRITES: allC(interior) - } "Calculate Z4c Constraints" -} - -SCHEDULE Z4coo_RHS IN Z4coo_RHSGroup -{ - LANG: C - READS: chi(everywhere) - READS: gamma_tilde(everywhere) - READS: K_hat(everywhere) - READS: A_tilde(everywhere) - READS: Gam_tilde(everywhere) - READS: Theta(everywhere) - READS: alphaG(everywhere) - READS: betaG(everywhere) - READS: TmunuBaseX::eTtt(interior) - READS: TmunuBaseX::eTti(interior) - READS: TmunuBaseX::eTij(interior) - WRITES: chi_rhs(interior) - WRITES: gamma_tilde_rhs(interior) - WRITES: K_hat_rhs(interior) - WRITES: A_tilde_rhs(interior) - WRITES: Gam_tilde_rhs(interior) - WRITES: Theta_rhs(interior) - WRITES: alphaG_rhs(interior) - WRITES: betaG_rhs(interior) - SYNC: chi_rhs - SYNC: gamma_tilde_rhs - SYNC: K_hat_rhs - SYNC: A_tilde_rhs - SYNC: Gam_tilde_rhs - SYNC: Theta_rhs - SYNC: alphaG_rhs - SYNC: betaG_rhs -} "Calculate Z4c RHS" - -if (CCTK_Equals(boundary_conditions, "NewRadX")) { - SCHEDULE Z4coo_Apply_NewRadX_BC IN Z4coo_RHSGroup AFTER Z4c_RHS - { - LANG: C - READS: chi(everywhere) - READS: gamma_tilde(everywhere) - READS: K_hat(everywhere) - READS: A_tilde(everywhere) - READS: Gam_tilde(everywhere) - READS: Theta(everywhere) - READS: alphaG(everywhere) - READS: betaG(everywhere) - READS: chi_rhs(everywhere) - READS: gamma_tilde_rhs(everywhere) - READS: K_hat_rhs(everywhere) - READS: A_tilde_rhs(everywhere) - READS: Gam_tilde_rhs(everywhere) - READS: Theta_rhs(everywhere) - READS: alphaG_rhs(everywhere) - READS: betaG_rhs(everywhere) - WRITES: chi_rhs(interior) - WRITES: gamma_tilde_rhs(interior) - WRITES: K_hat_rhs(interior) - WRITES: A_tilde_rhs(interior) - WRITES: Gam_tilde_rhs(interior) - WRITES: Theta_rhs(interior) - WRITES: alphaG_rhs(interior) - WRITES: betaG_rhs(interior) - } "Apply radiative boundary conditions to Z4c RHS variables using NewRadX" -} diff --git a/Z4coo/src/adm.cxx b/Z4coo/src/adm.cxx deleted file mode 100644 index 2c732321..00000000 --- a/Z4coo/src/adm.cxx +++ /dev/null @@ -1,90 +0,0 @@ -#include -#include -#include -#include - -#include -#include -#include - -#ifdef __CUDACC__ -#include -#endif - -#include - -namespace Z4coo { -using namespace Arith; -using namespace Loop; -using namespace std; - -extern "C" void Z4coo_ADM(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS_Z4coo_ADM; - DECLARE_CCTK_PARAMETERS; - - const array indextype = {0, 0, 0}; - const GF3D2layout layout2(cctkGH, indextype); - - const GF3D2 gf_chi(layout2, chi); - - const smat, 3> gf_gamt{ - GF3D2(layout2, gammatxx), - GF3D2(layout2, gammatxy), - GF3D2(layout2, gammatxz), - GF3D2(layout2, gammatyy), - GF3D2(layout2, gammatyz), - GF3D2(layout2, gammatzz)}; - - const GF3D2 gf_exKh(layout2, Kh); - - const smat, 3> gf_exAt{ - GF3D2(layout2, Atxx), - GF3D2(layout2, Atxy), - GF3D2(layout2, Atxz), - GF3D2(layout2, Atyy), - GF3D2(layout2, Atyz), - GF3D2(layout2, Atzz)}; - - const GF3D2 gf_Theta(layout2, Theta); - - const GF3D2 gf_alpha(layout2, alphaG); - - const vec, 3> gf_beta{ - GF3D2(layout2, betaGx), - GF3D2(layout2, betaGy), - GF3D2(layout2, betaGz)}; - - const smat, 3> gf_ADMgam{ - GF3D2(layout2, gxx), GF3D2(layout2, gxy), - GF3D2(layout2, gxz), GF3D2(layout2, gyy), - GF3D2(layout2, gyz), GF3D2(layout2, gzz)}; - - const smat, 3> gf_ADMK{ - GF3D2(layout2, kxx), GF3D2(layout2, kxy), - GF3D2(layout2, kxz), GF3D2(layout2, kyy), - GF3D2(layout2, kyz), GF3D2(layout2, kzz)}; - - const GF3D2 gf_ADMalpha(layout2, alp); - - const vec, 3> gf_ADMbeta{GF3D2(layout2, betax), - GF3D2(layout2, betay), - GF3D2(layout2, betaz)}; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const Loop::GridDescBaseDevice grid(cctkGH); - -#ifdef __CUDACC__ - const nvtxRangeId_t range = nvtxRangeStartA("Z4coo_ADM::adm"); -#endif - -#include "../wolfram/Z4coo_set_ADM.hxx" - -#ifdef __CUDACC__ - nvtxRangeEnd(range); -#endif -} - -} // namespace Z4coo diff --git a/Z4coo/src/constraint.cxx b/Z4coo/src/constraint.cxx deleted file mode 100644 index af4fd4f2..00000000 --- a/Z4coo/src/constraint.cxx +++ /dev/null @@ -1,192 +0,0 @@ -#include - -#ifdef __CUDACC__ -// Disable CCTK_DEBUG since the debug information takes too much -// parameter space to launch the kernels -#ifdef CCTK_DEBUG -#undef CCTK_DEBUG -#endif -#endif - -#include "derivs.hxx" - -#include -#include - -#include -#include -#include - -#ifdef __CUDACC__ -#include -#endif - -#include - -namespace Z4coo { -using namespace Arith; -using namespace Loop; -using namespace std; - -template inline T Power(T x, int y) { - return (y == 2) ? Arith::pow2(x) : Arith::pown(x, y); -} - -extern "C" void Z4coo_Constraints(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS_Z4coo_Constraints; - DECLARE_CCTK_PARAMETERS; - - for (int d = 0; d < 3; ++d) - if (cctk_nghostzones[d] < deriv_order / 2 + 1) - CCTK_VERROR("Need at least %d ghost zones", deriv_order / 2 + 1); - - // - - const array indextype = {0, 0, 0}; - const array nghostzones = {cctk_nghostzones[0], cctk_nghostzones[1], - cctk_nghostzones[2]}; - vect imin, imax; - GridDescBase(cctkGH).box_int<0, 0, 0>(nghostzones, imin, imax); - // Suffix 1: with ghost zones, suffix 0: without ghost zones - const GF3D2layout layout2(cctkGH, indextype); - const GF3D5layout layout5(imin, imax); - - const GF3D2 gf_chi(layout2, chi); - - const smat, 3> gf_gamt{ - GF3D2(layout2, gammatxx), - GF3D2(layout2, gammatxy), - GF3D2(layout2, gammatxz), - GF3D2(layout2, gammatyy), - GF3D2(layout2, gammatyz), - GF3D2(layout2, gammatzz)}; - - const GF3D2 gf_exKh(layout2, Kh); - - const smat, 3> gf_exAt{ - GF3D2(layout2, Atxx), - GF3D2(layout2, Atxy), - GF3D2(layout2, Atxz), - GF3D2(layout2, Atyy), - GF3D2(layout2, Atyz), - GF3D2(layout2, Atzz)}; - - const vec, 3> gf_trGt{ - GF3D2(layout2, Gamtx), - GF3D2(layout2, Gamty), - GF3D2(layout2, Gamtz)}; - - const GF3D2 gf_Theta(layout2, Theta); - - const GF3D2 gf_alpha(layout2, alphaG); - - const vec, 3> gf_beta{ - GF3D2(layout2, betaGx), - GF3D2(layout2, betaGy), - GF3D2(layout2, betaGz)}; - - // Tile variables for derivatives and so on - const int ntmps = 118; - GF3D5vector tmps(layout5, ntmps); - int itmp = 0; - - const auto make_gf = [&]() { return GF3D5(tmps(itmp++)); }; - const auto make_vec = [&](const auto &f) { - return vec, 3>([&](int) { return f(); }); - }; - const auto make_mat = [&](const auto &f) { - return smat, 3>([&](int, int) { return f(); }); - }; - const auto make_vec_gf = [&]() { return make_vec(make_gf); }; - const auto make_mat_gf = [&]() { return make_mat(make_gf); }; - const auto make_vec_vec_gf = [&]() { return make_vec(make_vec_gf); }; - const auto make_mat_vec_gf = [&]() { return make_mat(make_vec_gf); }; - const auto make_mat_mat_gf = [&]() { return make_mat(make_mat_gf); }; - - const GF3D5 tl_chi(make_gf()); - const vec, 3> tl_dchi(make_vec_gf()); - const smat, 3> tl_ddchi(make_mat_gf()); - calc_derivs2(cctkGH, gf_chi, tl_chi, tl_dchi, tl_ddchi, layout5); - - const smat, 3> tl_gamt(make_mat_gf()); - const smat, 3>, 3> tl_dgamt(make_mat_vec_gf()); - const smat, 3>, 3> tl_ddgamt(make_mat_mat_gf()); - calc_derivs2(cctkGH, gf_gamt, tl_gamt, tl_dgamt, tl_ddgamt, layout5); - - const GF3D5 tl_exKh(make_gf()); - const vec, 3> tl_dexKh(make_vec_gf()); - calc_derivs(cctkGH, gf_exKh, tl_exKh, tl_dexKh, layout5); - - const smat, 3> tl_exAt(make_mat_gf()); - const smat, 3>, 3> tl_dexAt(make_mat_vec_gf()); - calc_derivs(cctkGH, gf_exAt, tl_exAt, tl_dexAt, layout5); - - const vec, 3> tl_trGt(make_vec_gf()); - const vec, 3>, 3> tl_dtrGt(make_vec_vec_gf()); - calc_derivs(cctkGH, gf_trGt, tl_trGt, tl_dtrGt, layout5); - - const GF3D5 tl_Theta(make_gf()); - const vec, 3> tl_dTheta(make_vec_gf()); - calc_derivs(cctkGH, gf_Theta, tl_Theta, tl_dTheta, layout5); - - const GF3D5 tl_alpha(make_gf()); - calc_copy(cctkGH, gf_alpha, tl_alpha, layout5); - - const vec, 3> tl_beta(make_vec_gf()); - calc_copy(cctkGH, gf_beta, tl_beta, layout5); - - if (itmp != ntmps) - CCTK_VERROR("Wrong number of temporary variables: ntmps=%d itmp=%d", ntmps, - itmp); - itmp = -1; - - // - - const GF3D2 gf_eTtt(layout2, eTtt); - - const vec, 3> gf_eTt{ - GF3D2(layout2, eTtx), - GF3D2(layout2, eTty), - GF3D2(layout2, eTtz)}; - - const smat, 3> gf_eT{ - GF3D2(layout2, eTxx), - GF3D2(layout2, eTxy), - GF3D2(layout2, eTxz), - GF3D2(layout2, eTyy), - GF3D2(layout2, eTyz), - GF3D2(layout2, eTzz)}; - - // - - const vec, 3> gf_ZtC{GF3D2(layout2, ZtCx), - GF3D2(layout2, ZtCy), - GF3D2(layout2, ZtCz)}; - - const GF3D2 gf_HC(layout2, HC); - - const vec, 3> gf_MtC{GF3D2(layout2, MtCx), - GF3D2(layout2, MtCy), - GF3D2(layout2, MtCz)}; - - // - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const CCTK_REAL cpi = acos(-1.0); - - const Loop::GridDescBaseDevice grid(cctkGH); -#ifdef __CUDACC__ - const nvtxRangeId_t range = nvtxRangeStartA("Z4coo_Constraints::constraints"); -#endif - -#include "../wolfram/Z4coo_set_constraint.hxx" - -#ifdef __CUDACC__ - nvtxRangeEnd(range); -#endif -} - -} // namespace Z4coo diff --git a/Z4coo/src/derivs.hxx b/Z4coo/src/derivs.hxx deleted file mode 100644 index d149769f..00000000 --- a/Z4coo/src/derivs.hxx +++ /dev/null @@ -1,492 +0,0 @@ -#ifndef DERIVS_HXX -#define DERIVS_HXX - -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -namespace Z4coo { -using namespace Arith; -using namespace Loop; -using namespace std; - -//////////////////////////////////////////////////////////////////////////////// - -constexpr int deriv_order = 4; - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv1d(const simdl &mask, const T *restrict const var, const ptrdiff_t di, - const T dx) { - if constexpr (deriv_order == 2) - return -1 / T(2) * - (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[+di])) / dx; - if constexpr (deriv_order == 4) - return (1 / T(12) * - (maskz_loadu(mask, &var[-2 * di]) - - maskz_loadu(mask, &var[+2 * di])) // - - - 2 / T(3) * - (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[+di]))) / - dx; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv1d_upwind(const simdl &mask, const T *restrict const var, - const ptrdiff_t di, const simd &vel, const T dx) { - // arXiv:1111.2177 [gr-qc], (71) - if constexpr (deriv_order == 2) { - // if (sign) - // // + [ 0 -1 +1 0 0] - // // + 1/2 [+1 -2 +1 0 0] - // // [+1/2 -2 +3/2 0 0] - // return (1 / T(2) * var[-2 * di] // - // - 2 * var[-di] // - // + 3 / T(2) * var[0]) / - // dx; - // else - // // + [ 0 0 -1 +1 0 ] - // // - 1/2 [ 0 0 +1 -2 +1 ] - // // [ 0 0 -3/2 +2 -1/2] - // return (-3 / T(2) * var[0] // - // + 2 * var[+di] // - // - 1 / T(2) * var[+2 * di]) / - // dx; - - // + 1/2 [+1/2 -2 +3/2 0 0 ] - // + 1/2 [ 0 0 -3/2 +2 -1/2] - // [+1/4 -1 0 +1 -1/4] - const simd symm = - 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) - - maskz_loadu(mask, &var[2 * di])) // - - (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[di])); - // + 1/2 [+1/2 -2 +3/2 0 0 ] - // - 1/2 [ 0 0 -3/2 +2 -1/2] - // [+1/4 -1 +3/2 -1 +1/4] - const simd anti = - 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[2 * di])) // - - (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[di])) // - + 3 / T(2) * maskz_loadu(mask, &var[0]); - return (vel * symm - fabs(vel) * anti) / dx; - } - if constexpr (deriv_order == 4) { - // A fourth order stencil for a first derivative, shifted by one grid point - - // if (sign) - // return (-1 / T(12) * var[-3 * di] // - // + 1 / T(2) * var[-2 * di] // - // - 3 / T(2) * var[-di] // - // + 5 / T(6) * var[0] // - // + 1 / T(4) * var[+di]) / - // dx; - // else - // return (-1 / T(4) * var[-di] // - // - 5 / T(6) * var[0] // - // + 3 / T(2) * var[+di] // - // - 1 / T(2) * var[+2 * di] // - // + 1 / T(12) * var[+3 * di]) / - // dx; - - const simd symm = - -1 / T(24) * - (maskz_loadu(mask, &var[-3 * di]) - - maskz_loadu(mask, &var[3 * di])) // - + 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) - - maskz_loadu(mask, &var[2 * di])) // - - - 7 / T(8) * (maskz_loadu(mask, &var[-di]) - maskz_loadu(mask, &var[di])); - const simd anti = - -1 / T(24) * - (maskz_loadu(mask, &var[-3 * di]) + - maskz_loadu(mask, &var[3 * di])) // - + 1 / T(4) * - (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[2 * di])) // - - 5 / T(8) * - (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[di])) // - + 5 / T(6) * maskz_loadu(mask, &var[0]); - return (vel * symm - fabs(vel) * anti) / dx; - } -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv2_1d(const simdl &mask, const T *restrict const var, const ptrdiff_t di, - const T dx) { - if constexpr (deriv_order == 2) - return ((maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[+di])) // - - 2 * maskz_loadu(mask, &var[0])) / - pow2(dx); - if constexpr (deriv_order == 4) - return (-1 / T(12) * - (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[+2 * di])) // - + - 4 / T(3) * - (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[+di])) // - - 5 / T(2) * maskz_loadu(mask, &var[0])) / - pow2(dx); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv2_2d(const int vavail, const simdl &mask, const T *restrict const var, - const ptrdiff_t di, const ptrdiff_t dj, const T dx, const T dy) { - constexpr size_t vsize = tuple_size_v>; - if (di == 1) { - assert(vavail > 0); - constexpr int maxnpoints = deriv_order + 1 + vsize - 1; - const int npoints = deriv_order + 1 + min(int(vsize), vavail) - 1; - array, div_ceil(maxnpoints, int(vsize))> arrx; - for (int i = 0; i < maxnpoints; i += vsize) { - if (i < npoints) { - const simdl mask1 = mask_for_loop_tail>(i, npoints); - arrx[div_floor(i, int(vsize))] = - deriv1d(mask1, &var[i - deriv_order / 2], dj, dy); - } - } -#ifdef CCTK_DEBUG - for (int i = npoints; i < align_ceil(maxnpoints, int(vsize)); ++i) - ((T *)&arrx[0])[i] = Arith::nan()(); // unused -#endif - const T *const varx = (T *)&arrx[0] + deriv_order / 2; - return deriv1d(mask, varx, 1, dx); - } else { - assert(dj != 1); - array, deriv_order + 1> arrx; - for (int j = -deriv_order / 2; j <= deriv_order / 2; ++j) - if (j == 0) { -#ifdef CCTK_DEBUG - arrx[deriv_order / 2 + j] = Arith::nan>()(); // unused -#endif - } else { - arrx[deriv_order / 2 + j] = deriv1d(mask, &var[j * dj], di, dx); - } - const T *const varx = (T *)(&arrx[deriv_order / 2]); - return deriv1d(mask, varx, vsize, dy); - } -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv1d_diss(const simdl &mask, const T *restrict const var, - const ptrdiff_t di, const T dx) { - if constexpr (deriv_order == 2) - return ((maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[+2 * di])) // - - - 4 * (maskz_loadu(mask, &var[-di]) + maskz_loadu(mask, &var[+di])) // - + 6 * maskz_loadu(mask, &var[0])) / - dx; - if constexpr (deriv_order == 4) - return ((maskz_loadu(mask, &var[-3 * di]) + - maskz_loadu(mask, &var[+3 * di])) // - - 6 * (maskz_loadu(mask, &var[-2 * di]) + - maskz_loadu(mask, &var[+2 * di])) // - + 15 * (maskz_loadu(mask, &var[-di]) + - maskz_loadu(mask, &var[+di])) // - - 20 * maskz_loadu(mask, &var[0])) / - dx; -} - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - static_assert(dir >= 0 && dir < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir)); - return deriv1d(mask, &gf_(I), di, dx(dir)); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv_upwind(const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec, D> &vel, - const vec &dx) { - static_assert(dir >= 0 && dir < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir)); - return deriv1d_upwind(mask, &gf_(I), di, vel(dir), dx(dir)); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST enable_if_t<(dir1 == dir2), simd> -deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir1 >= 0 && dir1 < D, ""); - static_assert(dir2 >= 0 && dir2 < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir1)); - return deriv2_1d(mask, &gf_(I), di, dx(dir1)); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST enable_if_t<(dir1 != dir2), simd> -deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir1 >= 0 && dir1 < D, ""); - static_assert(dir2 >= 0 && dir2 < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir1)); - const ptrdiff_t dj = gf_.delta(DI(dir2)); - return deriv2_2d(vavail, mask, &gf_(I), di, dj, dx(dir1), dx(dir2)); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv_diss(const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - static_assert(dir >= 0 && dir < D, ""); - const auto &DI = vect::unit; - const ptrdiff_t di = gf_.delta(DI(dir)); - return deriv1d_diss(mask, &gf_(I), di, dx(dir)); -} - -//////////////////////////////////////////////////////////////////////////////// - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST vec, dim> -deriv(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - return {deriv<0>(mask, gf_, I, dx), deriv<1>(mask, gf_, I, dx), - deriv<2>(mask, gf_, I, dx)}; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -deriv_upwind(const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec, dim> &vel, - const vec &dx) { - return deriv_upwind<0>(mask, gf_, I, vel, dx) + - deriv_upwind<1>(mask, gf_, I, vel, dx) + - deriv_upwind<2>(mask, gf_, I, vel, dx); -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST smat, dim> -deriv2(const int vavail, const simdl &mask, const GF3D2 &gf_, - const vect &I, const vec &dx) { - return {deriv2<0, 0>(vavail, mask, gf_, I, dx), - deriv2<0, 1>(vavail, mask, gf_, I, dx), - deriv2<0, 2>(vavail, mask, gf_, I, dx), - deriv2<1, 1>(vavail, mask, gf_, I, dx), - deriv2<1, 2>(vavail, mask, gf_, I, dx), - deriv2<2, 2>(vavail, mask, gf_, I, dx)}; -} - -template -inline ARITH_INLINE ARITH_DEVICE ARITH_HOST simd -diss(const simdl &mask, const GF3D2 &gf_, const vect &I, - const vec &dx) { - // arXiv:gr-qc/0610128, (63), with r=2 - constexpr int diss_order = deriv_order + 2; - constexpr int sign = diss_order % 4 == 0 ? -1 : +1; - return sign / T(pown(2, deriv_order + 2)) * - (deriv_diss<0>(mask, gf_, I, dx) // - + deriv_diss<1>(mask, gf_, I, dx) // - + deriv_diss<2>(mask, gf_, I, dx)); -} - -//////////////////////////////////////////////////////////////////////////////// - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_copy(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - }); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const vec, dim> &dgf0, - const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - const auto dval = deriv(mask, gf1, p.I, dx); - dgf0.store(mask, index0, dval); - }); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs2(const cGH *restrict const cctkGH, const GF3D2 &gf1, - const GF3D5 &gf0, const vec, dim> &dgf0, - const smat, dim> &ddgf0, const GF3D5layout &layout0) { - DECLARE_CCTK_ARGUMENTS; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const int vavail = p.imax - p.i; - const GF3D5index index0(layout0, p.I); - const auto val = gf1(mask, p.I); - gf0.store(mask, index0, val); - const auto dval = deriv(mask, gf1, p.I, dx); - dgf0.store(mask, index0, dval); - const auto ddval = deriv2(vavail, mask, gf1, p.I, dx); - ddgf0.store(mask, index0, ddval); - }); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_copy(const cGH *restrict const cctkGH, - const vec, dim> &gf0_, - const vec, dim> &gf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_copy(cctkGH, gf0_(a), gf_(a), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -calc_derivs(const cGH *restrict const cctkGH, - const vec, dim> &gf0_, const vec, dim> &gf_, - const vec, dim>, dim> &dgf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_derivs(cctkGH, gf0_(a), gf_(a), dgf_(a), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_derivs2( - const cGH *restrict const cctkGH, const vec, dim> &gf0_, - const vec, dim> &gf_, const vec, dim>, dim> &dgf_, - const vec, dim>, dim> &ddgf_, const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - calc_derivs2(cctkGH, gf0_(a), gf_(a), dgf_(a), ddgf_(a), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_copy(const cGH *restrict const cctkGH, - const smat, dim> &gf0_, - const smat, dim> &gf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_copy(cctkGH, gf0_(a, b), gf_(a, b), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_derivs( - const cGH *restrict const cctkGH, const smat, dim> &gf0_, - const smat, dim> &gf_, const smat, dim>, dim> &dgf_, - const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_derivs(cctkGH, gf0_(a, b), gf_(a, b), dgf_(a, b), layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void calc_derivs2( - const cGH *restrict const cctkGH, const smat, dim> &gf0_, - const smat, dim> &gf_, const smat, dim>, dim> &dgf_, - const smat, dim>, dim> &ddgf_, const GF3D5layout &layout) { - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - calc_derivs2(cctkGH, gf0_(a, b), gf_(a, b), dgf_(a, b), ddgf_(a, b), - layout); -} - -template -CCTK_ATTRIBUTE_NOINLINE void -apply_upwind_diss(const cGH *restrict const cctkGH, const GF3D2 &gf_, - const vec, dim> &gf_betaG_, - const GF3D2 &gf_rhs_) { - DECLARE_CCTK_ARGUMENTS; - DECLARE_CCTK_PARAMETERS; - - const vec dx([&](int a) { return CCTK_DELTA_SPACE(a); }); - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - if (epsdiss == 0) { - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const vec betaG = gf_betaG_(mask, p.I); - const vreal rhs_old = gf_rhs_(mask, p.I); - const vreal rhs_new = - rhs_old + deriv_upwind(mask, gf_, p.I, betaG, dx); - gf_rhs_.store(mask, p.I, rhs_new); - }); - - } else { - - const Loop::GridDescBaseDevice grid(cctkGH); - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const vec betaG = gf_betaG_(mask, p.I); - const vreal rhs_old = gf_rhs_(mask, p.I); - const vreal rhs_new = rhs_old + - deriv_upwind(mask, gf_, p.I, betaG, dx) + - epsdiss * diss(mask, gf_, p.I, dx); - gf_rhs_.store(mask, p.I, rhs_new); - }); - } -} - -} // namespace Z4coo - -#endif // #ifndef DERIVS_HXX diff --git a/Z4coo/src/enforce.cxx b/Z4coo/src/enforce.cxx deleted file mode 100644 index c43cc157..00000000 --- a/Z4coo/src/enforce.cxx +++ /dev/null @@ -1,128 +0,0 @@ -#include -#include -#include -#include - -#include -#include -#include - -#ifdef __CUDACC__ -#include -#endif - -#include -#include - -namespace Z4coo { -using namespace Arith; -using namespace Loop; -using namespace std; - -extern "C" void Z4coo_Enforce(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_Z4coo_Enforce; - DECLARE_CCTK_PARAMETERS; - - const array indextype = {0, 0, 0}; - const GF3D2layout layout1(cctkGH, indextype); - - const GF3D2 &gf_chi = chi; - - const GF3D2 &gf_alphaG = alphaG; - - const smat, 3> gf_gammat{ - gammatxx, gammatxy, gammatxz, gammatyy, gammatyz, gammatzz, - }; - - const smat, 3> gf_At{ - Atxx, Atxy, Atxz, Atyy, Atyz, Atzz, - }; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const auto delta3 = one >()(); - -#ifdef __CUDACC__ - const nvtxRangeId_t range = nvtxRangeStartA("Z4coo_Enforce::enforce"); -#endif - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D2index index1(layout1, p.I); - - // Load - const vreal chi_old = gf_chi(mask, index1); - const vreal alphaG_old = gf_alphaG(mask, index1); - - const smat gammat_old = gf_gammat(mask, index1); - const smat At_old = gf_At(mask, index1); - - // Enforce floors - - const vreal chi = fmax(vreal(chi_floor - 1), chi_old); - const vreal alphaG = fmax(vreal(alphaG_floor - 1), alphaG_old); - - // Enforce algebraic constraints - // See arXiv:1212.2901 [gr-qc]. - - const vreal detgammat_old = calc_det(delta3 + gammat_old); - const vreal chi1_old = 1 / cbrt(detgammat_old) - 1; - const smat gammat([&](int a, int b) ARITH_INLINE { - return (1 + chi1_old) * (delta3(a, b) + gammat_old(a, b)) - - delta3(a, b); - }); -#ifdef CCTK_DEBUG - const vreal detgammat = calc_det(delta3 + gammat); - const vreal gammat_norm = maxabs(delta3 + gammat); - const vreal gammat_scale = gammat_norm; -#if !defined __CUDACC__ && !defined __HIPCC__ - if (!(all(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale))) { - ostringstream buf; - buf << "det gammat is not one: gammat=" << gammat - << " det(gammat)=" << detgammat; - CCTK_VERROR("%s", buf.str().c_str()); - } -#endif - assert(all(fabs(detgammat - 1) <= 1.0e-12 * gammat_scale)); -#endif - - const smat gammatu = - calc_inv(delta3 + gammat, vreal(1)) - delta3; - - const vreal traceAt_old = sum_symm<3>([&](int x, int y) ARITH_INLINE { - return (delta3(x, y) + gammatu(x, y)) * At_old(x, y); - }); - const smat At([&](int a, int b) ARITH_INLINE { - return At_old(a, b) - traceAt_old / 3 * (delta3(a, b) + gammat(a, b)); - }); -#ifdef CCTK_DEBUG - const vreal traceAt = sum_symm<3>([&](int x, int y) ARITH_INLINE { - return (delta3(x, y) + gammatu(x, y)) * At(x, y); - }); - const vreal gammatu_norm = maxabs(delta3 + gammatu); - const vreal At_norm = maxabs(At); - const vreal At_scale = fmax(fmax(gammat_norm, gammatu_norm), At_norm); -#if !defined __CUDACC__ && !defined __HIPCC__ - if (!(all(fabs(traceAt) <= 1.0e-12 * At_scale))) { - ostringstream buf; - buf << "tr At: At=" << At << " tr(At)=" << traceAt; - CCTK_VERROR("%s", buf.str().c_str()); - } -#endif - assert(all(fabs(traceAt) <= 1.0e-12 * At_scale)); -#endif - - // Store - gf_chi.store(mask, index1, chi); - gf_gammat.store(mask, index1, gammat); - gf_At.store(mask, index1, At); - gf_alphaG.store(mask, index1, alphaG); - }); -#ifdef __CUDACC__ - nvtxRangeEnd(range); -#endif -} - -} // namespace Z4coo diff --git a/Z4coo/src/initial1.cxx b/Z4coo/src/initial1.cxx deleted file mode 100644 index 3160f630..00000000 --- a/Z4coo/src/initial1.cxx +++ /dev/null @@ -1,131 +0,0 @@ -#include -#include -#include -#include - -#include -#include - -#ifdef __CUDACC__ -#include -#endif - -#include - -namespace Z4coo { -using namespace Arith; -using namespace Loop; -using namespace std; - -extern "C" void Z4coo_Initial1(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS_Z4coo_Initial1; - - const array indextype = {0, 0, 0}; - const GF3D2layout layout1(cctkGH, indextype); - - const smat, 3> gf_g1{ - GF3D2(layout1, gxx), - GF3D2(layout1, gxy), - GF3D2(layout1, gxz), - GF3D2(layout1, gyy), - GF3D2(layout1, gyz), - GF3D2(layout1, gzz)}; - - const smat, 3> gf_K1{ - GF3D2(layout1, kxx), - GF3D2(layout1, kxy), - GF3D2(layout1, kxz), - GF3D2(layout1, kyy), - GF3D2(layout1, kyz), - GF3D2(layout1, kzz)}; - - const GF3D2 gf_alp1(layout1, alp); - - const vec, 3> gf_beta1{ - GF3D2(layout1, betax), - GF3D2(layout1, betay), - GF3D2(layout1, betaz)}; - - const GF3D2 gf_chi1(layout1, chi); - - const smat, 3> gf_gammat1{ - GF3D2(layout1, gammatxx), GF3D2(layout1, gammatxy), - GF3D2(layout1, gammatxz), GF3D2(layout1, gammatyy), - GF3D2(layout1, gammatyz), GF3D2(layout1, gammatzz)}; - - const GF3D2 gf_Kh1(layout1, Kh); - - const smat, 3> gf_At1{ - GF3D2(layout1, Atxx), GF3D2(layout1, Atxy), - GF3D2(layout1, Atxz), GF3D2(layout1, Atyy), - GF3D2(layout1, Atyz), GF3D2(layout1, Atzz)}; - - const GF3D2 gf_Theta1(layout1, Theta); - - const GF3D2 gf_alphaG1(layout1, alphaG); - - const vec, 3> gf_betaG1{GF3D2(layout1, betaGx), - GF3D2(layout1, betaGy), - GF3D2(layout1, betaGz)}; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const auto delta3 = one >()(); - - const Loop::GridDescBaseDevice grid(cctkGH); -#ifdef __CUDACC__ - const nvtxRangeId_t range = nvtxRangeStartA("Z4coo_Initial1::initial1"); -#endif - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D2index index1(layout1, p.I); - - // Load - const smat g = gf_g1(mask, index1, one >()()); - const smat K = gf_K1(mask, index1); - const vreal alp = gf_alp1(mask, index1, 1); - const vec beta = gf_beta1(mask, index1); - - // Calculate Z4c variables (all except Gammat) - const vreal detg = calc_det(g); - const smat gu = calc_inv(g, detg); - - const vreal chi = 1 / cbrt(detg) - 1; - - const smat gammat([&](int a, int b) ARITH_INLINE { - return (1 + chi) * g(a, b) - delta3(a, b); - }); - - const vreal trK = sum_symm<3>( - [&](int x, int y) ARITH_INLINE { return gu(x, y) * K(x, y); }); - - const vreal Theta = 0; - - const vreal Kh = trK - 2 * Theta; - - const smat At([&](int a, int b) ARITH_INLINE { - return (1 + chi) * (K(a, b) - trK / 3 * g(a, b)); - }); - - const vreal alphaG = alp - 1; - - const vec betaG([&](int a) ARITH_INLINE { return beta(a); }); - - // Store - gf_chi1.store(mask, index1, chi); - gf_gammat1.store(mask, index1, gammat); - gf_Kh1.store(mask, index1, Kh); - gf_At1.store(mask, index1, At); - gf_Theta1.store(mask, index1, Theta); - gf_alphaG1.store(mask, index1, alphaG); - gf_betaG1.store(mask, index1, betaG); - }); -#ifdef __CUDACC__ - nvtxRangeEnd(range); -#endif -} - -} // namespace Z4coo diff --git a/Z4coo/src/initial2.cxx b/Z4coo/src/initial2.cxx deleted file mode 100644 index 141724cc..00000000 --- a/Z4coo/src/initial2.cxx +++ /dev/null @@ -1,108 +0,0 @@ -#include "derivs.hxx" - -#include -#include - -#include -#include - -#ifdef __CUDACC__ -#include -#endif - -namespace Z4coo { -using namespace Arith; -using namespace Loop; -using namespace std; - -template -ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr vec, D> -calc_gammal(const smat, D> &dg) { - // Gammal_abc - return vec, D>([&](int a) ARITH_INLINE { - return smat([&](int b, int c) ARITH_INLINE { - return (dg(a, b)(c) + dg(a, c)(b) - dg(b, c)(a)) / 2; - }); - }); -} - -template -ARITH_INLINE ARITH_DEVICE ARITH_HOST constexpr vec, D> -calc_gamma(const smat &gu, const vec, D> &Gammal) { - // Gamma^a_bc - return vec, D>([&](int a) ARITH_INLINE { - return smat([&](int b, int c) ARITH_INLINE { - return sum([&](int x) - ARITH_INLINE { return gu(a, x) * Gammal(x)(b, c); }); - }); - }); -} - -extern "C" void Z4coo_Initial2(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS_Z4coo_Initial2; - - const vec dx{ - CCTK_DELTA_SPACE(0), - CCTK_DELTA_SPACE(1), - CCTK_DELTA_SPACE(2), - }; - - const array indextype = {0, 0, 0}; - const GF3D2layout layout1(cctkGH, indextype); - - const smat, 3> gf_gammat1{ - GF3D2(layout1, gammatxx), - GF3D2(layout1, gammatxy), - GF3D2(layout1, gammatxz), - GF3D2(layout1, gammatyy), - GF3D2(layout1, gammatyz), - GF3D2(layout1, gammatzz)}; - - const vec, 3> gf_Gamt1{GF3D2(layout1, Gamtx), - GF3D2(layout1, Gamty), - GF3D2(layout1, Gamtz)}; - - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - const auto delta3 = one >()(); - - const Loop::GridDescBaseDevice grid(cctkGH); -#ifdef __CUDACC__ - const nvtxRangeId_t range = nvtxRangeStartA("Z4coo_Initial2::initial2"); -#endif - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D2index index1(layout1, p.I); - - // Load - const smat gammat = gf_gammat1(mask, index1); - - // Calculate Z4c variables (only Gamt) - const smat gammatu = - calc_inv(delta3 + gammat, vreal(1)) - delta3; - - const smat, 3> dgammat([&](int a, int b) { - return deriv(mask, gf_gammat1(a, b), p.I, dx); - }); - - const vec, 3> Gammatl = calc_gammal(dgammat); - const vec, 3> Gammat = - calc_gamma(delta3 + gammatu, Gammatl); - const vec Gamt([&](int a) ARITH_INLINE { - return sum_symm<3>([&](int x, int y) ARITH_INLINE { - return (delta3(x, y) + gammatu(x, y)) * Gammat(a)(x, y); - }); - }); - - // Store - gf_Gamt1.store(mask, index1, Gamt); - }); -#ifdef __CUDACC__ - nvtxRangeEnd(range); -#endif -} - -} // namespace Z4coo diff --git a/Z4coo/src/make.code.defn b/Z4coo/src/make.code.defn deleted file mode 100644 index da24e9d2..00000000 --- a/Z4coo/src/make.code.defn +++ /dev/null @@ -1,14 +0,0 @@ -# Main make.code.defn file for thorn Z4c - -# Source files in this directory -SRCS = \ - adm.cxx \ - constraint.cxx \ - enforce.cxx \ - initial1.cxx \ - initial2.cxx \ - newradx.cxx \ - rhs.cxx - -# Subdirectories containing source files -SUBDIRS = diff --git a/Z4coo/src/newradx.cxx b/Z4coo/src/newradx.cxx deleted file mode 100644 index 151f4ea7..00000000 --- a/Z4coo/src/newradx.cxx +++ /dev/null @@ -1,36 +0,0 @@ -#include -#include -#include -#include - -namespace Z4coo { - -extern "C" void Z4coo_Apply_NewRadX_BC(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTSX_Z4coo_Apply_NewRadX_BC; - DECLARE_CCTK_PARAMETERS; - - NewRadX::NewRadX_Apply(cctkGH, chi, chi_rhs, 1, 1, n_chi); - NewRadX::NewRadX_Apply(cctkGH, gammatxx, gammatxx_rhs, 1, 1, n_gammat); - NewRadX::NewRadX_Apply(cctkGH, gammatxy, gammatxy_rhs, 0, 1, n_gammat); - NewRadX::NewRadX_Apply(cctkGH, gammatxz, gammatxz_rhs, 0, 1, n_gammat); - NewRadX::NewRadX_Apply(cctkGH, gammatyy, gammatyy_rhs, 1, 1, n_gammat); - NewRadX::NewRadX_Apply(cctkGH, gammatyz, gammatyz_rhs, 0, 1, n_gammat); - NewRadX::NewRadX_Apply(cctkGH, gammatzz, gammatzz_rhs, 1, 1, n_gammat); - NewRadX::NewRadX_Apply(cctkGH, Kh, Kh_rhs, 0, 1, n_Kh); - NewRadX::NewRadX_Apply(cctkGH, Atxx, Atxx_rhs, 0, 1, n_At); - NewRadX::NewRadX_Apply(cctkGH, Atxy, Atxy_rhs, 0, 1, n_At); - NewRadX::NewRadX_Apply(cctkGH, Atxz, Atxz_rhs, 0, 1, n_At); - NewRadX::NewRadX_Apply(cctkGH, Atyy, Atyy_rhs, 0, 1, n_At); - NewRadX::NewRadX_Apply(cctkGH, Atyz, Atyz_rhs, 0, 1, n_At); - NewRadX::NewRadX_Apply(cctkGH, Atzz, Atzz_rhs, 0, 1, n_At); - NewRadX::NewRadX_Apply(cctkGH, Gamtx, Gamtx_rhs, 0, 1, n_Gamt); - NewRadX::NewRadX_Apply(cctkGH, Gamty, Gamty_rhs, 0, 1, n_Gamt); - NewRadX::NewRadX_Apply(cctkGH, Gamtz, Gamtz_rhs, 0, 1, n_Gamt); - NewRadX::NewRadX_Apply(cctkGH, Theta, Theta_rhs, 0, 1, n_Theta); - NewRadX::NewRadX_Apply(cctkGH, alphaG, alphaG_rhs, 1, 1, n_alphaG); - NewRadX::NewRadX_Apply(cctkGH, betaGx, betaGx_rhs, 0, 1, n_betaG); - NewRadX::NewRadX_Apply(cctkGH, betaGy, betaGy_rhs, 0, 1, n_betaG); - NewRadX::NewRadX_Apply(cctkGH, betaGz, betaGz_rhs, 0, 1, n_betaG); -} - -} // namespace Z4coo diff --git a/Z4coo/src/rhs.cxx b/Z4coo/src/rhs.cxx deleted file mode 100644 index bd466fac..00000000 --- a/Z4coo/src/rhs.cxx +++ /dev/null @@ -1,253 +0,0 @@ -#include - -#ifdef __CUDACC__ -// Disable CCTK_DEBUG since the debug information takes too much -// parameter space to launch the kernels -#ifdef CCTK_DEBUG -#undef CCTK_DEBUG -#endif -#endif - -// #define Power(x, y) (Arith::pown((x), (y))) - -#include "derivs.hxx" - -#include -#include -#include -#include - -#include -#include -#include - -#ifdef __CUDACC__ -#include -#endif - -#include - -namespace Z4coo { -using namespace Arith; -using namespace Loop; -using namespace std; - -template inline T Power(T x, int y) { - return (y == 2) ? Arith::pow2(x) : Arith::pown(x, y); -} - -extern "C" void Z4coo_RHS(CCTK_ARGUMENTS) { - DECLARE_CCTK_ARGUMENTS_Z4coo_RHS; - DECLARE_CCTK_PARAMETERS; - - for (int d = 0; d < 3; ++d) - if (cctk_nghostzones[d] < deriv_order / 2 + 1) - CCTK_VERROR("Need at least %d ghost zones", deriv_order / 2 + 1); - - const array indextype = {0, 0, 0}; - const array nghostzones = {cctk_nghostzones[0], cctk_nghostzones[1], - cctk_nghostzones[2]}; - vect imin, imax; - GridDescBase(cctkGH).box_int<0, 0, 0>(nghostzones, imin, imax); - - // Suffix 1: with ghost zones, suffix 0: without ghost zones - const GF3D2layout layout2(cctkGH, indextype); - const GF3D5layout layout5(imin, imax); - - // Input grid functions - const GF3D2 gf_chi(layout2, chi); - - const smat, 3> gf_gamt{ - GF3D2(layout2, gammatxx), - GF3D2(layout2, gammatxy), - GF3D2(layout2, gammatxz), - GF3D2(layout2, gammatyy), - GF3D2(layout2, gammatyz), - GF3D2(layout2, gammatzz)}; - - const GF3D2 gf_exKh(layout2, Kh); - - const smat, 3> gf_exAt{ - GF3D2(layout2, Atxx), - GF3D2(layout2, Atxy), - GF3D2(layout2, Atxz), - GF3D2(layout2, Atyy), - GF3D2(layout2, Atyz), - GF3D2(layout2, Atzz)}; - - const vec, 3> gf_trGt{ - GF3D2(layout2, Gamtx), - GF3D2(layout2, Gamty), - GF3D2(layout2, Gamtz)}; - - const GF3D2 gf_Theta(layout2, Theta); - - const GF3D2 gf_alpha(layout2, alphaG); - - const vec, 3> gf_beta{ - GF3D2(layout2, betaGx), - GF3D2(layout2, betaGy), - GF3D2(layout2, betaGz)}; - - // Tile variables for derivatives and so on - const int ntmps = 136; - GF3D5vector tmps(layout5, ntmps); - int itmp = 0; - - const auto make_gf = [&]() { return GF3D5(tmps(itmp++)); }; - const auto make_vec = [&](const auto &f) { - return vec, 3>([&](int) { return f(); }); - }; - const auto make_mat = [&](const auto &f) { - return smat, 3>([&](int, int) { return f(); }); - }; - const auto make_vec_gf = [&]() { return make_vec(make_gf); }; - const auto make_mat_gf = [&]() { return make_mat(make_gf); }; - const auto make_vec_vec_gf = [&]() { return make_vec(make_vec_gf); }; - const auto make_vec_mat_gf = [&]() { return make_vec(make_mat_gf); }; - const auto make_mat_vec_gf = [&]() { return make_mat(make_vec_gf); }; - const auto make_mat_mat_gf = [&]() { return make_mat(make_mat_gf); }; - - const GF3D5 tl_chi(make_gf()); - const vec, 3> tl_dchi(make_vec_gf()); - const smat, 3> tl_ddchi(make_mat_gf()); - calc_derivs2(cctkGH, gf_chi, tl_chi, tl_dchi, tl_ddchi, layout5); - - const smat, 3> tl_gamt(make_mat_gf()); - const smat, 3>, 3> tl_dgamt(make_mat_vec_gf()); - const smat, 3>, 3> tl_ddgamt(make_mat_mat_gf()); - calc_derivs2(cctkGH, gf_gamt, tl_gamt, tl_dgamt, tl_ddgamt, layout5); - - const GF3D5 tl_exKh(make_gf()); - const vec, 3> tl_dexKh(make_vec_gf()); - calc_derivs(cctkGH, gf_exKh, tl_exKh, tl_dexKh, layout5); - - const smat, 3> tl_exAt(make_mat_gf()); - calc_copy(cctkGH, gf_exAt, tl_exAt, layout5); - - const vec, 3> tl_trGt(make_vec_gf()); - const vec, 3>, 3> tl_dtrGt(make_vec_vec_gf()); - calc_derivs(cctkGH, gf_trGt, tl_trGt, tl_dtrGt, layout5); - - const GF3D5 tl_Theta(make_gf()); - const vec, 3> tl_dTheta(make_vec_gf()); - calc_derivs(cctkGH, gf_Theta, tl_Theta, tl_dTheta, layout5); - - const GF3D5 tl_alpha(make_gf()); - const vec, 3> tl_dalpha(make_vec_gf()); - const smat, 3> tl_ddalpha(make_mat_gf()); - calc_derivs2(cctkGH, gf_alpha, tl_alpha, tl_dalpha, tl_ddalpha, layout5); - - const vec, 3> tl_beta(make_vec_gf()); - const vec, 3>, 3> tl_dbeta(make_vec_vec_gf()); - const vec, 3>, 3> tl_ddbeta(make_vec_mat_gf()); - calc_derivs2(cctkGH, gf_beta, tl_beta, tl_dbeta, tl_ddbeta, layout5); - - if (itmp != ntmps) - CCTK_VERROR("Wrong number of temporary variables: ntmps=%d itmp=%d", ntmps, - itmp); - itmp = -1; - - // More input grid functions - const GF3D2 gf_eTtt(layout2, eTtt); - - const vec, 3> gf_eTt{ - GF3D2(layout2, eTtx), - GF3D2(layout2, eTty), - GF3D2(layout2, eTtz)}; - - const smat, 3> gf_eT{ - GF3D2(layout2, eTxx), - GF3D2(layout2, eTxy), - GF3D2(layout2, eTxz), - GF3D2(layout2, eTyy), - GF3D2(layout2, eTyz), - GF3D2(layout2, eTzz)}; - - // Output grid functions - const GF3D2 gf_dtchi(layout2, chi_rhs); - - const smat, 3> gf_dtgamt{ - GF3D2(layout2, gammatxx_rhs), - GF3D2(layout2, gammatxy_rhs), - GF3D2(layout2, gammatxz_rhs), - GF3D2(layout2, gammatyy_rhs), - GF3D2(layout2, gammatyz_rhs), - GF3D2(layout2, gammatzz_rhs)}; - - const GF3D2 gf_dtexKh(layout2, Kh_rhs); - - const smat, 3> gf_dtexAt{ - GF3D2(layout2, Atxx_rhs), GF3D2(layout2, Atxy_rhs), - GF3D2(layout2, Atxz_rhs), GF3D2(layout2, Atyy_rhs), - GF3D2(layout2, Atyz_rhs), GF3D2(layout2, Atzz_rhs)}; - - const vec, 3> gf_dttrGt{ - GF3D2(layout2, Gamtx_rhs), - GF3D2(layout2, Gamty_rhs), - GF3D2(layout2, Gamtz_rhs)}; - - const GF3D2 gf_dtTheta(layout2, Theta_rhs); - - const GF3D2 gf_dtalpha(layout2, alphaG_rhs); - - const vec, 3> gf_dtbeta{ - GF3D2(layout2, betaGx_rhs), - GF3D2(layout2, betaGy_rhs), - GF3D2(layout2, betaGz_rhs)}; - - // simd types - typedef simd vreal; - typedef simdl vbool; - constexpr size_t vsize = tuple_size_v; - - // parameters - const CCTK_REAL cpi = acos(-1.0); - const CCTK_REAL ckappa1 = kappa1; - const CCTK_REAL ckappa2 = kappa2; - const CCTK_REAL cmuL = f_mu_L; - const CCTK_REAL cmuS = f_mu_S; - const CCTK_REAL ceta = eta; - - // Loop - const Loop::GridDescBaseDevice grid(cctkGH); - -#ifdef __CUDACC__ - const nvtxRangeId_t range = nvtxRangeStartA("Z4coo_RHS::rhs"); -#endif - -#include "../wolfram/Z4coo_set_rhs.hxx" - -#ifdef __CUDACC__ - nvtxRangeEnd(range); -#endif - - // Upwind and dissipation terms - - // TODO: Consider fusing the loops to reduce memory bandwidth - - apply_upwind_diss(cctkGH, gf_chi, gf_beta, gf_dtchi); - - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - apply_upwind_diss(cctkGH, gf_gamt(a, b), gf_beta, gf_dtgamt(a, b)); - - apply_upwind_diss(cctkGH, gf_exKh, gf_beta, gf_dtexKh); - - for (int a = 0; a < 3; ++a) - for (int b = a; b < 3; ++b) - apply_upwind_diss(cctkGH, gf_exAt(a, b), gf_beta, gf_dtexAt(a, b)); - - for (int a = 0; a < 3; ++a) - apply_upwind_diss(cctkGH, gf_trGt(a), gf_beta, gf_dttrGt(a)); - - if (!set_Theta_zero) - apply_upwind_diss(cctkGH, gf_Theta, gf_beta, gf_dtTheta); - - apply_upwind_diss(cctkGH, gf_alpha, gf_beta, gf_dtalpha); - - for (int a = 0; a < 3; ++a) - apply_upwind_diss(cctkGH, gf_beta(a), gf_beta, gf_dtbeta(a)); -} - -} // namespace Z4coo diff --git a/Z4coo/wolfram/Z4coo_set_ADM.hxx b/Z4coo/wolfram/Z4coo_set_ADM.hxx deleted file mode 100644 index 1a5ad5cd..00000000 --- a/Z4coo/wolfram/Z4coo_set_ADM.hxx +++ /dev/null @@ -1,126 +0,0 @@ -/* Z4coo_set_ADM.hxx */ -/* Produced with Mathematica */ - -#ifndef Z4COO_SET_ADM_HXX -#define Z4COO_SET_ADM_HXX - -const GF3D2 &local_ADMgam11 = gf_ADMgam(0,0); -const GF3D2 &local_ADMgam12 = gf_ADMgam(0,1); -const GF3D2 &local_ADMgam13 = gf_ADMgam(0,2); -const GF3D2 &local_ADMgam22 = gf_ADMgam(1,1); -const GF3D2 &local_ADMgam23 = gf_ADMgam(1,2); -const GF3D2 &local_ADMgam33 = gf_ADMgam(2,2); -const GF3D2 &local_ADMK11 = gf_ADMK(0,0); -const GF3D2 &local_ADMK12 = gf_ADMK(0,1); -const GF3D2 &local_ADMK13 = gf_ADMK(0,2); -const GF3D2 &local_ADMK22 = gf_ADMK(1,1); -const GF3D2 &local_ADMK23 = gf_ADMK(1,2); -const GF3D2 &local_ADMK33 = gf_ADMK(2,2); -const GF3D2 &local_ADMalpha = gf_ADMalpha; -const GF3D2 &local_ADMbeta1 = gf_ADMbeta(0); -const GF3D2 &local_ADMbeta2 = gf_ADMbeta(1); -const GF3D2 &local_ADMbeta3 = gf_ADMbeta(2); - -grid.loop_all_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D2index index2(layout2, p.I); - -const auto &tmp_chi = gf_chi(mask, index2); -const auto &tmp_gamt = gf_gamt(mask, index2); -const auto &tmp_exKh = gf_exKh(mask, index2); -const auto &tmp_exAt = gf_exAt(mask, index2); -const auto &tmp_Theta = gf_Theta(mask, index2); -const auto &tmp_alpha = gf_alpha(mask, index2); -const auto &tmp_beta = gf_beta(mask, index2); - -const vreal chi = tmp_chi + 1; -const vreal gamt11 = tmp_gamt(0,0) + 1; -const vreal gamt12 = tmp_gamt(0,1); -const vreal gamt13 = tmp_gamt(0,2); -const vreal gamt22 = tmp_gamt(1,1) + 1; -const vreal gamt23 = tmp_gamt(1,2); -const vreal gamt33 = tmp_gamt(2,2) + 1; -const vreal exKh = tmp_exKh; -const vreal exAt11 = tmp_exAt(0,0); -const vreal exAt12 = tmp_exAt(0,1); -const vreal exAt13 = tmp_exAt(0,2); -const vreal exAt22 = tmp_exAt(1,1); -const vreal exAt23 = tmp_exAt(1,2); -const vreal exAt33 = tmp_exAt(2,2); -const vreal Theta = tmp_Theta; -const vreal alpha = tmp_alpha + 1; -const vreal beta1 = tmp_beta(0); -const vreal beta2 = tmp_beta(1); -const vreal beta3 = tmp_beta(2); - -local_ADMgam11.store(mask, index2, -gamt11/chi -); - -local_ADMgam12.store(mask, index2, -gamt12/chi -); - -local_ADMgam13.store(mask, index2, -gamt13/chi -); - -local_ADMgam22.store(mask, index2, -gamt22/chi -); - -local_ADMgam23.store(mask, index2, -gamt23/chi -); - -local_ADMgam33.store(mask, index2, -gamt33/chi -); - -local_ADMK11.store(mask, index2, -(exAt11 + (gamt11*(exKh + 2*Theta))/3.)/chi -); - -local_ADMK12.store(mask, index2, -(exAt12 + (gamt12*(exKh + 2*Theta))/3.)/chi -); - -local_ADMK13.store(mask, index2, -(exAt13 + (gamt13*(exKh + 2*Theta))/3.)/chi -); - -local_ADMK22.store(mask, index2, -(exAt22 + (gamt22*(exKh + 2*Theta))/3.)/chi -); - -local_ADMK23.store(mask, index2, -(exAt23 + (gamt23*(exKh + 2*Theta))/3.)/chi -); - -local_ADMK33.store(mask, index2, -(exAt33 + (gamt33*(exKh + 2*Theta))/3.)/chi -); - -local_ADMalpha.store(mask, index2, -alpha -); - -local_ADMbeta1.store(mask, index2, -beta1 -); - -local_ADMbeta2.store(mask, index2, -beta2 -); - -local_ADMbeta3.store(mask, index2, -beta3 -); - - -}); - -#endif // #ifndef Z4COO_SET_ADM_HXX - -/* Z4coo_set_ADM.hxx */ diff --git a/Z4coo/wolfram/Z4coo_set_constraint.hxx b/Z4coo/wolfram/Z4coo_set_constraint.hxx deleted file mode 100644 index d55de695..00000000 --- a/Z4coo/wolfram/Z4coo_set_constraint.hxx +++ /dev/null @@ -1,1034 +0,0 @@ -/* Z4coo_set_constraint.hxx */ -/* Produced with Mathematica */ - -#ifndef Z4COO_SET_CONSTRAINT_HXX -#define Z4COO_SET_CONSTRAINT_HXX - -const GF3D2 &local_ZtC1 = gf_ZtC(0); -const GF3D2 &local_ZtC2 = gf_ZtC(1); -const GF3D2 &local_ZtC3 = gf_ZtC(2); -const GF3D2 &local_HC = gf_HC; -const GF3D2 &local_MtC1 = gf_MtC(0); -const GF3D2 &local_MtC2 = gf_MtC(1); -const GF3D2 &local_MtC3 = gf_MtC(2); - -noinline([&]() __attribute__((__flatten__, __hot__)) { - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D2index index2(layout2, p.I); - const GF3D5index index5(layout5, p.I); - -const auto &tmp_eTtt = gf_eTtt(mask, index2); -const auto &tmp_eTt = gf_eTt(mask, index2); -const auto &tmp_eT = gf_eT(mask, index2); -const auto &tmp_chi = tl_chi(mask, index5); -const auto &tmp_gamt = tl_gamt(mask, index5); -const auto &tmp_exKh = tl_exKh(mask, index5); -const auto &tmp_exAt = tl_exAt(mask, index5); -const auto &tmp_trGt = tl_trGt(mask, index5); -const auto &tmp_Theta = tl_Theta(mask, index5); -const auto &tmp_alpha = tl_alpha(mask, index5); -const auto &tmp_beta = tl_beta(mask, index5); -const auto &tmp_dchi = tl_dchi(mask, index5); -const auto &tmp_dgamt = tl_dgamt(mask, index5); -const auto &tmp_dexKh = tl_dexKh(mask, index5); -const auto &tmp_dexAt = tl_dexAt(mask, index5); -const auto &tmp_dtrGt = tl_dtrGt(mask, index5); -const auto &tmp_dTheta = tl_dTheta(mask, index5); -const auto &tmp_ddchi = tl_ddchi(mask, index5); -const auto &tmp_ddgamt = tl_ddgamt(mask, index5); - -const vreal eTtt = tmp_eTtt; -const vreal eTt1 = tmp_eTt(0); -const vreal eTt2 = tmp_eTt(1); -const vreal eTt3 = tmp_eTt(2); -const vreal eT11 = tmp_eT(0,0); -const vreal eT12 = tmp_eT(0,1); -const vreal eT13 = tmp_eT(0,2); -const vreal eT22 = tmp_eT(1,1); -const vreal eT23 = tmp_eT(1,2); -const vreal eT33 = tmp_eT(2,2); -const vreal chi = tmp_chi + 1; -const vreal gamt11 = tmp_gamt(0,0) + 1; -const vreal gamt12 = tmp_gamt(0,1); -const vreal gamt13 = tmp_gamt(0,2); -const vreal gamt22 = tmp_gamt(1,1) + 1; -const vreal gamt23 = tmp_gamt(1,2); -const vreal gamt33 = tmp_gamt(2,2) + 1; -const vreal exKh = tmp_exKh; -const vreal exAt11 = tmp_exAt(0,0); -const vreal exAt12 = tmp_exAt(0,1); -const vreal exAt13 = tmp_exAt(0,2); -const vreal exAt22 = tmp_exAt(1,1); -const vreal exAt23 = tmp_exAt(1,2); -const vreal exAt33 = tmp_exAt(2,2); -const vreal trGt1 = tmp_trGt(0); -const vreal trGt2 = tmp_trGt(1); -const vreal trGt3 = tmp_trGt(2); -const vreal Theta = tmp_Theta; -const vreal alpha = tmp_alpha + 1; -const vreal beta1 = tmp_beta(0); -const vreal beta2 = tmp_beta(1); -const vreal beta3 = tmp_beta(2); -const vreal dchi1 = tmp_dchi(0); -const vreal dchi2 = tmp_dchi(1); -const vreal dchi3 = tmp_dchi(2); -const vreal dgamt111 = tmp_dgamt(0,0)(0); -const vreal dgamt112 = tmp_dgamt(0,1)(0); -const vreal dgamt113 = tmp_dgamt(0,2)(0); -const vreal dgamt122 = tmp_dgamt(1,1)(0); -const vreal dgamt123 = tmp_dgamt(1,2)(0); -const vreal dgamt133 = tmp_dgamt(2,2)(0); -const vreal dgamt211 = tmp_dgamt(0,0)(1); -const vreal dgamt212 = tmp_dgamt(0,1)(1); -const vreal dgamt213 = tmp_dgamt(0,2)(1); -const vreal dgamt222 = tmp_dgamt(1,1)(1); -const vreal dgamt223 = tmp_dgamt(1,2)(1); -const vreal dgamt233 = tmp_dgamt(2,2)(1); -const vreal dgamt311 = tmp_dgamt(0,0)(2); -const vreal dgamt312 = tmp_dgamt(0,1)(2); -const vreal dgamt313 = tmp_dgamt(0,2)(2); -const vreal dgamt322 = tmp_dgamt(1,1)(2); -const vreal dgamt323 = tmp_dgamt(1,2)(2); -const vreal dgamt333 = tmp_dgamt(2,2)(2); -const vreal dexKh1 = tmp_dexKh(0); -const vreal dexKh2 = tmp_dexKh(1); -const vreal dexKh3 = tmp_dexKh(2); -const vreal dexAt111 = tmp_dexAt(0,0)(0); -const vreal dexAt112 = tmp_dexAt(0,1)(0); -const vreal dexAt113 = tmp_dexAt(0,2)(0); -const vreal dexAt122 = tmp_dexAt(1,1)(0); -const vreal dexAt123 = tmp_dexAt(1,2)(0); -const vreal dexAt133 = tmp_dexAt(2,2)(0); -const vreal dexAt211 = tmp_dexAt(0,0)(1); -const vreal dexAt212 = tmp_dexAt(0,1)(1); -const vreal dexAt213 = tmp_dexAt(0,2)(1); -const vreal dexAt222 = tmp_dexAt(1,1)(1); -const vreal dexAt223 = tmp_dexAt(1,2)(1); -const vreal dexAt233 = tmp_dexAt(2,2)(1); -const vreal dexAt311 = tmp_dexAt(0,0)(2); -const vreal dexAt312 = tmp_dexAt(0,1)(2); -const vreal dexAt313 = tmp_dexAt(0,2)(2); -const vreal dexAt322 = tmp_dexAt(1,1)(2); -const vreal dexAt323 = tmp_dexAt(1,2)(2); -const vreal dexAt333 = tmp_dexAt(2,2)(2); -const vreal dtrGt11 = tmp_dtrGt(0)(0); -const vreal dtrGt12 = tmp_dtrGt(1)(0); -const vreal dtrGt13 = tmp_dtrGt(2)(0); -const vreal dtrGt21 = tmp_dtrGt(0)(1); -const vreal dtrGt22 = tmp_dtrGt(1)(1); -const vreal dtrGt23 = tmp_dtrGt(2)(1); -const vreal dtrGt31 = tmp_dtrGt(0)(2); -const vreal dtrGt32 = tmp_dtrGt(1)(2); -const vreal dtrGt33 = tmp_dtrGt(2)(2); -const vreal dTheta1 = tmp_dTheta(0); -const vreal dTheta2 = tmp_dTheta(1); -const vreal dTheta3 = tmp_dTheta(2); -const vreal ddchi11 = tmp_ddchi(0,0); -const vreal ddchi12 = tmp_ddchi(0,1); -const vreal ddchi13 = tmp_ddchi(0,2); -const vreal ddchi22 = tmp_ddchi(1,1); -const vreal ddchi23 = tmp_ddchi(1,2); -const vreal ddchi33 = tmp_ddchi(2,2); -const vreal ddgamt1111 = tmp_ddgamt(0,0)(0,0); -const vreal ddgamt1112 = tmp_ddgamt(0,1)(0,0); -const vreal ddgamt1113 = tmp_ddgamt(0,2)(0,0); -const vreal ddgamt1122 = tmp_ddgamt(1,1)(0,0); -const vreal ddgamt1123 = tmp_ddgamt(1,2)(0,0); -const vreal ddgamt1133 = tmp_ddgamt(2,2)(0,0); -const vreal ddgamt1211 = tmp_ddgamt(0,0)(0,1); -const vreal ddgamt1212 = tmp_ddgamt(0,1)(0,1); -const vreal ddgamt1213 = tmp_ddgamt(0,2)(0,1); -const vreal ddgamt1222 = tmp_ddgamt(1,1)(0,1); -const vreal ddgamt1223 = tmp_ddgamt(1,2)(0,1); -const vreal ddgamt1233 = tmp_ddgamt(2,2)(0,1); -const vreal ddgamt1311 = tmp_ddgamt(0,0)(0,2); -const vreal ddgamt1312 = tmp_ddgamt(0,1)(0,2); -const vreal ddgamt1313 = tmp_ddgamt(0,2)(0,2); -const vreal ddgamt1322 = tmp_ddgamt(1,1)(0,2); -const vreal ddgamt1323 = tmp_ddgamt(1,2)(0,2); -const vreal ddgamt1333 = tmp_ddgamt(2,2)(0,2); -const vreal ddgamt2211 = tmp_ddgamt(0,0)(1,1); -const vreal ddgamt2212 = tmp_ddgamt(0,1)(1,1); -const vreal ddgamt2213 = tmp_ddgamt(0,2)(1,1); -const vreal ddgamt2222 = tmp_ddgamt(1,1)(1,1); -const vreal ddgamt2223 = tmp_ddgamt(1,2)(1,1); -const vreal ddgamt2233 = tmp_ddgamt(2,2)(1,1); -const vreal ddgamt2311 = tmp_ddgamt(0,0)(1,2); -const vreal ddgamt2312 = tmp_ddgamt(0,1)(1,2); -const vreal ddgamt2313 = tmp_ddgamt(0,2)(1,2); -const vreal ddgamt2322 = tmp_ddgamt(1,1)(1,2); -const vreal ddgamt2323 = tmp_ddgamt(1,2)(1,2); -const vreal ddgamt2333 = tmp_ddgamt(2,2)(1,2); -const vreal ddgamt3311 = tmp_ddgamt(0,0)(2,2); -const vreal ddgamt3312 = tmp_ddgamt(0,1)(2,2); -const vreal ddgamt3313 = tmp_ddgamt(0,2)(2,2); -const vreal ddgamt3322 = tmp_ddgamt(1,1)(2,2); -const vreal ddgamt3323 = tmp_ddgamt(1,2)(2,2); -const vreal ddgamt3333 = tmp_ddgamt(2,2)(2,2); - -vreal invgamt11 -= --Power(gamt23,2) + gamt22*gamt33 -; - -vreal invgamt12 -= -gamt13*gamt23 - gamt12*gamt33 -; - -vreal invgamt13 -= --(gamt13*gamt22) + gamt12*gamt23 -; - -vreal invgamt22 -= --Power(gamt13,2) + gamt11*gamt33 -; - -vreal invgamt23 -= -gamt12*gamt13 - gamt11*gamt23 -; - -vreal invgamt33 -= --Power(gamt12,2) + gamt11*gamt22 -; - -vreal invgam11 -= -chi*invgamt11 -; - -vreal invgam12 -= -chi*invgamt12 -; - -vreal invgam13 -= -chi*invgamt13 -; - -vreal invgam22 -= -chi*invgamt22 -; - -vreal invgam23 -= -chi*invgamt23 -; - -vreal invgam33 -= -chi*invgamt33 -; - -vreal GtDDD111 -= -dgamt111/2. -; - -vreal GtDDD112 -= -dgamt211/2. -; - -vreal GtDDD113 -= -dgamt311/2. -; - -vreal GtDDD122 -= --0.5*dgamt122 + dgamt212 -; - -vreal GtDDD123 -= -(-dgamt123 + dgamt213 + dgamt312)/2. -; - -vreal GtDDD133 -= --0.5*dgamt133 + dgamt313 -; - -vreal GtDDD211 -= -dgamt112 - dgamt211/2. -; - -vreal GtDDD212 -= -dgamt122/2. -; - -vreal GtDDD213 -= -(dgamt123 - dgamt213 + dgamt312)/2. -; - -vreal GtDDD222 -= -dgamt222/2. -; - -vreal GtDDD223 -= -dgamt322/2. -; - -vreal GtDDD233 -= --0.5*dgamt233 + dgamt323 -; - -vreal GtDDD311 -= -dgamt113 - dgamt311/2. -; - -vreal GtDDD312 -= -(dgamt123 + dgamt213 - dgamt312)/2. -; - -vreal GtDDD313 -= -dgamt133/2. -; - -vreal GtDDD322 -= -dgamt223 - dgamt322/2. -; - -vreal GtDDD323 -= -dgamt233/2. -; - -vreal GtDDD333 -= -dgamt333/2. -; - -vreal GtDDU111 -= -GtDDD111*invgamt11 + GtDDD112*invgamt12 + GtDDD113*invgamt13 -; - -vreal GtDDU112 -= -GtDDD111*invgamt12 + GtDDD112*invgamt22 + GtDDD113*invgamt23 -; - -vreal GtDDU113 -= -GtDDD111*invgamt13 + GtDDD112*invgamt23 + GtDDD113*invgamt33 -; - -vreal GtDDU121 -= -GtDDD112*invgamt11 + GtDDD122*invgamt12 + GtDDD123*invgamt13 -; - -vreal GtDDU122 -= -GtDDD112*invgamt12 + GtDDD122*invgamt22 + GtDDD123*invgamt23 -; - -vreal GtDDU123 -= -GtDDD112*invgamt13 + GtDDD122*invgamt23 + GtDDD123*invgamt33 -; - -vreal GtDDU131 -= -GtDDD113*invgamt11 + GtDDD123*invgamt12 + GtDDD133*invgamt13 -; - -vreal GtDDU132 -= -GtDDD113*invgamt12 + GtDDD123*invgamt22 + GtDDD133*invgamt23 -; - -vreal GtDDU133 -= -GtDDD113*invgamt13 + GtDDD123*invgamt23 + GtDDD133*invgamt33 -; - -vreal GtDDU211 -= -GtDDD211*invgamt11 + GtDDD212*invgamt12 + GtDDD213*invgamt13 -; - -vreal GtDDU212 -= -GtDDD211*invgamt12 + GtDDD212*invgamt22 + GtDDD213*invgamt23 -; - -vreal GtDDU213 -= -GtDDD211*invgamt13 + GtDDD212*invgamt23 + GtDDD213*invgamt33 -; - -vreal GtDDU221 -= -GtDDD212*invgamt11 + GtDDD222*invgamt12 + GtDDD223*invgamt13 -; - -vreal GtDDU222 -= -GtDDD212*invgamt12 + GtDDD222*invgamt22 + GtDDD223*invgamt23 -; - -vreal GtDDU223 -= -GtDDD212*invgamt13 + GtDDD222*invgamt23 + GtDDD223*invgamt33 -; - -vreal GtDDU231 -= -GtDDD213*invgamt11 + GtDDD223*invgamt12 + GtDDD233*invgamt13 -; - -vreal GtDDU232 -= -GtDDD213*invgamt12 + GtDDD223*invgamt22 + GtDDD233*invgamt23 -; - -vreal GtDDU233 -= -GtDDD213*invgamt13 + GtDDD223*invgamt23 + GtDDD233*invgamt33 -; - -vreal GtDDU311 -= -GtDDD311*invgamt11 + GtDDD312*invgamt12 + GtDDD313*invgamt13 -; - -vreal GtDDU312 -= -GtDDD311*invgamt12 + GtDDD312*invgamt22 + GtDDD313*invgamt23 -; - -vreal GtDDU313 -= -GtDDD311*invgamt13 + GtDDD312*invgamt23 + GtDDD313*invgamt33 -; - -vreal GtDDU321 -= -GtDDD312*invgamt11 + GtDDD322*invgamt12 + GtDDD323*invgamt13 -; - -vreal GtDDU322 -= -GtDDD312*invgamt12 + GtDDD322*invgamt22 + GtDDD323*invgamt23 -; - -vreal GtDDU323 -= -GtDDD312*invgamt13 + GtDDD322*invgamt23 + GtDDD323*invgamt33 -; - -vreal GtDDU331 -= -GtDDD313*invgamt11 + GtDDD323*invgamt12 + GtDDD333*invgamt13 -; - -vreal GtDDU332 -= -GtDDD313*invgamt12 + GtDDD323*invgamt22 + GtDDD333*invgamt23 -; - -vreal GtDDU333 -= -GtDDD313*invgamt13 + GtDDD323*invgamt23 + GtDDD333*invgamt33 -; - -vreal Gt111 -= -GtDDD111*invgamt11 + GtDDD211*invgamt12 + GtDDD311*invgamt13 -; - -vreal Gt112 -= -GtDDD112*invgamt11 + GtDDD212*invgamt12 + GtDDD312*invgamt13 -; - -vreal Gt113 -= -GtDDD113*invgamt11 + GtDDD213*invgamt12 + GtDDD313*invgamt13 -; - -vreal Gt122 -= -GtDDD122*invgamt11 + GtDDD222*invgamt12 + GtDDD322*invgamt13 -; - -vreal Gt123 -= -GtDDD123*invgamt11 + GtDDD223*invgamt12 + GtDDD323*invgamt13 -; - -vreal Gt133 -= -GtDDD133*invgamt11 + GtDDD233*invgamt12 + GtDDD333*invgamt13 -; - -vreal Gt211 -= -GtDDD111*invgamt12 + GtDDD211*invgamt22 + GtDDD311*invgamt23 -; - -vreal Gt212 -= -GtDDD112*invgamt12 + GtDDD212*invgamt22 + GtDDD312*invgamt23 -; - -vreal Gt213 -= -GtDDD113*invgamt12 + GtDDD213*invgamt22 + GtDDD313*invgamt23 -; - -vreal Gt222 -= -GtDDD122*invgamt12 + GtDDD222*invgamt22 + GtDDD322*invgamt23 -; - -vreal Gt223 -= -GtDDD123*invgamt12 + GtDDD223*invgamt22 + GtDDD323*invgamt23 -; - -vreal Gt233 -= -GtDDD133*invgamt12 + GtDDD233*invgamt22 + GtDDD333*invgamt23 -; - -vreal Gt311 -= -GtDDD111*invgamt13 + GtDDD211*invgamt23 + GtDDD311*invgamt33 -; - -vreal Gt312 -= -GtDDD112*invgamt13 + GtDDD212*invgamt23 + GtDDD312*invgamt33 -; - -vreal Gt313 -= -GtDDD113*invgamt13 + GtDDD213*invgamt23 + GtDDD313*invgamt33 -; - -vreal Gt322 -= -GtDDD122*invgamt13 + GtDDD222*invgamt23 + GtDDD322*invgamt33 -; - -vreal Gt323 -= -GtDDD123*invgamt13 + GtDDD223*invgamt23 + GtDDD323*invgamt33 -; - -vreal Gt333 -= -GtDDD133*invgamt13 + GtDDD233*invgamt23 + GtDDD333*invgamt33 -; - -vreal trGtd1 -= -Gt111*invgamt11 + 2*Gt112*invgamt12 + 2*Gt113*invgamt13 + Gt122*invgamt22 + - 2*Gt123*invgamt23 + Gt133*invgamt33 -; - -vreal trGtd2 -= -Gt211*invgamt11 + 2*Gt212*invgamt12 + 2*Gt213*invgamt13 + Gt222*invgamt22 + - 2*Gt223*invgamt23 + Gt233*invgamt33 -; - -vreal trGtd3 -= -Gt311*invgamt11 + 2*Gt312*invgamt12 + 2*Gt313*invgamt13 + Gt322*invgamt22 + - 2*Gt323*invgamt23 + Gt333*invgamt33 -; - -vreal exAtUU11 -= -exAt11*Power(invgamt11,2) + 2*exAt12*invgamt11*invgamt12 + - exAt22*Power(invgamt12,2) + 2*exAt13*invgamt11*invgamt13 + - 2*exAt23*invgamt12*invgamt13 + exAt33*Power(invgamt13,2) -; - -vreal exAtUU12 -= -exAt11*invgamt11*invgamt12 + exAt13*invgamt12*invgamt13 + - exAt22*invgamt12*invgamt22 + exAt23*invgamt13*invgamt22 + - exAt12*(Power(invgamt12,2) + invgamt11*invgamt22) + - exAt13*invgamt11*invgamt23 + exAt23*invgamt12*invgamt23 + - exAt33*invgamt13*invgamt23 -; - -vreal exAtUU13 -= -exAt11*invgamt11*invgamt13 + exAt12*invgamt12*invgamt13 + - exAt13*Power(invgamt13,2) + exAt12*invgamt11*invgamt23 + - exAt22*invgamt12*invgamt23 + exAt23*invgamt13*invgamt23 + - exAt13*invgamt11*invgamt33 + exAt23*invgamt12*invgamt33 + - exAt33*invgamt13*invgamt33 -; - -vreal exAtUU22 -= -exAt11*Power(invgamt12,2) + 2*exAt12*invgamt12*invgamt22 + - exAt22*Power(invgamt22,2) + 2*exAt13*invgamt12*invgamt23 + - 2*exAt23*invgamt22*invgamt23 + exAt33*Power(invgamt23,2) -; - -vreal exAtUU23 -= -exAt11*invgamt12*invgamt13 + exAt12*invgamt13*invgamt22 + - exAt12*invgamt12*invgamt23 + exAt13*invgamt13*invgamt23 + - exAt22*invgamt22*invgamt23 + exAt23*Power(invgamt23,2) + - exAt13*invgamt12*invgamt33 + exAt23*invgamt22*invgamt33 + - exAt33*invgamt23*invgamt33 -; - -vreal exAtUU33 -= -exAt11*Power(invgamt13,2) + 2*exAt12*invgamt13*invgamt23 + - exAt22*Power(invgamt23,2) + 2*exAt13*invgamt13*invgamt33 + - 2*exAt23*invgamt23*invgamt33 + exAt33*Power(invgamt33,2) -; - -vreal tDtDchi11 -= -ddchi11 - dchi1*Gt111 - dchi2*Gt211 - dchi3*Gt311 -; - -vreal tDtDchi12 -= -ddchi12 - dchi1*Gt112 - dchi2*Gt212 - dchi3*Gt312 -; - -vreal tDtDchi13 -= -ddchi13 - dchi1*Gt113 - dchi2*Gt213 - dchi3*Gt313 -; - -vreal tDtDchi22 -= -ddchi22 - dchi1*Gt122 - dchi2*Gt222 - dchi3*Gt322 -; - -vreal tDtDchi23 -= -ddchi23 - dchi1*Gt123 - dchi2*Gt223 - dchi3*Gt323 -; - -vreal tDtDchi33 -= -ddchi33 - dchi1*Gt133 - dchi2*Gt233 - dchi3*Gt333 -; - -vreal Rtchi11 -= -(-(Power(dchi1,2)*(1 + 3*gamt11*invgamt11)) - - 6*dchi1*gamt11*(dchi2*invgamt12 + dchi3*invgamt13) - - 3*Power(dchi2,2)*gamt11*invgamt22 - 6*dchi2*dchi3*gamt11*invgamt23 - - 3*Power(dchi3,2)*gamt11*invgamt33 + 2*chi*tDtDchi11 + - 2*chi*gamt11*invgamt11*tDtDchi11 + 4*chi*gamt11*invgamt12*tDtDchi12 + - 4*chi*gamt11*invgamt13*tDtDchi13 + 2*chi*gamt11*invgamt22*tDtDchi22 + - 4*chi*gamt11*invgamt23*tDtDchi23 + 2*chi*gamt11*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi12 -= -(-3*Power(dchi1,2)*gamt12*invgamt11 - - dchi1*(dchi2 + 6*dchi2*gamt12*invgamt12 + 6*dchi3*gamt12*invgamt13) - - 3*Power(dchi2,2)*gamt12*invgamt22 - 6*dchi2*dchi3*gamt12*invgamt23 - - 3*Power(dchi3,2)*gamt12*invgamt33 + 2*chi*gamt12*invgamt11*tDtDchi11 + - 2*chi*tDtDchi12 + 4*chi*gamt12*invgamt12*tDtDchi12 + - 4*chi*gamt12*invgamt13*tDtDchi13 + 2*chi*gamt12*invgamt22*tDtDchi22 + - 4*chi*gamt12*invgamt23*tDtDchi23 + 2*chi*gamt12*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi13 -= -(-3*Power(dchi1,2)*gamt13*invgamt11 - - dchi1*(dchi3 + 6*dchi2*gamt13*invgamt12 + 6*dchi3*gamt13*invgamt13) - - 3*Power(dchi2,2)*gamt13*invgamt22 - 6*dchi2*dchi3*gamt13*invgamt23 - - 3*Power(dchi3,2)*gamt13*invgamt33 + 2*chi*gamt13*invgamt11*tDtDchi11 + - 4*chi*gamt13*invgamt12*tDtDchi12 + 2*chi*tDtDchi13 + - 4*chi*gamt13*invgamt13*tDtDchi13 + 2*chi*gamt13*invgamt22*tDtDchi22 + - 4*chi*gamt13*invgamt23*tDtDchi23 + 2*chi*gamt13*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi22 -= -(-3*Power(dchi1,2)*gamt22*invgamt11 - 6*dchi1*dchi3*gamt22*invgamt13 - - Power(dchi2,2)*(1 + 3*gamt22*invgamt22) - - 6*dchi2*gamt22*(dchi1*invgamt12 + dchi3*invgamt23) - - 3*Power(dchi3,2)*gamt22*invgamt33 + 2*chi*gamt22*invgamt11*tDtDchi11 + - 4*chi*gamt22*invgamt12*tDtDchi12 + 4*chi*gamt22*invgamt13*tDtDchi13 + - 2*chi*tDtDchi22 + 2*chi*gamt22*invgamt22*tDtDchi22 + - 4*chi*gamt22*invgamt23*tDtDchi23 + 2*chi*gamt22*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi23 -= -(-3*Power(dchi1,2)*gamt23*invgamt11 - 6*dchi1*dchi3*gamt23*invgamt13 - - 3*Power(dchi2,2)*gamt23*invgamt22 - - dchi2*(dchi3 + 6*dchi1*gamt23*invgamt12 + 6*dchi3*gamt23*invgamt23) - - 3*Power(dchi3,2)*gamt23*invgamt33 + 2*chi*gamt23*invgamt11*tDtDchi11 + - 4*chi*gamt23*invgamt12*tDtDchi12 + 4*chi*gamt23*invgamt13*tDtDchi13 + - 2*chi*gamt23*invgamt22*tDtDchi22 + 2*chi*tDtDchi23 + - 4*chi*gamt23*invgamt23*tDtDchi23 + 2*chi*gamt23*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi33 -= -(-3*Power(dchi1,2)*gamt33*invgamt11 - 6*dchi1*dchi2*gamt33*invgamt12 - - 3*Power(dchi2,2)*gamt33*invgamt22 - - 6*dchi3*gamt33*(dchi1*invgamt13 + dchi2*invgamt23) - - Power(dchi3,2)*(1 + 3*gamt33*invgamt33) + - 2*chi*gamt33*invgamt11*tDtDchi11 + 4*chi*gamt33*invgamt12*tDtDchi12 + - 4*chi*gamt33*invgamt13*tDtDchi13 + 2*chi*gamt33*invgamt22*tDtDchi22 + - 4*chi*gamt33*invgamt23*tDtDchi23 + 2*chi*tDtDchi33 + - 2*chi*gamt33*invgamt33*tDtDchi33)/(4.*Power(chi,2)) -; - -vreal Rt11 -= -dtrGt11*gamt11 + dtrGt12*gamt12 + dtrGt13*gamt13 + 3*Gt111*GtDDU111 + - 3*Gt112*GtDDU112 + 3*Gt113*GtDDU113 + 2*Gt211*GtDDU121 + - 2*Gt212*GtDDU122 + 2*Gt213*GtDDU123 + 2*Gt311*GtDDU131 + - 2*Gt312*GtDDU132 + 2*Gt313*GtDDU133 + Gt211*GtDDU211 + Gt212*GtDDU212 + - Gt213*GtDDU213 + Gt311*GtDDU311 + Gt312*GtDDU312 + Gt313*GtDDU313 - - (ddgamt1111*invgamt11)/2. - ddgamt1211*invgamt12 - ddgamt1311*invgamt13 - - (ddgamt2211*invgamt22)/2. - ddgamt2311*invgamt23 - - (ddgamt3311*invgamt33)/2. + GtDDD111*trGtd1 + GtDDD112*trGtd2 + - GtDDD113*trGtd3 -; - -vreal Rt12 -= -(dtrGt21*gamt11 + dtrGt11*gamt12 + dtrGt22*gamt12 + dtrGt23*gamt13 + - dtrGt12*gamt22 + dtrGt13*gamt23 + 2*Gt112*GtDDU111 + 2*Gt122*GtDDU112 + - 2*Gt123*GtDDU113 + 2*Gt111*GtDDU121 + 2*Gt212*GtDDU121 + - 2*Gt112*GtDDU122 + 2*Gt222*GtDDU122 + 2*Gt113*GtDDU123 + - 2*Gt223*GtDDU123 + 2*Gt312*GtDDU131 + 2*Gt322*GtDDU132 + - 2*Gt323*GtDDU133 + 2*Gt111*GtDDU211 + 2*Gt112*GtDDU212 + - 2*Gt113*GtDDU213 + 4*Gt211*GtDDU221 + 4*Gt212*GtDDU222 + - 4*Gt213*GtDDU223 + 2*Gt311*GtDDU231 + 2*Gt312*GtDDU232 + - 2*Gt313*GtDDU233 + 2*Gt311*GtDDU321 + 2*Gt312*GtDDU322 + - 2*Gt313*GtDDU323 - ddgamt1112*invgamt11 - 2*ddgamt1212*invgamt12 - - 2*ddgamt1312*invgamt13 - ddgamt2212*invgamt22 - - 2*ddgamt2312*invgamt23 - ddgamt3312*invgamt33 + GtDDD112*trGtd1 + - GtDDD211*trGtd1 + GtDDD122*trGtd2 + GtDDD212*trGtd2 + GtDDD123*trGtd3 + - GtDDD213*trGtd3)/2. -; - -vreal Rt13 -= -(dtrGt31*gamt11 + dtrGt32*gamt12 + dtrGt11*gamt13 + dtrGt33*gamt13 + - dtrGt12*gamt23 + dtrGt13*gamt33 + 2*Gt113*GtDDU111 + 2*Gt123*GtDDU112 + - 2*Gt133*GtDDU113 + 2*Gt213*GtDDU121 + 2*Gt223*GtDDU122 + - 2*Gt233*GtDDU123 + 2*Gt111*GtDDU131 + 2*Gt313*GtDDU131 + - 2*Gt112*GtDDU132 + 2*Gt323*GtDDU132 + 2*Gt113*GtDDU133 + - 2*Gt333*GtDDU133 + 2*Gt211*GtDDU231 + 2*Gt212*GtDDU232 + - 2*Gt213*GtDDU233 + 2*Gt111*GtDDU311 + 2*Gt112*GtDDU312 + - 2*Gt113*GtDDU313 + 2*Gt211*GtDDU321 + 2*Gt212*GtDDU322 + - 2*Gt213*GtDDU323 + 4*Gt311*GtDDU331 + 4*Gt312*GtDDU332 + - 4*Gt313*GtDDU333 - ddgamt1113*invgamt11 - 2*ddgamt1213*invgamt12 - - 2*ddgamt1313*invgamt13 - ddgamt2213*invgamt22 - - 2*ddgamt2313*invgamt23 - ddgamt3313*invgamt33 + GtDDD113*trGtd1 + - GtDDD311*trGtd1 + GtDDD123*trGtd2 + GtDDD312*trGtd2 + GtDDD133*trGtd3 + - GtDDD313*trGtd3)/2. -; - -vreal Rt22 -= -dtrGt21*gamt12 + dtrGt22*gamt22 + dtrGt23*gamt23 + Gt112*GtDDU121 + - Gt122*GtDDU122 + Gt123*GtDDU123 + 2*Gt112*GtDDU211 + 2*Gt122*GtDDU212 + - 2*Gt123*GtDDU213 + 3*Gt212*GtDDU221 + 3*Gt222*GtDDU222 + - 3*Gt223*GtDDU223 + 2*Gt312*GtDDU231 + 2*Gt322*GtDDU232 + - 2*Gt323*GtDDU233 + Gt312*GtDDU321 + Gt322*GtDDU322 + Gt323*GtDDU323 - - (ddgamt1122*invgamt11)/2. - ddgamt1222*invgamt12 - ddgamt1322*invgamt13 - - (ddgamt2222*invgamt22)/2. - ddgamt2322*invgamt23 - - (ddgamt3322*invgamt33)/2. + GtDDD212*trGtd1 + GtDDD222*trGtd2 + - GtDDD223*trGtd3 -; - -vreal Rt23 -= -(dtrGt31*gamt12 + dtrGt21*gamt13 + dtrGt32*gamt22 + dtrGt22*gamt23 + - dtrGt33*gamt23 + dtrGt23*gamt33 + 2*Gt112*GtDDU131 + 2*Gt122*GtDDU132 + - 2*Gt123*GtDDU133 + 2*Gt113*GtDDU211 + 2*Gt123*GtDDU212 + - 2*Gt133*GtDDU213 + 2*Gt213*GtDDU221 + 2*Gt223*GtDDU222 + - 2*Gt233*GtDDU223 + 2*Gt212*GtDDU231 + 2*Gt313*GtDDU231 + - 2*Gt222*GtDDU232 + 2*Gt323*GtDDU232 + 2*Gt223*GtDDU233 + - 2*Gt333*GtDDU233 + 2*Gt112*GtDDU311 + 2*Gt122*GtDDU312 + - 2*Gt123*GtDDU313 + 2*Gt212*GtDDU321 + 2*Gt222*GtDDU322 + - 2*Gt223*GtDDU323 + 4*Gt312*GtDDU331 + 4*Gt322*GtDDU332 + - 4*Gt323*GtDDU333 - ddgamt1123*invgamt11 - 2*ddgamt1223*invgamt12 - - 2*ddgamt1323*invgamt13 - ddgamt2223*invgamt22 - - 2*ddgamt2323*invgamt23 - ddgamt3323*invgamt33 + GtDDD213*trGtd1 + - GtDDD312*trGtd1 + GtDDD223*trGtd2 + GtDDD322*trGtd2 + GtDDD233*trGtd3 + - GtDDD323*trGtd3)/2. -; - -vreal Rt33 -= -dtrGt31*gamt13 + dtrGt32*gamt23 + dtrGt33*gamt33 + Gt113*GtDDU131 + - Gt123*GtDDU132 + Gt133*GtDDU133 + Gt213*GtDDU231 + Gt223*GtDDU232 + - Gt233*GtDDU233 + 2*Gt113*GtDDU311 + 2*Gt123*GtDDU312 + 2*Gt133*GtDDU313 + - 2*Gt213*GtDDU321 + 2*Gt223*GtDDU322 + 2*Gt233*GtDDU323 + - 3*Gt313*GtDDU331 + 3*Gt323*GtDDU332 + 3*Gt333*GtDDU333 - - (ddgamt1133*invgamt11)/2. - ddgamt1233*invgamt12 - ddgamt1333*invgamt13 - - (ddgamt2233*invgamt22)/2. - ddgamt2333*invgamt23 - - (ddgamt3333*invgamt33)/2. + GtDDD313*trGtd1 + GtDDD323*trGtd2 + - GtDDD333*trGtd3 -; - -vreal R11 -= -Rt11 + Rtchi11 -; - -vreal R12 -= -Rt12 + Rtchi12 -; - -vreal R13 -= -Rt13 + Rtchi13 -; - -vreal R22 -= -Rt22 + Rtchi22 -; - -vreal R23 -= -Rt23 + Rtchi23 -; - -vreal R33 -= -Rt33 + Rtchi33 -; - -vreal trR -= -invgam11*R11 + 2*invgam12*R12 + 2*invgam13*R13 + invgam22*R22 + - 2*invgam23*R23 + invgam33*R33 -; - -vreal rho -= -(Power(beta1,2)*eT11 + Power(beta2,2)*eT22 + 2*beta2*beta3*eT23 + - Power(beta3,2)*eT33 + 2*beta1*(beta2*eT12 + beta3*eT13 - eTt1) - - 2*beta2*eTt2 - 2*beta3*eTt3 + eTtt)/Power(alpha,2) -; - -vreal Sm1 -= -(beta1*eT11 + beta2*eT12 + beta3*eT13 - eTt1)/alpha -; - -vreal Sm2 -= -(beta1*eT12 + beta2*eT22 + beta3*eT23 - eTt2)/alpha -; - -vreal Sm3 -= -(beta1*eT13 + beta2*eT23 + beta3*eT33 - eTt3)/alpha -; - -vreal trdexAtUU1 -= --2*dgamt111*exAtUU11*invgamt11 - dgamt211*exAtUU12*invgamt11 - - 2*dgamt113*exAtUU13*invgamt11 - dgamt311*exAtUU13*invgamt11 - - dgamt212*exAtUU22*invgamt11 - dgamt213*exAtUU23*invgamt11 - - dgamt312*exAtUU23*invgamt11 - dgamt313*exAtUU33*invgamt11 + - dexAt111*Power(invgamt11,2) - dgamt211*exAtUU11*invgamt12 - - 2*dgamt122*exAtUU12*invgamt12 - 2*dgamt212*exAtUU12*invgamt12 - - 2*dgamt123*exAtUU13*invgamt12 - dgamt213*exAtUU13*invgamt12 - - dgamt312*exAtUU13*invgamt12 - dgamt222*exAtUU22*invgamt12 - - dgamt223*exAtUU23*invgamt12 - dgamt322*exAtUU23*invgamt12 - - dgamt323*exAtUU33*invgamt12 + 2*dexAt112*invgamt11*invgamt12 + - dexAt211*invgamt11*invgamt12 + dexAt122*Power(invgamt12,2) + - dexAt212*Power(invgamt12,2) - - 2*dgamt112*(exAtUU12*invgamt11 + exAtUU11*invgamt12) - - 2*dgamt113*exAtUU11*invgamt13 - dgamt311*exAtUU11*invgamt13 - - 2*dgamt123*exAtUU12*invgamt13 - dgamt213*exAtUU12*invgamt13 - - dgamt312*exAtUU12*invgamt13 - 2*dgamt133*exAtUU13*invgamt13 - - 2*dgamt313*exAtUU13*invgamt13 - dgamt223*exAtUU22*invgamt13 - - dgamt233*exAtUU23*invgamt13 - dgamt323*exAtUU23*invgamt13 - - dgamt333*exAtUU33*invgamt13 + 2*dexAt113*invgamt11*invgamt13 + - dexAt311*invgamt11*invgamt13 + 2*dexAt123*invgamt12*invgamt13 + - dexAt213*invgamt12*invgamt13 + dexAt312*invgamt12*invgamt13 + - dexAt133*Power(invgamt13,2) + dexAt313*Power(invgamt13,2) - - dgamt212*exAtUU11*invgamt22 - dgamt222*exAtUU12*invgamt22 - - dgamt223*exAtUU13*invgamt22 + dexAt212*invgamt11*invgamt22 + - dexAt222*invgamt12*invgamt22 + dexAt223*invgamt13*invgamt22 - - dgamt213*exAtUU11*invgamt23 - dgamt312*exAtUU11*invgamt23 - - dgamt223*exAtUU12*invgamt23 - dgamt322*exAtUU12*invgamt23 - - dgamt233*exAtUU13*invgamt23 - dgamt323*exAtUU13*invgamt23 + - dexAt213*invgamt11*invgamt23 + dexAt312*invgamt11*invgamt23 + - dexAt223*invgamt12*invgamt23 + dexAt322*invgamt12*invgamt23 + - dexAt233*invgamt13*invgamt23 + dexAt323*invgamt13*invgamt23 - - dgamt313*exAtUU11*invgamt33 - dgamt323*exAtUU12*invgamt33 - - dgamt333*exAtUU13*invgamt33 + dexAt313*invgamt11*invgamt33 + - dexAt323*invgamt12*invgamt33 + dexAt333*invgamt13*invgamt33 -; - -vreal trdexAtUU2 -= --(dgamt113*exAtUU23*invgamt11) - 2*dgamt211*exAtUU12*invgamt12 - - dgamt113*exAtUU13*invgamt12 - dgamt311*exAtUU13*invgamt12 - - dgamt122*exAtUU22*invgamt12 - 2*dgamt212*exAtUU22*invgamt12 - - dgamt123*exAtUU23*invgamt12 - 2*dgamt213*exAtUU23*invgamt12 - - dgamt312*exAtUU23*invgamt12 - dgamt313*exAtUU33*invgamt12 + - dexAt111*invgamt11*invgamt12 + dexAt112*Power(invgamt12,2) + - dexAt211*Power(invgamt12,2) - - dgamt111*(exAtUU12*invgamt11 + exAtUU11*invgamt12) - - dgamt113*exAtUU12*invgamt13 - dgamt311*exAtUU12*invgamt13 - - dgamt123*exAtUU22*invgamt13 - dgamt312*exAtUU22*invgamt13 - - dgamt133*exAtUU23*invgamt13 - dgamt313*exAtUU23*invgamt13 + - dexAt113*invgamt12*invgamt13 + dexAt311*invgamt12*invgamt13 - - dgamt122*exAtUU12*invgamt22 - 2*dgamt212*exAtUU12*invgamt22 - - dgamt123*exAtUU13*invgamt22 - dgamt312*exAtUU13*invgamt22 - - 2*dgamt222*exAtUU22*invgamt22 - 2*dgamt223*exAtUU23*invgamt22 - - dgamt322*exAtUU23*invgamt22 - dgamt323*exAtUU33*invgamt22 + - dexAt112*invgamt11*invgamt22 + dexAt122*invgamt12*invgamt22 + - 2*dexAt212*invgamt12*invgamt22 + dexAt123*invgamt13*invgamt22 + - dexAt312*invgamt13*invgamt22 + dexAt222*Power(invgamt22,2) - - dgamt112*(exAtUU22*invgamt11 + 2*exAtUU12*invgamt12 + - exAtUU11*invgamt22) - dgamt113*exAtUU11*invgamt23 - - dgamt123*exAtUU12*invgamt23 - 2*dgamt213*exAtUU12*invgamt23 - - dgamt312*exAtUU12*invgamt23 - dgamt133*exAtUU13*invgamt23 - - dgamt313*exAtUU13*invgamt23 - 2*dgamt223*exAtUU22*invgamt23 - - dgamt322*exAtUU22*invgamt23 - 2*dgamt233*exAtUU23*invgamt23 - - 2*dgamt323*exAtUU23*invgamt23 - dgamt333*exAtUU33*invgamt23 + - dexAt113*invgamt11*invgamt23 + dexAt123*invgamt12*invgamt23 + - 2*dexAt213*invgamt12*invgamt23 + dexAt312*invgamt12*invgamt23 + - dexAt133*invgamt13*invgamt23 + dexAt313*invgamt13*invgamt23 + - 2*dexAt223*invgamt22*invgamt23 + dexAt322*invgamt22*invgamt23 + - dexAt233*Power(invgamt23,2) + dexAt323*Power(invgamt23,2) - - dgamt313*exAtUU12*invgamt33 - dgamt323*exAtUU22*invgamt33 - - dgamt333*exAtUU23*invgamt33 + dexAt313*invgamt12*invgamt33 + - dexAt323*invgamt22*invgamt33 + dexAt333*invgamt23*invgamt33 -; - -vreal trdexAtUU3 -= --(dgamt113*exAtUU33*invgamt11) - dgamt211*exAtUU13*invgamt12 - - dgamt122*exAtUU23*invgamt12 - dgamt212*exAtUU23*invgamt12 - - dgamt123*exAtUU33*invgamt12 - dgamt213*exAtUU33*invgamt12 - - dgamt211*exAtUU12*invgamt13 - 2*dgamt113*exAtUU13*invgamt13 - - 2*dgamt311*exAtUU13*invgamt13 - dgamt212*exAtUU22*invgamt13 - - dgamt123*exAtUU23*invgamt13 - dgamt213*exAtUU23*invgamt13 - - 2*dgamt312*exAtUU23*invgamt13 - dgamt133*exAtUU33*invgamt13 - - 2*dgamt313*exAtUU33*invgamt13 + dexAt111*invgamt11*invgamt13 + - dexAt112*invgamt12*invgamt13 + dexAt211*invgamt12*invgamt13 + - dexAt113*Power(invgamt13,2) + dexAt311*Power(invgamt13,2) - - dgamt111*(exAtUU13*invgamt11 + exAtUU11*invgamt13) - - dgamt212*exAtUU13*invgamt22 - dgamt222*exAtUU23*invgamt22 - - dgamt223*exAtUU33*invgamt22 + dexAt212*invgamt13*invgamt22 - - dgamt122*exAtUU12*invgamt23 - dgamt212*exAtUU12*invgamt23 - - dgamt123*exAtUU13*invgamt23 - dgamt213*exAtUU13*invgamt23 - - 2*dgamt312*exAtUU13*invgamt23 - dgamt222*exAtUU22*invgamt23 - - 2*dgamt223*exAtUU23*invgamt23 - 2*dgamt322*exAtUU23*invgamt23 - - dgamt233*exAtUU33*invgamt23 - 2*dgamt323*exAtUU33*invgamt23 + - dexAt112*invgamt11*invgamt23 + dexAt122*invgamt12*invgamt23 + - dexAt212*invgamt12*invgamt23 + dexAt123*invgamt13*invgamt23 + - dexAt213*invgamt13*invgamt23 + 2*dexAt312*invgamt13*invgamt23 + - dexAt222*invgamt22*invgamt23 + dexAt223*Power(invgamt23,2) + - dexAt322*Power(invgamt23,2) - - dgamt112*(exAtUU23*invgamt11 + exAtUU13*invgamt12 + exAtUU12*invgamt13 + - exAtUU11*invgamt23) - dgamt113*exAtUU11*invgamt33 - - dgamt123*exAtUU12*invgamt33 - dgamt213*exAtUU12*invgamt33 - - dgamt133*exAtUU13*invgamt33 - 2*dgamt313*exAtUU13*invgamt33 - - dgamt223*exAtUU22*invgamt33 - dgamt233*exAtUU23*invgamt33 - - 2*dgamt323*exAtUU23*invgamt33 - 2*dgamt333*exAtUU33*invgamt33 + - dexAt113*invgamt11*invgamt33 + dexAt123*invgamt12*invgamt33 + - dexAt213*invgamt12*invgamt33 + dexAt133*invgamt13*invgamt33 + - 2*dexAt313*invgamt13*invgamt33 + dexAt223*invgamt22*invgamt33 + - dexAt233*invgamt23*invgamt33 + 2*dexAt323*invgamt23*invgamt33 + - dexAt333*Power(invgamt33,2) -; - - -local_ZtC1.store(mask, index2, -(trGt1 - trGtd1)/2. -); - -local_ZtC2.store(mask, index2, -(trGt2 - trGtd2)/2. -); - -local_ZtC3.store(mask, index2, -(trGt3 - trGtd3)/2. -); - -local_HC.store(mask, index2, -exAt11*exAtUU11 + 2*exAt12*exAtUU12 + 2*exAt13*exAtUU13 + exAt22*exAtUU22 + - 2*exAt23*exAtUU23 + exAt33*exAtUU33 - (2*Power(exKh,2))/3. - 16*cpi*rho - - (8*exKh*Theta)/3. - (8*Power(Theta,2))/3. + trR -); - -local_MtC1.store(mask, index2, -((-2*(dchi1*exAtUU11 + dchi2*exAtUU12 + dchi3*exAtUU13))/chi + - 3*exAtUU11*Gt111 + 6*exAtUU12*Gt112 + 6*exAtUU13*Gt113 + - 3*exAtUU22*Gt122 + 6*exAtUU23*Gt123 + 3*exAtUU33*Gt133 - - 2*dexKh1*invgamt11 - 4*dTheta1*invgamt11 - 2*dexKh2*invgamt12 - - 4*dTheta2*invgamt12 - 2*dexKh3*invgamt13 - 4*dTheta3*invgamt13 - - 24*cpi*invgamt11*Sm1 - 24*cpi*invgamt12*Sm2 - 24*cpi*invgamt13*Sm3 + - 3*trdexAtUU1)/3. -); - -local_MtC2.store(mask, index2, -((-2*(dchi1*exAtUU12 + dchi2*exAtUU22 + dchi3*exAtUU23))/chi + - 3*exAtUU11*Gt211 + 6*exAtUU12*Gt212 + 6*exAtUU13*Gt213 + - 3*exAtUU22*Gt222 + 6*exAtUU23*Gt223 + 3*exAtUU33*Gt233 - - 2*dexKh1*invgamt12 - 4*dTheta1*invgamt12 - 2*dexKh2*invgamt22 - - 4*dTheta2*invgamt22 - 2*dexKh3*invgamt23 - 4*dTheta3*invgamt23 - - 24*cpi*invgamt12*Sm1 - 24*cpi*invgamt22*Sm2 - 24*cpi*invgamt23*Sm3 + - 3*trdexAtUU2)/3. -); - -local_MtC3.store(mask, index2, -((-2*(dchi1*exAtUU13 + dchi2*exAtUU23 + dchi3*exAtUU33))/chi + - 3*exAtUU11*Gt311 + 6*exAtUU12*Gt312 + 6*exAtUU13*Gt313 + - 3*exAtUU22*Gt322 + 6*exAtUU23*Gt323 + 3*exAtUU33*Gt333 - - 2*dexKh1*invgamt13 - 4*dTheta1*invgamt13 - 2*dexKh2*invgamt23 - - 4*dTheta2*invgamt23 - 2*dexKh3*invgamt33 - 4*dTheta3*invgamt33 - - 24*cpi*invgamt13*Sm1 - 24*cpi*invgamt23*Sm2 - 24*cpi*invgamt33*Sm3 + - 3*trdexAtUU3)/3. -); - - - }); -}); - -#endif // #ifndef Z4COO_SET_CONSTRAINT_HXX - -/* Z4coo_set_constraint.hxx */ diff --git a/Z4coo/wolfram/Z4coo_set_rhs.hxx b/Z4coo/wolfram/Z4coo_set_rhs.hxx deleted file mode 100644 index 093ace92..00000000 --- a/Z4coo/wolfram/Z4coo_set_rhs.hxx +++ /dev/null @@ -1,1523 +0,0 @@ -/* Z4coo_set_rhs.hxx */ -/* Produced with Mathematica */ - -#ifndef Z4COO_SET_RHS_HXX -#define Z4COO_SET_RHS_HXX - -const GF3D2 &local_dtchi = gf_dtchi; -const GF3D2 &local_dtgamt11 = gf_dtgamt(0,0); -const GF3D2 &local_dtgamt12 = gf_dtgamt(0,1); -const GF3D2 &local_dtgamt13 = gf_dtgamt(0,2); -const GF3D2 &local_dtgamt22 = gf_dtgamt(1,1); -const GF3D2 &local_dtgamt23 = gf_dtgamt(1,2); -const GF3D2 &local_dtgamt33 = gf_dtgamt(2,2); -const GF3D2 &local_dtexKh = gf_dtexKh; -const GF3D2 &local_dtexAt11 = gf_dtexAt(0,0); -const GF3D2 &local_dtexAt12 = gf_dtexAt(0,1); -const GF3D2 &local_dtexAt13 = gf_dtexAt(0,2); -const GF3D2 &local_dtexAt22 = gf_dtexAt(1,1); -const GF3D2 &local_dtexAt23 = gf_dtexAt(1,2); -const GF3D2 &local_dtexAt33 = gf_dtexAt(2,2); -const GF3D2 &local_dttrGt1 = gf_dttrGt(0); -const GF3D2 &local_dttrGt2 = gf_dttrGt(1); -const GF3D2 &local_dttrGt3 = gf_dttrGt(2); -const GF3D2 &local_dtTheta = gf_dtTheta; -const GF3D2 &local_dtalpha = gf_dtalpha; -const GF3D2 &local_dtbeta1 = gf_dtbeta(0); -const GF3D2 &local_dtbeta2 = gf_dtbeta(1); -const GF3D2 &local_dtbeta3 = gf_dtbeta(2); - -noinline([&]() __attribute__((__flatten__, __hot__)) { - grid.loop_int_device<0, 0, 0, vsize>( - grid.nghostzones, [=] ARITH_DEVICE(const PointDesc &p) ARITH_INLINE { - const vbool mask = mask_for_loop_tail(p.i, p.imax); - const GF3D2index index2(layout2, p.I); - const GF3D5index index5(layout5, p.I); - -const auto &tmp_eTtt = gf_eTtt(mask, index2); -const auto &tmp_eTt = gf_eTt(mask, index2); -const auto &tmp_eT = gf_eT(mask, index2); -const auto &tmp_chi = tl_chi(mask, index5); -const auto &tmp_gamt = tl_gamt(mask, index5); -const auto &tmp_exKh = tl_exKh(mask, index5); -const auto &tmp_exAt = tl_exAt(mask, index5); -const auto &tmp_trGt = tl_trGt(mask, index5); -const auto &tmp_Theta = tl_Theta(mask, index5); -const auto &tmp_alpha = tl_alpha(mask, index5); -const auto &tmp_beta = tl_beta(mask, index5); -const auto &tmp_dchi = tl_dchi(mask, index5); -const auto &tmp_dgamt = tl_dgamt(mask, index5); -const auto &tmp_dexKh = tl_dexKh(mask, index5); -const auto &tmp_dtrGt = tl_dtrGt(mask, index5); -const auto &tmp_dTheta = tl_dTheta(mask, index5); -const auto &tmp_dalpha = tl_dalpha(mask, index5); -const auto &tmp_dbeta = tl_dbeta(mask, index5); -const auto &tmp_ddchi = tl_ddchi(mask, index5); -const auto &tmp_ddgamt = tl_ddgamt(mask, index5); -const auto &tmp_ddalpha = tl_ddalpha(mask, index5); -const auto &tmp_ddbeta = tl_ddbeta(mask, index5); - -const vreal eTtt = tmp_eTtt; -const vreal eTt1 = tmp_eTt(0); -const vreal eTt2 = tmp_eTt(1); -const vreal eTt3 = tmp_eTt(2); -const vreal eT11 = tmp_eT(0,0); -const vreal eT12 = tmp_eT(0,1); -const vreal eT13 = tmp_eT(0,2); -const vreal eT22 = tmp_eT(1,1); -const vreal eT23 = tmp_eT(1,2); -const vreal eT33 = tmp_eT(2,2); -const vreal chi = tmp_chi + 1; -const vreal gamt11 = tmp_gamt(0,0) + 1; -const vreal gamt12 = tmp_gamt(0,1); -const vreal gamt13 = tmp_gamt(0,2); -const vreal gamt22 = tmp_gamt(1,1) + 1; -const vreal gamt23 = tmp_gamt(1,2); -const vreal gamt33 = tmp_gamt(2,2) + 1; -const vreal exKh = tmp_exKh; -const vreal exAt11 = tmp_exAt(0,0); -const vreal exAt12 = tmp_exAt(0,1); -const vreal exAt13 = tmp_exAt(0,2); -const vreal exAt22 = tmp_exAt(1,1); -const vreal exAt23 = tmp_exAt(1,2); -const vreal exAt33 = tmp_exAt(2,2); -const vreal trGt1 = tmp_trGt(0); -const vreal trGt2 = tmp_trGt(1); -const vreal trGt3 = tmp_trGt(2); -const vreal Theta = tmp_Theta; -const vreal alpha = tmp_alpha + 1; -const vreal beta1 = tmp_beta(0); -const vreal beta2 = tmp_beta(1); -const vreal beta3 = tmp_beta(2); -const vreal dchi1 = tmp_dchi(0); -const vreal dchi2 = tmp_dchi(1); -const vreal dchi3 = tmp_dchi(2); -const vreal dgamt111 = tmp_dgamt(0,0)(0); -const vreal dgamt112 = tmp_dgamt(0,1)(0); -const vreal dgamt113 = tmp_dgamt(0,2)(0); -const vreal dgamt122 = tmp_dgamt(1,1)(0); -const vreal dgamt123 = tmp_dgamt(1,2)(0); -const vreal dgamt133 = tmp_dgamt(2,2)(0); -const vreal dgamt211 = tmp_dgamt(0,0)(1); -const vreal dgamt212 = tmp_dgamt(0,1)(1); -const vreal dgamt213 = tmp_dgamt(0,2)(1); -const vreal dgamt222 = tmp_dgamt(1,1)(1); -const vreal dgamt223 = tmp_dgamt(1,2)(1); -const vreal dgamt233 = tmp_dgamt(2,2)(1); -const vreal dgamt311 = tmp_dgamt(0,0)(2); -const vreal dgamt312 = tmp_dgamt(0,1)(2); -const vreal dgamt313 = tmp_dgamt(0,2)(2); -const vreal dgamt322 = tmp_dgamt(1,1)(2); -const vreal dgamt323 = tmp_dgamt(1,2)(2); -const vreal dgamt333 = tmp_dgamt(2,2)(2); -const vreal dexKh1 = tmp_dexKh(0); -const vreal dexKh2 = tmp_dexKh(1); -const vreal dexKh3 = tmp_dexKh(2); -const vreal dtrGt11 = tmp_dtrGt(0)(0); -const vreal dtrGt12 = tmp_dtrGt(1)(0); -const vreal dtrGt13 = tmp_dtrGt(2)(0); -const vreal dtrGt21 = tmp_dtrGt(0)(1); -const vreal dtrGt22 = tmp_dtrGt(1)(1); -const vreal dtrGt23 = tmp_dtrGt(2)(1); -const vreal dtrGt31 = tmp_dtrGt(0)(2); -const vreal dtrGt32 = tmp_dtrGt(1)(2); -const vreal dtrGt33 = tmp_dtrGt(2)(2); -const vreal dTheta1 = tmp_dTheta(0); -const vreal dTheta2 = tmp_dTheta(1); -const vreal dTheta3 = tmp_dTheta(2); -const vreal dalpha1 = tmp_dalpha(0); -const vreal dalpha2 = tmp_dalpha(1); -const vreal dalpha3 = tmp_dalpha(2); -const vreal dbeta11 = tmp_dbeta(0)(0); -const vreal dbeta12 = tmp_dbeta(1)(0); -const vreal dbeta13 = tmp_dbeta(2)(0); -const vreal dbeta21 = tmp_dbeta(0)(1); -const vreal dbeta22 = tmp_dbeta(1)(1); -const vreal dbeta23 = tmp_dbeta(2)(1); -const vreal dbeta31 = tmp_dbeta(0)(2); -const vreal dbeta32 = tmp_dbeta(1)(2); -const vreal dbeta33 = tmp_dbeta(2)(2); -const vreal ddchi11 = tmp_ddchi(0,0); -const vreal ddchi12 = tmp_ddchi(0,1); -const vreal ddchi13 = tmp_ddchi(0,2); -const vreal ddchi22 = tmp_ddchi(1,1); -const vreal ddchi23 = tmp_ddchi(1,2); -const vreal ddchi33 = tmp_ddchi(2,2); -const vreal ddgamt1111 = tmp_ddgamt(0,0)(0,0); -const vreal ddgamt1112 = tmp_ddgamt(0,1)(0,0); -const vreal ddgamt1113 = tmp_ddgamt(0,2)(0,0); -const vreal ddgamt1122 = tmp_ddgamt(1,1)(0,0); -const vreal ddgamt1123 = tmp_ddgamt(1,2)(0,0); -const vreal ddgamt1133 = tmp_ddgamt(2,2)(0,0); -const vreal ddgamt1211 = tmp_ddgamt(0,0)(0,1); -const vreal ddgamt1212 = tmp_ddgamt(0,1)(0,1); -const vreal ddgamt1213 = tmp_ddgamt(0,2)(0,1); -const vreal ddgamt1222 = tmp_ddgamt(1,1)(0,1); -const vreal ddgamt1223 = tmp_ddgamt(1,2)(0,1); -const vreal ddgamt1233 = tmp_ddgamt(2,2)(0,1); -const vreal ddgamt1311 = tmp_ddgamt(0,0)(0,2); -const vreal ddgamt1312 = tmp_ddgamt(0,1)(0,2); -const vreal ddgamt1313 = tmp_ddgamt(0,2)(0,2); -const vreal ddgamt1322 = tmp_ddgamt(1,1)(0,2); -const vreal ddgamt1323 = tmp_ddgamt(1,2)(0,2); -const vreal ddgamt1333 = tmp_ddgamt(2,2)(0,2); -const vreal ddgamt2211 = tmp_ddgamt(0,0)(1,1); -const vreal ddgamt2212 = tmp_ddgamt(0,1)(1,1); -const vreal ddgamt2213 = tmp_ddgamt(0,2)(1,1); -const vreal ddgamt2222 = tmp_ddgamt(1,1)(1,1); -const vreal ddgamt2223 = tmp_ddgamt(1,2)(1,1); -const vreal ddgamt2233 = tmp_ddgamt(2,2)(1,1); -const vreal ddgamt2311 = tmp_ddgamt(0,0)(1,2); -const vreal ddgamt2312 = tmp_ddgamt(0,1)(1,2); -const vreal ddgamt2313 = tmp_ddgamt(0,2)(1,2); -const vreal ddgamt2322 = tmp_ddgamt(1,1)(1,2); -const vreal ddgamt2323 = tmp_ddgamt(1,2)(1,2); -const vreal ddgamt2333 = tmp_ddgamt(2,2)(1,2); -const vreal ddgamt3311 = tmp_ddgamt(0,0)(2,2); -const vreal ddgamt3312 = tmp_ddgamt(0,1)(2,2); -const vreal ddgamt3313 = tmp_ddgamt(0,2)(2,2); -const vreal ddgamt3322 = tmp_ddgamt(1,1)(2,2); -const vreal ddgamt3323 = tmp_ddgamt(1,2)(2,2); -const vreal ddgamt3333 = tmp_ddgamt(2,2)(2,2); -const vreal ddalpha11 = tmp_ddalpha(0,0); -const vreal ddalpha12 = tmp_ddalpha(0,1); -const vreal ddalpha13 = tmp_ddalpha(0,2); -const vreal ddalpha22 = tmp_ddalpha(1,1); -const vreal ddalpha23 = tmp_ddalpha(1,2); -const vreal ddalpha33 = tmp_ddalpha(2,2); -const vreal ddbeta111 = tmp_ddbeta(0)(0,0); -const vreal ddbeta112 = tmp_ddbeta(1)(0,0); -const vreal ddbeta113 = tmp_ddbeta(2)(0,0); -const vreal ddbeta121 = tmp_ddbeta(0)(0,1); -const vreal ddbeta122 = tmp_ddbeta(1)(0,1); -const vreal ddbeta123 = tmp_ddbeta(2)(0,1); -const vreal ddbeta131 = tmp_ddbeta(0)(0,2); -const vreal ddbeta132 = tmp_ddbeta(1)(0,2); -const vreal ddbeta133 = tmp_ddbeta(2)(0,2); -const vreal ddbeta221 = tmp_ddbeta(0)(1,1); -const vreal ddbeta222 = tmp_ddbeta(1)(1,1); -const vreal ddbeta223 = tmp_ddbeta(2)(1,1); -const vreal ddbeta231 = tmp_ddbeta(0)(1,2); -const vreal ddbeta232 = tmp_ddbeta(1)(1,2); -const vreal ddbeta233 = tmp_ddbeta(2)(1,2); -const vreal ddbeta331 = tmp_ddbeta(0)(2,2); -const vreal ddbeta332 = tmp_ddbeta(1)(2,2); -const vreal ddbeta333 = tmp_ddbeta(2)(2,2); - -vreal invgamt11 -= --Power(gamt23,2) + gamt22*gamt33 -; - -vreal invgamt12 -= -gamt13*gamt23 - gamt12*gamt33 -; - -vreal invgamt13 -= --(gamt13*gamt22) + gamt12*gamt23 -; - -vreal invgamt22 -= --Power(gamt13,2) + gamt11*gamt33 -; - -vreal invgamt23 -= -gamt12*gamt13 - gamt11*gamt23 -; - -vreal invgamt33 -= --Power(gamt12,2) + gamt11*gamt22 -; - -vreal invgam11 -= -chi*invgamt11 -; - -vreal invgam12 -= -chi*invgamt12 -; - -vreal invgam13 -= -chi*invgamt13 -; - -vreal invgam22 -= -chi*invgamt22 -; - -vreal invgam23 -= -chi*invgamt23 -; - -vreal invgam33 -= -chi*invgamt33 -; - -vreal gam11 -= -gamt11/chi -; - -vreal gam12 -= -gamt12/chi -; - -vreal gam13 -= -gamt13/chi -; - -vreal gam22 -= -gamt22/chi -; - -vreal gam23 -= -gamt23/chi -; - -vreal gam33 -= -gamt33/chi -; - -vreal GtDDD111 -= -dgamt111/2. -; - -vreal GtDDD112 -= -dgamt211/2. -; - -vreal GtDDD113 -= -dgamt311/2. -; - -vreal GtDDD122 -= --0.5*dgamt122 + dgamt212 -; - -vreal GtDDD123 -= -(-dgamt123 + dgamt213 + dgamt312)/2. -; - -vreal GtDDD133 -= --0.5*dgamt133 + dgamt313 -; - -vreal GtDDD211 -= -dgamt112 - dgamt211/2. -; - -vreal GtDDD212 -= -dgamt122/2. -; - -vreal GtDDD213 -= -(dgamt123 - dgamt213 + dgamt312)/2. -; - -vreal GtDDD222 -= -dgamt222/2. -; - -vreal GtDDD223 -= -dgamt322/2. -; - -vreal GtDDD233 -= --0.5*dgamt233 + dgamt323 -; - -vreal GtDDD311 -= -dgamt113 - dgamt311/2. -; - -vreal GtDDD312 -= -(dgamt123 + dgamt213 - dgamt312)/2. -; - -vreal GtDDD313 -= -dgamt133/2. -; - -vreal GtDDD322 -= -dgamt223 - dgamt322/2. -; - -vreal GtDDD323 -= -dgamt233/2. -; - -vreal GtDDD333 -= -dgamt333/2. -; - -vreal GtDDU111 -= -GtDDD111*invgamt11 + GtDDD112*invgamt12 + GtDDD113*invgamt13 -; - -vreal GtDDU112 -= -GtDDD111*invgamt12 + GtDDD112*invgamt22 + GtDDD113*invgamt23 -; - -vreal GtDDU113 -= -GtDDD111*invgamt13 + GtDDD112*invgamt23 + GtDDD113*invgamt33 -; - -vreal GtDDU121 -= -GtDDD112*invgamt11 + GtDDD122*invgamt12 + GtDDD123*invgamt13 -; - -vreal GtDDU122 -= -GtDDD112*invgamt12 + GtDDD122*invgamt22 + GtDDD123*invgamt23 -; - -vreal GtDDU123 -= -GtDDD112*invgamt13 + GtDDD122*invgamt23 + GtDDD123*invgamt33 -; - -vreal GtDDU131 -= -GtDDD113*invgamt11 + GtDDD123*invgamt12 + GtDDD133*invgamt13 -; - -vreal GtDDU132 -= -GtDDD113*invgamt12 + GtDDD123*invgamt22 + GtDDD133*invgamt23 -; - -vreal GtDDU133 -= -GtDDD113*invgamt13 + GtDDD123*invgamt23 + GtDDD133*invgamt33 -; - -vreal GtDDU211 -= -GtDDD211*invgamt11 + GtDDD212*invgamt12 + GtDDD213*invgamt13 -; - -vreal GtDDU212 -= -GtDDD211*invgamt12 + GtDDD212*invgamt22 + GtDDD213*invgamt23 -; - -vreal GtDDU213 -= -GtDDD211*invgamt13 + GtDDD212*invgamt23 + GtDDD213*invgamt33 -; - -vreal GtDDU221 -= -GtDDD212*invgamt11 + GtDDD222*invgamt12 + GtDDD223*invgamt13 -; - -vreal GtDDU222 -= -GtDDD212*invgamt12 + GtDDD222*invgamt22 + GtDDD223*invgamt23 -; - -vreal GtDDU223 -= -GtDDD212*invgamt13 + GtDDD222*invgamt23 + GtDDD223*invgamt33 -; - -vreal GtDDU231 -= -GtDDD213*invgamt11 + GtDDD223*invgamt12 + GtDDD233*invgamt13 -; - -vreal GtDDU232 -= -GtDDD213*invgamt12 + GtDDD223*invgamt22 + GtDDD233*invgamt23 -; - -vreal GtDDU233 -= -GtDDD213*invgamt13 + GtDDD223*invgamt23 + GtDDD233*invgamt33 -; - -vreal GtDDU311 -= -GtDDD311*invgamt11 + GtDDD312*invgamt12 + GtDDD313*invgamt13 -; - -vreal GtDDU312 -= -GtDDD311*invgamt12 + GtDDD312*invgamt22 + GtDDD313*invgamt23 -; - -vreal GtDDU313 -= -GtDDD311*invgamt13 + GtDDD312*invgamt23 + GtDDD313*invgamt33 -; - -vreal GtDDU321 -= -GtDDD312*invgamt11 + GtDDD322*invgamt12 + GtDDD323*invgamt13 -; - -vreal GtDDU322 -= -GtDDD312*invgamt12 + GtDDD322*invgamt22 + GtDDD323*invgamt23 -; - -vreal GtDDU323 -= -GtDDD312*invgamt13 + GtDDD322*invgamt23 + GtDDD323*invgamt33 -; - -vreal GtDDU331 -= -GtDDD313*invgamt11 + GtDDD323*invgamt12 + GtDDD333*invgamt13 -; - -vreal GtDDU332 -= -GtDDD313*invgamt12 + GtDDD323*invgamt22 + GtDDD333*invgamt23 -; - -vreal GtDDU333 -= -GtDDD313*invgamt13 + GtDDD323*invgamt23 + GtDDD333*invgamt33 -; - -vreal Gt111 -= -GtDDD111*invgamt11 + GtDDD211*invgamt12 + GtDDD311*invgamt13 -; - -vreal Gt112 -= -GtDDD112*invgamt11 + GtDDD212*invgamt12 + GtDDD312*invgamt13 -; - -vreal Gt113 -= -GtDDD113*invgamt11 + GtDDD213*invgamt12 + GtDDD313*invgamt13 -; - -vreal Gt122 -= -GtDDD122*invgamt11 + GtDDD222*invgamt12 + GtDDD322*invgamt13 -; - -vreal Gt123 -= -GtDDD123*invgamt11 + GtDDD223*invgamt12 + GtDDD323*invgamt13 -; - -vreal Gt133 -= -GtDDD133*invgamt11 + GtDDD233*invgamt12 + GtDDD333*invgamt13 -; - -vreal Gt211 -= -GtDDD111*invgamt12 + GtDDD211*invgamt22 + GtDDD311*invgamt23 -; - -vreal Gt212 -= -GtDDD112*invgamt12 + GtDDD212*invgamt22 + GtDDD312*invgamt23 -; - -vreal Gt213 -= -GtDDD113*invgamt12 + GtDDD213*invgamt22 + GtDDD313*invgamt23 -; - -vreal Gt222 -= -GtDDD122*invgamt12 + GtDDD222*invgamt22 + GtDDD322*invgamt23 -; - -vreal Gt223 -= -GtDDD123*invgamt12 + GtDDD223*invgamt22 + GtDDD323*invgamt23 -; - -vreal Gt233 -= -GtDDD133*invgamt12 + GtDDD233*invgamt22 + GtDDD333*invgamt23 -; - -vreal Gt311 -= -GtDDD111*invgamt13 + GtDDD211*invgamt23 + GtDDD311*invgamt33 -; - -vreal Gt312 -= -GtDDD112*invgamt13 + GtDDD212*invgamt23 + GtDDD312*invgamt33 -; - -vreal Gt313 -= -GtDDD113*invgamt13 + GtDDD213*invgamt23 + GtDDD313*invgamt33 -; - -vreal Gt322 -= -GtDDD122*invgamt13 + GtDDD222*invgamt23 + GtDDD322*invgamt33 -; - -vreal Gt323 -= -GtDDD123*invgamt13 + GtDDD223*invgamt23 + GtDDD323*invgamt33 -; - -vreal Gt333 -= -GtDDD133*invgamt13 + GtDDD233*invgamt23 + GtDDD333*invgamt33 -; - -vreal trGtd1 -= -Gt111*invgamt11 + 2*Gt112*invgamt12 + 2*Gt113*invgamt13 + Gt122*invgamt22 + - 2*Gt123*invgamt23 + Gt133*invgamt33 -; - -vreal trGtd2 -= -Gt211*invgamt11 + 2*Gt212*invgamt12 + 2*Gt213*invgamt13 + Gt222*invgamt22 + - 2*Gt223*invgamt23 + Gt233*invgamt33 -; - -vreal trGtd3 -= -Gt311*invgamt11 + 2*Gt312*invgamt12 + 2*Gt313*invgamt13 + Gt322*invgamt22 + - 2*Gt323*invgamt23 + Gt333*invgamt33 -; - -vreal dgam111 -= -(chi*dgamt111 - dchi1*gamt11)/Power(chi,2) -; - -vreal dgam112 -= -(chi*dgamt112 - dchi1*gamt12)/Power(chi,2) -; - -vreal dgam113 -= -(chi*dgamt113 - dchi1*gamt13)/Power(chi,2) -; - -vreal dgam122 -= -(chi*dgamt122 - dchi1*gamt22)/Power(chi,2) -; - -vreal dgam123 -= -(chi*dgamt123 - dchi1*gamt23)/Power(chi,2) -; - -vreal dgam133 -= -(chi*dgamt133 - dchi1*gamt33)/Power(chi,2) -; - -vreal dgam211 -= -(chi*dgamt211 - dchi2*gamt11)/Power(chi,2) -; - -vreal dgam212 -= -(chi*dgamt212 - dchi2*gamt12)/Power(chi,2) -; - -vreal dgam213 -= -(chi*dgamt213 - dchi2*gamt13)/Power(chi,2) -; - -vreal dgam222 -= -(chi*dgamt222 - dchi2*gamt22)/Power(chi,2) -; - -vreal dgam223 -= -(chi*dgamt223 - dchi2*gamt23)/Power(chi,2) -; - -vreal dgam233 -= -(chi*dgamt233 - dchi2*gamt33)/Power(chi,2) -; - -vreal dgam311 -= -(chi*dgamt311 - dchi3*gamt11)/Power(chi,2) -; - -vreal dgam312 -= -(chi*dgamt312 - dchi3*gamt12)/Power(chi,2) -; - -vreal dgam313 -= -(chi*dgamt313 - dchi3*gamt13)/Power(chi,2) -; - -vreal dgam322 -= -(chi*dgamt322 - dchi3*gamt22)/Power(chi,2) -; - -vreal dgam323 -= -(chi*dgamt323 - dchi3*gamt23)/Power(chi,2) -; - -vreal dgam333 -= -(chi*dgamt333 - dchi3*gamt33)/Power(chi,2) -; - -vreal GamDDD111 -= -dgam111/2. -; - -vreal GamDDD112 -= -dgam211/2. -; - -vreal GamDDD113 -= -dgam311/2. -; - -vreal GamDDD122 -= --0.5*dgam122 + dgam212 -; - -vreal GamDDD123 -= -(-dgam123 + dgam213 + dgam312)/2. -; - -vreal GamDDD133 -= --0.5*dgam133 + dgam313 -; - -vreal GamDDD211 -= -dgam112 - dgam211/2. -; - -vreal GamDDD212 -= -dgam122/2. -; - -vreal GamDDD213 -= -(dgam123 - dgam213 + dgam312)/2. -; - -vreal GamDDD222 -= -dgam222/2. -; - -vreal GamDDD223 -= -dgam322/2. -; - -vreal GamDDD233 -= --0.5*dgam233 + dgam323 -; - -vreal GamDDD311 -= -dgam113 - dgam311/2. -; - -vreal GamDDD312 -= -(dgam123 + dgam213 - dgam312)/2. -; - -vreal GamDDD313 -= -dgam133/2. -; - -vreal GamDDD322 -= -dgam223 - dgam322/2. -; - -vreal GamDDD323 -= -dgam233/2. -; - -vreal GamDDD333 -= -dgam333/2. -; - -vreal Gam111 -= -GamDDD111*invgam11 + GamDDD211*invgam12 + GamDDD311*invgam13 -; - -vreal Gam112 -= -GamDDD112*invgam11 + GamDDD212*invgam12 + GamDDD312*invgam13 -; - -vreal Gam113 -= -GamDDD113*invgam11 + GamDDD213*invgam12 + GamDDD313*invgam13 -; - -vreal Gam122 -= -GamDDD122*invgam11 + GamDDD222*invgam12 + GamDDD322*invgam13 -; - -vreal Gam123 -= -GamDDD123*invgam11 + GamDDD223*invgam12 + GamDDD323*invgam13 -; - -vreal Gam133 -= -GamDDD133*invgam11 + GamDDD233*invgam12 + GamDDD333*invgam13 -; - -vreal Gam211 -= -GamDDD111*invgam12 + GamDDD211*invgam22 + GamDDD311*invgam23 -; - -vreal Gam212 -= -GamDDD112*invgam12 + GamDDD212*invgam22 + GamDDD312*invgam23 -; - -vreal Gam213 -= -GamDDD113*invgam12 + GamDDD213*invgam22 + GamDDD313*invgam23 -; - -vreal Gam222 -= -GamDDD122*invgam12 + GamDDD222*invgam22 + GamDDD322*invgam23 -; - -vreal Gam223 -= -GamDDD123*invgam12 + GamDDD223*invgam22 + GamDDD323*invgam23 -; - -vreal Gam233 -= -GamDDD133*invgam12 + GamDDD233*invgam22 + GamDDD333*invgam23 -; - -vreal Gam311 -= -GamDDD111*invgam13 + GamDDD211*invgam23 + GamDDD311*invgam33 -; - -vreal Gam312 -= -GamDDD112*invgam13 + GamDDD212*invgam23 + GamDDD312*invgam33 -; - -vreal Gam313 -= -GamDDD113*invgam13 + GamDDD213*invgam23 + GamDDD313*invgam33 -; - -vreal Gam322 -= -GamDDD122*invgam13 + GamDDD222*invgam23 + GamDDD322*invgam33 -; - -vreal Gam323 -= -GamDDD123*invgam13 + GamDDD223*invgam23 + GamDDD323*invgam33 -; - -vreal Gam333 -= -GamDDD133*invgam13 + GamDDD233*invgam23 + GamDDD333*invgam33 -; - -vreal exAtUU11 -= -exAt11*Power(invgamt11,2) + 2*exAt12*invgamt11*invgamt12 + - exAt22*Power(invgamt12,2) + 2*exAt13*invgamt11*invgamt13 + - 2*exAt23*invgamt12*invgamt13 + exAt33*Power(invgamt13,2) -; - -vreal exAtUU12 -= -exAt11*invgamt11*invgamt12 + exAt13*invgamt12*invgamt13 + - exAt22*invgamt12*invgamt22 + exAt23*invgamt13*invgamt22 + - exAt12*(Power(invgamt12,2) + invgamt11*invgamt22) + - exAt13*invgamt11*invgamt23 + exAt23*invgamt12*invgamt23 + - exAt33*invgamt13*invgamt23 -; - -vreal exAtUU13 -= -exAt11*invgamt11*invgamt13 + exAt12*invgamt12*invgamt13 + - exAt13*Power(invgamt13,2) + exAt12*invgamt11*invgamt23 + - exAt22*invgamt12*invgamt23 + exAt23*invgamt13*invgamt23 + - exAt13*invgamt11*invgamt33 + exAt23*invgamt12*invgamt33 + - exAt33*invgamt13*invgamt33 -; - -vreal exAtUU22 -= -exAt11*Power(invgamt12,2) + 2*exAt12*invgamt12*invgamt22 + - exAt22*Power(invgamt22,2) + 2*exAt13*invgamt12*invgamt23 + - 2*exAt23*invgamt22*invgamt23 + exAt33*Power(invgamt23,2) -; - -vreal exAtUU23 -= -exAt11*invgamt12*invgamt13 + exAt12*invgamt13*invgamt22 + - exAt12*invgamt12*invgamt23 + exAt13*invgamt13*invgamt23 + - exAt22*invgamt22*invgamt23 + exAt23*Power(invgamt23,2) + - exAt13*invgamt12*invgamt33 + exAt23*invgamt22*invgamt33 + - exAt33*invgamt23*invgamt33 -; - -vreal exAtUU33 -= -exAt11*Power(invgamt13,2) + 2*exAt12*invgamt13*invgamt23 + - exAt22*Power(invgamt23,2) + 2*exAt13*invgamt13*invgamt33 + - 2*exAt23*invgamt23*invgamt33 + exAt33*Power(invgamt33,2) -; - -vreal tDtDchi11 -= -ddchi11 - dchi1*Gt111 - dchi2*Gt211 - dchi3*Gt311 -; - -vreal tDtDchi12 -= -ddchi12 - dchi1*Gt112 - dchi2*Gt212 - dchi3*Gt312 -; - -vreal tDtDchi13 -= -ddchi13 - dchi1*Gt113 - dchi2*Gt213 - dchi3*Gt313 -; - -vreal tDtDchi22 -= -ddchi22 - dchi1*Gt122 - dchi2*Gt222 - dchi3*Gt322 -; - -vreal tDtDchi23 -= -ddchi23 - dchi1*Gt123 - dchi2*Gt223 - dchi3*Gt323 -; - -vreal tDtDchi33 -= -ddchi33 - dchi1*Gt133 - dchi2*Gt233 - dchi3*Gt333 -; - -vreal DDalpha11 -= -ddalpha11 - dalpha1*Gam111 - dalpha2*Gam211 - dalpha3*Gam311 -; - -vreal DDalpha12 -= -ddalpha12 - dalpha1*Gam112 - dalpha2*Gam212 - dalpha3*Gam312 -; - -vreal DDalpha13 -= -ddalpha13 - dalpha1*Gam113 - dalpha2*Gam213 - dalpha3*Gam313 -; - -vreal DDalpha22 -= -ddalpha22 - dalpha1*Gam122 - dalpha2*Gam222 - dalpha3*Gam322 -; - -vreal DDalpha23 -= -ddalpha23 - dalpha1*Gam123 - dalpha2*Gam223 - dalpha3*Gam323 -; - -vreal DDalpha33 -= -ddalpha33 - dalpha1*Gam133 - dalpha2*Gam233 - dalpha3*Gam333 -; - -vreal Rtchi11 -= -(-(Power(dchi1,2)*(1 + 3*gamt11*invgamt11)) - - 6*dchi1*gamt11*(dchi2*invgamt12 + dchi3*invgamt13) - - 3*Power(dchi2,2)*gamt11*invgamt22 - 6*dchi2*dchi3*gamt11*invgamt23 - - 3*Power(dchi3,2)*gamt11*invgamt33 + 2*chi*tDtDchi11 + - 2*chi*gamt11*invgamt11*tDtDchi11 + 4*chi*gamt11*invgamt12*tDtDchi12 + - 4*chi*gamt11*invgamt13*tDtDchi13 + 2*chi*gamt11*invgamt22*tDtDchi22 + - 4*chi*gamt11*invgamt23*tDtDchi23 + 2*chi*gamt11*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi12 -= -(-3*Power(dchi1,2)*gamt12*invgamt11 - - dchi1*(dchi2 + 6*dchi2*gamt12*invgamt12 + 6*dchi3*gamt12*invgamt13) - - 3*Power(dchi2,2)*gamt12*invgamt22 - 6*dchi2*dchi3*gamt12*invgamt23 - - 3*Power(dchi3,2)*gamt12*invgamt33 + 2*chi*gamt12*invgamt11*tDtDchi11 + - 2*chi*tDtDchi12 + 4*chi*gamt12*invgamt12*tDtDchi12 + - 4*chi*gamt12*invgamt13*tDtDchi13 + 2*chi*gamt12*invgamt22*tDtDchi22 + - 4*chi*gamt12*invgamt23*tDtDchi23 + 2*chi*gamt12*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi13 -= -(-3*Power(dchi1,2)*gamt13*invgamt11 - - dchi1*(dchi3 + 6*dchi2*gamt13*invgamt12 + 6*dchi3*gamt13*invgamt13) - - 3*Power(dchi2,2)*gamt13*invgamt22 - 6*dchi2*dchi3*gamt13*invgamt23 - - 3*Power(dchi3,2)*gamt13*invgamt33 + 2*chi*gamt13*invgamt11*tDtDchi11 + - 4*chi*gamt13*invgamt12*tDtDchi12 + 2*chi*tDtDchi13 + - 4*chi*gamt13*invgamt13*tDtDchi13 + 2*chi*gamt13*invgamt22*tDtDchi22 + - 4*chi*gamt13*invgamt23*tDtDchi23 + 2*chi*gamt13*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi22 -= -(-3*Power(dchi1,2)*gamt22*invgamt11 - 6*dchi1*dchi3*gamt22*invgamt13 - - Power(dchi2,2)*(1 + 3*gamt22*invgamt22) - - 6*dchi2*gamt22*(dchi1*invgamt12 + dchi3*invgamt23) - - 3*Power(dchi3,2)*gamt22*invgamt33 + 2*chi*gamt22*invgamt11*tDtDchi11 + - 4*chi*gamt22*invgamt12*tDtDchi12 + 4*chi*gamt22*invgamt13*tDtDchi13 + - 2*chi*tDtDchi22 + 2*chi*gamt22*invgamt22*tDtDchi22 + - 4*chi*gamt22*invgamt23*tDtDchi23 + 2*chi*gamt22*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi23 -= -(-3*Power(dchi1,2)*gamt23*invgamt11 - 6*dchi1*dchi3*gamt23*invgamt13 - - 3*Power(dchi2,2)*gamt23*invgamt22 - - dchi2*(dchi3 + 6*dchi1*gamt23*invgamt12 + 6*dchi3*gamt23*invgamt23) - - 3*Power(dchi3,2)*gamt23*invgamt33 + 2*chi*gamt23*invgamt11*tDtDchi11 + - 4*chi*gamt23*invgamt12*tDtDchi12 + 4*chi*gamt23*invgamt13*tDtDchi13 + - 2*chi*gamt23*invgamt22*tDtDchi22 + 2*chi*tDtDchi23 + - 4*chi*gamt23*invgamt23*tDtDchi23 + 2*chi*gamt23*invgamt33*tDtDchi33)/ - (4.*Power(chi,2)) -; - -vreal Rtchi33 -= -(-3*Power(dchi1,2)*gamt33*invgamt11 - 6*dchi1*dchi2*gamt33*invgamt12 - - 3*Power(dchi2,2)*gamt33*invgamt22 - - 6*dchi3*gamt33*(dchi1*invgamt13 + dchi2*invgamt23) - - Power(dchi3,2)*(1 + 3*gamt33*invgamt33) + - 2*chi*gamt33*invgamt11*tDtDchi11 + 4*chi*gamt33*invgamt12*tDtDchi12 + - 4*chi*gamt33*invgamt13*tDtDchi13 + 2*chi*gamt33*invgamt22*tDtDchi22 + - 4*chi*gamt33*invgamt23*tDtDchi23 + 2*chi*tDtDchi33 + - 2*chi*gamt33*invgamt33*tDtDchi33)/(4.*Power(chi,2)) -; - -vreal Rt11 -= -dtrGt11*gamt11 + dtrGt12*gamt12 + dtrGt13*gamt13 + 3*Gt111*GtDDU111 + - 3*Gt112*GtDDU112 + 3*Gt113*GtDDU113 + 2*Gt211*GtDDU121 + - 2*Gt212*GtDDU122 + 2*Gt213*GtDDU123 + 2*Gt311*GtDDU131 + - 2*Gt312*GtDDU132 + 2*Gt313*GtDDU133 + Gt211*GtDDU211 + Gt212*GtDDU212 + - Gt213*GtDDU213 + Gt311*GtDDU311 + Gt312*GtDDU312 + Gt313*GtDDU313 - - (ddgamt1111*invgamt11)/2. - ddgamt1211*invgamt12 - ddgamt1311*invgamt13 - - (ddgamt2211*invgamt22)/2. - ddgamt2311*invgamt23 - - (ddgamt3311*invgamt33)/2. + GtDDD111*trGtd1 + GtDDD112*trGtd2 + - GtDDD113*trGtd3 -; - -vreal Rt12 -= -(dtrGt21*gamt11 + dtrGt11*gamt12 + dtrGt22*gamt12 + dtrGt23*gamt13 + - dtrGt12*gamt22 + dtrGt13*gamt23 + 2*Gt112*GtDDU111 + 2*Gt122*GtDDU112 + - 2*Gt123*GtDDU113 + 2*Gt111*GtDDU121 + 2*Gt212*GtDDU121 + - 2*Gt112*GtDDU122 + 2*Gt222*GtDDU122 + 2*Gt113*GtDDU123 + - 2*Gt223*GtDDU123 + 2*Gt312*GtDDU131 + 2*Gt322*GtDDU132 + - 2*Gt323*GtDDU133 + 2*Gt111*GtDDU211 + 2*Gt112*GtDDU212 + - 2*Gt113*GtDDU213 + 4*Gt211*GtDDU221 + 4*Gt212*GtDDU222 + - 4*Gt213*GtDDU223 + 2*Gt311*GtDDU231 + 2*Gt312*GtDDU232 + - 2*Gt313*GtDDU233 + 2*Gt311*GtDDU321 + 2*Gt312*GtDDU322 + - 2*Gt313*GtDDU323 - ddgamt1112*invgamt11 - 2*ddgamt1212*invgamt12 - - 2*ddgamt1312*invgamt13 - ddgamt2212*invgamt22 - - 2*ddgamt2312*invgamt23 - ddgamt3312*invgamt33 + GtDDD112*trGtd1 + - GtDDD211*trGtd1 + GtDDD122*trGtd2 + GtDDD212*trGtd2 + GtDDD123*trGtd3 + - GtDDD213*trGtd3)/2. -; - -vreal Rt13 -= -(dtrGt31*gamt11 + dtrGt32*gamt12 + dtrGt11*gamt13 + dtrGt33*gamt13 + - dtrGt12*gamt23 + dtrGt13*gamt33 + 2*Gt113*GtDDU111 + 2*Gt123*GtDDU112 + - 2*Gt133*GtDDU113 + 2*Gt213*GtDDU121 + 2*Gt223*GtDDU122 + - 2*Gt233*GtDDU123 + 2*Gt111*GtDDU131 + 2*Gt313*GtDDU131 + - 2*Gt112*GtDDU132 + 2*Gt323*GtDDU132 + 2*Gt113*GtDDU133 + - 2*Gt333*GtDDU133 + 2*Gt211*GtDDU231 + 2*Gt212*GtDDU232 + - 2*Gt213*GtDDU233 + 2*Gt111*GtDDU311 + 2*Gt112*GtDDU312 + - 2*Gt113*GtDDU313 + 2*Gt211*GtDDU321 + 2*Gt212*GtDDU322 + - 2*Gt213*GtDDU323 + 4*Gt311*GtDDU331 + 4*Gt312*GtDDU332 + - 4*Gt313*GtDDU333 - ddgamt1113*invgamt11 - 2*ddgamt1213*invgamt12 - - 2*ddgamt1313*invgamt13 - ddgamt2213*invgamt22 - - 2*ddgamt2313*invgamt23 - ddgamt3313*invgamt33 + GtDDD113*trGtd1 + - GtDDD311*trGtd1 + GtDDD123*trGtd2 + GtDDD312*trGtd2 + GtDDD133*trGtd3 + - GtDDD313*trGtd3)/2. -; - -vreal Rt22 -= -dtrGt21*gamt12 + dtrGt22*gamt22 + dtrGt23*gamt23 + Gt112*GtDDU121 + - Gt122*GtDDU122 + Gt123*GtDDU123 + 2*Gt112*GtDDU211 + 2*Gt122*GtDDU212 + - 2*Gt123*GtDDU213 + 3*Gt212*GtDDU221 + 3*Gt222*GtDDU222 + - 3*Gt223*GtDDU223 + 2*Gt312*GtDDU231 + 2*Gt322*GtDDU232 + - 2*Gt323*GtDDU233 + Gt312*GtDDU321 + Gt322*GtDDU322 + Gt323*GtDDU323 - - (ddgamt1122*invgamt11)/2. - ddgamt1222*invgamt12 - ddgamt1322*invgamt13 - - (ddgamt2222*invgamt22)/2. - ddgamt2322*invgamt23 - - (ddgamt3322*invgamt33)/2. + GtDDD212*trGtd1 + GtDDD222*trGtd2 + - GtDDD223*trGtd3 -; - -vreal Rt23 -= -(dtrGt31*gamt12 + dtrGt21*gamt13 + dtrGt32*gamt22 + dtrGt22*gamt23 + - dtrGt33*gamt23 + dtrGt23*gamt33 + 2*Gt112*GtDDU131 + 2*Gt122*GtDDU132 + - 2*Gt123*GtDDU133 + 2*Gt113*GtDDU211 + 2*Gt123*GtDDU212 + - 2*Gt133*GtDDU213 + 2*Gt213*GtDDU221 + 2*Gt223*GtDDU222 + - 2*Gt233*GtDDU223 + 2*Gt212*GtDDU231 + 2*Gt313*GtDDU231 + - 2*Gt222*GtDDU232 + 2*Gt323*GtDDU232 + 2*Gt223*GtDDU233 + - 2*Gt333*GtDDU233 + 2*Gt112*GtDDU311 + 2*Gt122*GtDDU312 + - 2*Gt123*GtDDU313 + 2*Gt212*GtDDU321 + 2*Gt222*GtDDU322 + - 2*Gt223*GtDDU323 + 4*Gt312*GtDDU331 + 4*Gt322*GtDDU332 + - 4*Gt323*GtDDU333 - ddgamt1123*invgamt11 - 2*ddgamt1223*invgamt12 - - 2*ddgamt1323*invgamt13 - ddgamt2223*invgamt22 - - 2*ddgamt2323*invgamt23 - ddgamt3323*invgamt33 + GtDDD213*trGtd1 + - GtDDD312*trGtd1 + GtDDD223*trGtd2 + GtDDD322*trGtd2 + GtDDD233*trGtd3 + - GtDDD323*trGtd3)/2. -; - -vreal Rt33 -= -dtrGt31*gamt13 + dtrGt32*gamt23 + dtrGt33*gamt33 + Gt113*GtDDU131 + - Gt123*GtDDU132 + Gt133*GtDDU133 + Gt213*GtDDU231 + Gt223*GtDDU232 + - Gt233*GtDDU233 + 2*Gt113*GtDDU311 + 2*Gt123*GtDDU312 + 2*Gt133*GtDDU313 + - 2*Gt213*GtDDU321 + 2*Gt223*GtDDU322 + 2*Gt233*GtDDU323 + - 3*Gt313*GtDDU331 + 3*Gt323*GtDDU332 + 3*Gt333*GtDDU333 - - (ddgamt1133*invgamt11)/2. - ddgamt1233*invgamt12 - ddgamt1333*invgamt13 - - (ddgamt2233*invgamt22)/2. - ddgamt2333*invgamt23 - - (ddgamt3333*invgamt33)/2. + GtDDD313*trGtd1 + GtDDD323*trGtd2 + - GtDDD333*trGtd3 -; - -vreal R11 -= -Rt11 + Rtchi11 -; - -vreal R12 -= -Rt12 + Rtchi12 -; - -vreal R13 -= -Rt13 + Rtchi13 -; - -vreal R22 -= -Rt22 + Rtchi22 -; - -vreal R23 -= -Rt23 + Rtchi23 -; - -vreal R33 -= -Rt33 + Rtchi33 -; - -vreal trR -= -invgam11*R11 + 2*invgam12*R12 + 2*invgam13*R13 + invgam22*R22 + - 2*invgam23*R23 + invgam33*R33 -; - -vreal rho -= -(Power(beta1,2)*eT11 + Power(beta2,2)*eT22 + 2*beta2*beta3*eT23 + - Power(beta3,2)*eT33 + 2*beta1*(beta2*eT12 + beta3*eT13 - eTt1) - - 2*beta2*eTt2 - 2*beta3*eTt3 + eTtt)/Power(alpha,2) -; - -vreal Sm1 -= -(beta1*eT11 + beta2*eT12 + beta3*eT13 - eTt1)/alpha -; - -vreal Sm2 -= -(beta1*eT12 + beta2*eT22 + beta3*eT23 - eTt2)/alpha -; - -vreal Sm3 -= -(beta1*eT13 + beta2*eT23 + beta3*eT33 - eTt3)/alpha -; - -vreal Ss11 -= -eT11 -; - -vreal Ss12 -= -eT12 -; - -vreal Ss13 -= -eT13 -; - -vreal Ss22 -= -eT22 -; - -vreal Ss23 -= -eT23 -; - -vreal Ss33 -= -eT33 -; - -vreal trSs -= -invgam11*Ss11 + 2*invgam12*Ss12 + 2*invgam13*Ss13 + invgam22*Ss22 + - 2*invgam23*Ss23 + invgam33*Ss33 -; - - -local_dtchi.store(mask, index2, -(-2*chi*(dbeta11 + dbeta22 + dbeta33 - alpha*exKh - 2*alpha*Theta))/3. -); - -local_dtgamt11.store(mask, index2, --2*alpha*exAt11 + 2*dbeta11*gamt11 - - (2*(dbeta11 + dbeta22 + dbeta33)*gamt11)/3. + 2*dbeta12*gamt12 + - 2*dbeta13*gamt13 -); - -local_dtgamt12.store(mask, index2, --2*alpha*exAt12 + dbeta21*gamt11 + dbeta11*gamt12 + dbeta22*gamt12 - - (2*(dbeta11 + dbeta22 + dbeta33)*gamt12)/3. + dbeta23*gamt13 + - dbeta12*gamt22 + dbeta13*gamt23 -); - -local_dtgamt13.store(mask, index2, --2*alpha*exAt13 + dbeta31*gamt11 + dbeta32*gamt12 + dbeta11*gamt13 + - dbeta33*gamt13 - (2*(dbeta11 + dbeta22 + dbeta33)*gamt13)/3. + - dbeta12*gamt23 + dbeta13*gamt33 -); - -local_dtgamt22.store(mask, index2, --2*alpha*exAt22 + 2*dbeta21*gamt12 + 2*dbeta22*gamt22 - - (2*(dbeta11 + dbeta22 + dbeta33)*gamt22)/3. + 2*dbeta23*gamt23 -); - -local_dtgamt23.store(mask, index2, --2*alpha*exAt23 + dbeta31*gamt12 + dbeta21*gamt13 + dbeta32*gamt22 + - dbeta22*gamt23 + dbeta33*gamt23 - - (2*(dbeta11 + dbeta22 + dbeta33)*gamt23)/3. + dbeta23*gamt33 -); - -local_dtgamt33.store(mask, index2, --2*alpha*exAt33 + 2*dbeta31*gamt13 + 2*dbeta32*gamt23 + 2*dbeta33*gamt33 - - (2*(dbeta11 + dbeta22 + dbeta33)*gamt33)/3. -); - -local_dtexKh.store(mask, index2, --(DDalpha11*invgam11) - 2*DDalpha12*invgam12 - 2*DDalpha13*invgam13 - - DDalpha22*invgam22 - 2*DDalpha23*invgam23 - DDalpha33*invgam33 + - alpha*(exAt11*exAtUU11 + 2*exAt12*exAtUU12 + 2*exAt13*exAtUU13 + - exAt22*exAtUU22 + 2*exAt23*exAtUU23 + exAt33*exAtUU33 + - Power(exKh,2)/3. + 4*cpi*rho + ckappa1*Theta - ckappa1*ckappa2*Theta + - (4*exKh*Theta)/3. + (4*Power(Theta,2))/3. + 4*cpi*trSs) -); - -local_dtexAt11.store(mask, index2, -(4*dbeta11*exAt11 - 2*dbeta22*exAt11 - 2*dbeta33*exAt11 + - 6*dbeta12*exAt12 + 6*dbeta13*exAt13 + 3*alpha*exAt11*exKh - - 6*alpha*Power(exAt11,2)*invgamt11 - 12*alpha*exAt11*exAt12*invgamt12 - - 12*alpha*exAt11*exAt13*invgamt13 - 6*alpha*Power(exAt12,2)*invgamt22 - - 12*alpha*exAt12*exAt13*invgamt23 - 6*alpha*Power(exAt13,2)*invgamt33 + - chi*(DDalpha11*(-3 + gam11*invgam11) + 2*DDalpha12*gam11*invgam12 + - 2*DDalpha13*gam11*invgam13 + DDalpha22*gam11*invgam22 + - 2*DDalpha23*gam11*invgam23 + DDalpha33*gam11*invgam33 + - 3*alpha*R11 - alpha*gam11*invgam11*R11 - - 2*alpha*gam11*invgam12*R12 - 2*alpha*gam11*invgam13*R13 - - alpha*gam11*invgam22*R22 - 2*alpha*gam11*invgam23*R23 - - alpha*gam11*invgam33*R33 - 24*alpha*cpi*Ss11 + - 8*alpha*cpi*gam11*invgam11*Ss11 + - 16*alpha*cpi*gam11*invgam12*Ss12 + - 16*alpha*cpi*gam11*invgam13*Ss13 + - 8*alpha*cpi*gam11*invgam22*Ss22 + - 16*alpha*cpi*gam11*invgam23*Ss23 + 8*alpha*cpi*gam11*invgam33*Ss33) \ -+ 6*alpha*exAt11*Theta)/3. -); - -local_dtexAt12.store(mask, index2, -(3*dbeta21*exAt11 + dbeta11*exAt12 + dbeta22*exAt12 - 2*dbeta33*exAt12 + - 3*dbeta23*exAt13 + 3*dbeta12*exAt22 + 3*dbeta13*exAt23 + - 3*alpha*exAt12*exKh - 6*alpha*exAt11*exAt12*invgamt11 - - 6*alpha*Power(exAt12,2)*invgamt12 - 6*alpha*exAt11*exAt22*invgamt12 - - 6*alpha*exAt12*exAt13*invgamt13 - 6*alpha*exAt11*exAt23*invgamt13 - - 6*alpha*exAt12*exAt22*invgamt22 - 6*alpha*exAt13*exAt22*invgamt23 - - 6*alpha*exAt12*exAt23*invgamt23 - 6*alpha*exAt13*exAt23*invgamt33 + - chi*(DDalpha11*gam12*invgam11 + DDalpha12*(-3 + 2*gam12*invgam12) + - 2*DDalpha13*gam12*invgam13 + DDalpha22*gam12*invgam22 + - 2*DDalpha23*gam12*invgam23 + DDalpha33*gam12*invgam33 - - alpha*gam12*invgam11*R11 + 3*alpha*R12 - - 2*alpha*gam12*invgam12*R12 - 2*alpha*gam12*invgam13*R13 - - alpha*gam12*invgam22*R22 - 2*alpha*gam12*invgam23*R23 - - alpha*gam12*invgam33*R33 + 8*alpha*cpi*gam12*invgam11*Ss11 - - 24*alpha*cpi*Ss12 + 16*alpha*cpi*gam12*invgam12*Ss12 + - 16*alpha*cpi*gam12*invgam13*Ss13 + - 8*alpha*cpi*gam12*invgam22*Ss22 + - 16*alpha*cpi*gam12*invgam23*Ss23 + 8*alpha*cpi*gam12*invgam33*Ss33) \ -+ 6*alpha*exAt12*Theta)/3. -); - -local_dtexAt13.store(mask, index2, -(3*dbeta31*exAt11 + 3*dbeta32*exAt12 + dbeta11*exAt13 - 2*dbeta22*exAt13 + - dbeta33*exAt13 + 3*dbeta12*exAt23 + 3*dbeta13*exAt33 + - 3*alpha*exAt13*exKh - 6*alpha*exAt11*exAt13*invgamt11 - - 6*alpha*exAt12*exAt13*invgamt12 - 6*alpha*exAt11*exAt23*invgamt12 - - 6*alpha*Power(exAt13,2)*invgamt13 - 6*alpha*exAt11*exAt33*invgamt13 - - 6*alpha*exAt12*exAt23*invgamt22 - 6*alpha*exAt13*exAt23*invgamt23 - - 6*alpha*exAt12*exAt33*invgamt23 - 6*alpha*exAt13*exAt33*invgamt33 + - chi*(DDalpha11*gam13*invgam11 + 2*DDalpha12*gam13*invgam12 + - DDalpha13*(-3 + 2*gam13*invgam13) + DDalpha22*gam13*invgam22 + - 2*DDalpha23*gam13*invgam23 + DDalpha33*gam13*invgam33 - - alpha*gam13*invgam11*R11 - 2*alpha*gam13*invgam12*R12 + - 3*alpha*R13 - 2*alpha*gam13*invgam13*R13 - - alpha*gam13*invgam22*R22 - 2*alpha*gam13*invgam23*R23 - - alpha*gam13*invgam33*R33 + 8*alpha*cpi*gam13*invgam11*Ss11 + - 16*alpha*cpi*gam13*invgam12*Ss12 - 24*alpha*cpi*Ss13 + - 16*alpha*cpi*gam13*invgam13*Ss13 + - 8*alpha*cpi*gam13*invgam22*Ss22 + - 16*alpha*cpi*gam13*invgam23*Ss23 + 8*alpha*cpi*gam13*invgam33*Ss33) \ -+ 6*alpha*exAt13*Theta)/3. -); - -local_dtexAt22.store(mask, index2, -(6*dbeta21*exAt12 - 2*dbeta11*exAt22 + 4*dbeta22*exAt22 - - 2*dbeta33*exAt22 + 6*dbeta23*exAt23 + 3*alpha*exAt22*exKh - - 6*alpha*Power(exAt12,2)*invgamt11 - 12*alpha*exAt12*exAt22*invgamt12 - - 12*alpha*exAt12*exAt23*invgamt13 - 6*alpha*Power(exAt22,2)*invgamt22 - - 12*alpha*exAt22*exAt23*invgamt23 - 6*alpha*Power(exAt23,2)*invgamt33 + - chi*(DDalpha11*gam22*invgam11 + 2*DDalpha12*gam22*invgam12 + - 2*DDalpha13*gam22*invgam13 + DDalpha22*(-3 + gam22*invgam22) + - 2*DDalpha23*gam22*invgam23 + DDalpha33*gam22*invgam33 - - alpha*gam22*invgam11*R11 - 2*alpha*gam22*invgam12*R12 - - 2*alpha*gam22*invgam13*R13 + 3*alpha*R22 - - alpha*gam22*invgam22*R22 - 2*alpha*gam22*invgam23*R23 - - alpha*gam22*invgam33*R33 + 8*alpha*cpi*gam22*invgam11*Ss11 + - 16*alpha*cpi*gam22*invgam12*Ss12 + - 16*alpha*cpi*gam22*invgam13*Ss13 - 24*alpha*cpi*Ss22 + - 8*alpha*cpi*gam22*invgam22*Ss22 + - 16*alpha*cpi*gam22*invgam23*Ss23 + 8*alpha*cpi*gam22*invgam33*Ss33) \ -+ 6*alpha*exAt22*Theta)/3. -); - -local_dtexAt23.store(mask, index2, -(3*dbeta31*exAt12 + 3*dbeta21*exAt13 + 3*dbeta32*exAt22 - - 2*dbeta11*exAt23 + dbeta22*exAt23 + dbeta33*exAt23 + 3*dbeta23*exAt33 + - 3*alpha*exAt23*exKh - 6*alpha*exAt12*exAt13*invgamt11 - - 6*alpha*exAt13*exAt22*invgamt12 - 6*alpha*exAt12*exAt23*invgamt12 - - 6*alpha*exAt13*exAt23*invgamt13 - 6*alpha*exAt12*exAt33*invgamt13 - - 6*alpha*exAt22*exAt23*invgamt22 - 6*alpha*Power(exAt23,2)*invgamt23 - - 6*alpha*exAt22*exAt33*invgamt23 - 6*alpha*exAt23*exAt33*invgamt33 + - chi*(DDalpha11*gam23*invgam11 + 2*DDalpha12*gam23*invgam12 + - 2*DDalpha13*gam23*invgam13 + DDalpha22*gam23*invgam22 + - DDalpha23*(-3 + 2*gam23*invgam23) + DDalpha33*gam23*invgam33 - - alpha*gam23*invgam11*R11 - 2*alpha*gam23*invgam12*R12 - - 2*alpha*gam23*invgam13*R13 - alpha*gam23*invgam22*R22 + - 3*alpha*R23 - 2*alpha*gam23*invgam23*R23 - - alpha*gam23*invgam33*R33 + 8*alpha*cpi*gam23*invgam11*Ss11 + - 16*alpha*cpi*gam23*invgam12*Ss12 + - 16*alpha*cpi*gam23*invgam13*Ss13 + - 8*alpha*cpi*gam23*invgam22*Ss22 - 24*alpha*cpi*Ss23 + - 16*alpha*cpi*gam23*invgam23*Ss23 + 8*alpha*cpi*gam23*invgam33*Ss33) \ -+ 6*alpha*exAt23*Theta)/3. -); - -local_dtexAt33.store(mask, index2, -(6*dbeta31*exAt13 + 6*dbeta32*exAt23 - 2*dbeta11*exAt33 - - 2*dbeta22*exAt33 + 4*dbeta33*exAt33 + 3*alpha*exAt33*exKh - - 6*alpha*Power(exAt13,2)*invgamt11 - 12*alpha*exAt13*exAt23*invgamt12 - - 12*alpha*exAt13*exAt33*invgamt13 - 6*alpha*Power(exAt23,2)*invgamt22 - - 12*alpha*exAt23*exAt33*invgamt23 - 6*alpha*Power(exAt33,2)*invgamt33 + - chi*(DDalpha11*gam33*invgam11 + 2*DDalpha12*gam33*invgam12 + - 2*DDalpha13*gam33*invgam13 + DDalpha22*gam33*invgam22 + - 2*DDalpha23*gam33*invgam23 + DDalpha33*(-3 + gam33*invgam33) - - alpha*gam33*invgam11*R11 - 2*alpha*gam33*invgam12*R12 - - 2*alpha*gam33*invgam13*R13 - alpha*gam33*invgam22*R22 - - 2*alpha*gam33*invgam23*R23 + 3*alpha*R33 - - alpha*gam33*invgam33*R33 + 8*alpha*cpi*gam33*invgam11*Ss11 + - 16*alpha*cpi*gam33*invgam12*Ss12 + - 16*alpha*cpi*gam33*invgam13*Ss13 + - 8*alpha*cpi*gam33*invgam22*Ss22 + - 16*alpha*cpi*gam33*invgam23*Ss23 - 24*alpha*cpi*Ss33 + - 8*alpha*cpi*gam33*invgam33*Ss33) + 6*alpha*exAt33*Theta)/3. -); - -local_dttrGt1.store(mask, index2, -(-6*dalpha1*exAtUU11 - 6*dalpha2*exAtUU12 - 6*dalpha3*exAtUU13 - - (9*alpha*(dchi1*exAtUU11 + dchi2*exAtUU12 + dchi3*exAtUU13))/chi + - 4*ddbeta111*invgamt11 + ddbeta122*invgamt11 + ddbeta133*invgamt11 + - 7*ddbeta121*invgamt12 + ddbeta222*invgamt12 + ddbeta233*invgamt12 + - 7*ddbeta131*invgamt13 + ddbeta232*invgamt13 + ddbeta333*invgamt13 + - 3*ddbeta221*invgamt22 + 6*ddbeta231*invgamt23 + 3*ddbeta331*invgamt33 - - dbeta11*trGtd1 + 2*dbeta22*trGtd1 + 2*dbeta33*trGtd1 + - 2*alpha*(3*exAtUU11*Gt111 + 6*exAtUU12*Gt112 + 6*exAtUU13*Gt113 + - 3*exAtUU22*Gt122 + 6*exAtUU23*Gt123 + 3*exAtUU33*Gt133 - - 2*dexKh1*invgamt11 - dTheta1*invgamt11 - 2*dexKh2*invgamt12 - - dTheta2*invgamt12 - 2*dexKh3*invgamt13 - dTheta3*invgamt13 - - 24*cpi*invgamt11*Sm1 - 24*cpi*invgamt12*Sm2 - - 24*cpi*invgamt13*Sm3 - 3*ckappa1*trGt1 + 3*ckappa1*trGtd1) - - 3*dbeta21*trGtd2 - 3*dbeta31*trGtd3)/3. -); - -local_dttrGt2.store(mask, index2, -(-6*dalpha1*exAtUU12 - 6*dalpha2*exAtUU22 - 6*dalpha3*exAtUU23 - - (9*alpha*(dchi1*exAtUU12 + dchi2*exAtUU22 + dchi3*exAtUU23))/chi + - 3*ddbeta112*invgamt11 + ddbeta111*invgamt12 + 7*ddbeta122*invgamt12 + - ddbeta133*invgamt12 + 6*ddbeta132*invgamt13 + ddbeta121*invgamt22 + - 4*ddbeta222*invgamt22 + ddbeta233*invgamt22 + ddbeta131*invgamt23 + - 7*ddbeta232*invgamt23 + ddbeta333*invgamt23 + 3*ddbeta332*invgamt33 - - 3*dbeta12*trGtd1 + 2*dbeta11*trGtd2 - dbeta22*trGtd2 + - 2*dbeta33*trGtd2 + 2*alpha* - (3*exAtUU11*Gt211 + 6*exAtUU12*Gt212 + 6*exAtUU13*Gt213 + - 3*exAtUU22*Gt222 + 6*exAtUU23*Gt223 + 3*exAtUU33*Gt233 - - 2*dexKh1*invgamt12 - dTheta1*invgamt12 - 2*dexKh2*invgamt22 - - dTheta2*invgamt22 - 2*dexKh3*invgamt23 - dTheta3*invgamt23 - - 24*cpi*invgamt12*Sm1 - 24*cpi*invgamt22*Sm2 - - 24*cpi*invgamt23*Sm3 - 3*ckappa1*trGt2 + 3*ckappa1*trGtd2) - - 3*dbeta32*trGtd3)/3. -); - -local_dttrGt3.store(mask, index2, -(-6*dalpha1*exAtUU13 - 6*dalpha2*exAtUU23 - 6*dalpha3*exAtUU33 - - (9*alpha*(dchi1*exAtUU13 + dchi2*exAtUU23 + dchi3*exAtUU33))/chi + - 3*ddbeta113*invgamt11 + 6*ddbeta123*invgamt12 + ddbeta111*invgamt13 + - ddbeta122*invgamt13 + 7*ddbeta133*invgamt13 + 3*ddbeta223*invgamt22 + - ddbeta121*invgamt23 + ddbeta222*invgamt23 + 7*ddbeta233*invgamt23 + - ddbeta131*invgamt33 + ddbeta232*invgamt33 + 4*ddbeta333*invgamt33 - - 3*dbeta13*trGtd1 - 3*dbeta23*trGtd2 + 2*dbeta11*trGtd3 + - 2*dbeta22*trGtd3 - dbeta33*trGtd3 + - 2*alpha*(3*exAtUU11*Gt311 + 6*exAtUU12*Gt312 + 6*exAtUU13*Gt313 + - 3*exAtUU22*Gt322 + 6*exAtUU23*Gt323 + 3*exAtUU33*Gt333 - - 2*dexKh1*invgamt13 - dTheta1*invgamt13 - 2*dexKh2*invgamt23 - - dTheta2*invgamt23 - 2*dexKh3*invgamt33 - dTheta3*invgamt33 - - 24*cpi*invgamt13*Sm1 - 24*cpi*invgamt23*Sm2 - 24*cpi*invgamt33*Sm3 - - 3*ckappa1*trGt3 + 3*ckappa1*trGtd3))/3. -); - -local_dtTheta.store(mask, index2, --0.16666666666666666*(alpha*(3*exAt11*exAtUU11 + 6*exAt12*exAtUU12 + - 6*exAt13*exAtUU13 + 3*exAt22*exAtUU22 + 6*exAt23*exAtUU23 + - 3*exAt33*exAtUU33 - 2*Power(exKh,2) + 48*cpi*rho + 12*ckappa1*Theta + - 6*ckappa1*ckappa2*Theta - 8*exKh*Theta - 8*Power(Theta,2) - 3*trR)) -); - -local_dtalpha.store(mask, index2, --(alpha*cmuL*exKh) -); - -local_dtbeta1.store(mask, index2, --(beta1*ceta) + cmuS*trGt1 -); - -local_dtbeta2.store(mask, index2, --(beta2*ceta) + cmuS*trGt2 -); - -local_dtbeta3.store(mask, index2, --(beta3*ceta) + cmuS*trGt3 -); - - - }); -}); - -#endif // #ifndef Z4COO_SET_RHS_HXX - -/* Z4coo_set_rhs.hxx */