CGH0S7
632173ea10
Previously only GAUGE 0 and 1 were supported with a compile error for 2-7. Now supports all 8 gauge choices matching BSSN Fortran formulas: - GAUGE 2: variable-eta gamma-driver, chi-sqrt denominator - GAUGE 3: variable-eta gamma-driver, chi-linear denominator - GAUGE 4: first-order variable-eta, chi-sqrt denominator - GAUGE 5: first-order variable-eta, chi-linear denominator - GAUGE 6: Jason's rational position-dependent damping - GAUGE 7: Jason's exponential position-dependent damping Also fixes dtSf advection/dissipation guards for gauges where dtSf is active. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
902 lines
42 KiB
C
902 lines
42 KiB
C
#include "macrodef.h"
|
|
#include "bssn_rhs.h"
|
|
#include "fmisc.h"
|
|
#include "ricci_gamma.h"
|
|
#include "share_func.h"
|
|
#include "tool.h"
|
|
#include <vector>
|
|
|
|
#ifdef fortran1
|
|
#define f_constraint_bssn constraint_bssn
|
|
#define f_z4c_rhs_point z4c_rhs_point
|
|
#endif
|
|
#ifdef fortran2
|
|
#define f_constraint_bssn CONSTRAINT_BSSN
|
|
#define f_z4c_rhs_point Z4C_RHS_POINT
|
|
#endif
|
|
#ifdef fortran3
|
|
#define f_constraint_bssn constraint_bssn_
|
|
#define f_z4c_rhs_point z4c_rhs_point_
|
|
#endif
|
|
|
|
extern "C" void f_constraint_bssn(int *, double *, double *, double *,
|
|
double *, double *,
|
|
double *, double *, double *, double *, double *, double *,
|
|
double *, double *, double *, double *, double *, double *,
|
|
double *, double *, double *,
|
|
double *, double *, double *, double *, double *, double *, double *, double *, double *, double *,
|
|
double *, double *, double *, double *, double *, double *,
|
|
double *, double *, double *, double *, double *, double *,
|
|
double *, double *, double *, double *, double *, double *,
|
|
double *, double *, double *, double *, double *, double *, double *, double *,
|
|
double *, double *, double *,
|
|
int &);
|
|
|
|
extern "C" void f_z4c_rhs_point(
|
|
double &A11,
|
|
double &A12,
|
|
double &A13,
|
|
double &A22,
|
|
double &A23,
|
|
double &A33,
|
|
double &alpha,
|
|
double &B1,
|
|
double &B2,
|
|
double &B3,
|
|
double &beta1,
|
|
double &beta2,
|
|
double &beta3,
|
|
double &chi,
|
|
double &chiDivFloor,
|
|
double &da1,
|
|
double &dA111,
|
|
double &dA112,
|
|
double &dA113,
|
|
double &dA122,
|
|
double &dA123,
|
|
double &dA133,
|
|
double &da2,
|
|
double &dA211,
|
|
double &dA212,
|
|
double &dA213,
|
|
double &dA222,
|
|
double &dA223,
|
|
double &dA233,
|
|
double &da3,
|
|
double &dA311,
|
|
double &dA312,
|
|
double &dA313,
|
|
double &dA322,
|
|
double &dA323,
|
|
double &dA333,
|
|
double &db11,
|
|
double &dB11,
|
|
double &db12,
|
|
double &dB12,
|
|
double &db13,
|
|
double &dB13,
|
|
double &db21,
|
|
double &dB21,
|
|
double &db22,
|
|
double &dB22,
|
|
double &db23,
|
|
double &dB23,
|
|
double &db31,
|
|
double &dB31,
|
|
double &db32,
|
|
double &dB32,
|
|
double &db33,
|
|
double &dB33,
|
|
double &dchi1,
|
|
double &dchi2,
|
|
double &dchi3,
|
|
double &dda11,
|
|
double &dda12,
|
|
double &dda13,
|
|
double &dda22,
|
|
double &dda23,
|
|
double &dda33,
|
|
double &ddb111,
|
|
double &ddb112,
|
|
double &ddb113,
|
|
double &ddb121,
|
|
double &ddb122,
|
|
double &ddb123,
|
|
double &ddb131,
|
|
double &ddb132,
|
|
double &ddb133,
|
|
double &ddb221,
|
|
double &ddb222,
|
|
double &ddb223,
|
|
double &ddb231,
|
|
double &ddb232,
|
|
double &ddb233,
|
|
double &ddb331,
|
|
double &ddb332,
|
|
double &ddb333,
|
|
double &ddchi11,
|
|
double &ddchi12,
|
|
double &ddchi13,
|
|
double &ddchi22,
|
|
double &ddchi23,
|
|
double &ddchi33,
|
|
double &deldelg1111,
|
|
double &deldelg1112,
|
|
double &deldelg1113,
|
|
double &deldelg1122,
|
|
double &deldelg1123,
|
|
double &deldelg1133,
|
|
double &deldelg1211,
|
|
double &deldelg1212,
|
|
double &deldelg1213,
|
|
double &deldelg1222,
|
|
double &deldelg1223,
|
|
double &deldelg1233,
|
|
double &deldelg1311,
|
|
double &deldelg1312,
|
|
double &deldelg1313,
|
|
double &deldelg1322,
|
|
double &deldelg1323,
|
|
double &deldelg1333,
|
|
double &deldelg2211,
|
|
double &deldelg2212,
|
|
double &deldelg2213,
|
|
double &deldelg2222,
|
|
double &deldelg2223,
|
|
double &deldelg2233,
|
|
double &deldelg2311,
|
|
double &deldelg2312,
|
|
double &deldelg2313,
|
|
double &deldelg2322,
|
|
double &deldelg2323,
|
|
double &deldelg2333,
|
|
double &deldelg3311,
|
|
double &deldelg3312,
|
|
double &deldelg3313,
|
|
double &deldelg3322,
|
|
double &deldelg3323,
|
|
double &deldelg3333,
|
|
double &delG11,
|
|
double &delg111,
|
|
double &delg112,
|
|
double &delg113,
|
|
double &delG12,
|
|
double &delg122,
|
|
double &delg123,
|
|
double &delG13,
|
|
double &delg133,
|
|
double &delG21,
|
|
double &delg211,
|
|
double &delg212,
|
|
double &delg213,
|
|
double &delG22,
|
|
double &delg222,
|
|
double &delg223,
|
|
double &delG23,
|
|
double &delg233,
|
|
double &delG31,
|
|
double &delg311,
|
|
double &delg312,
|
|
double &delg313,
|
|
double &delG32,
|
|
double &delg322,
|
|
double &delg323,
|
|
double &delG33,
|
|
double &delg333,
|
|
double &dKhat1,
|
|
double &dKhat2,
|
|
double &dKhat3,
|
|
double &dTheta1,
|
|
double &dTheta2,
|
|
double &dTheta3,
|
|
double &G1,
|
|
double &g11,
|
|
double &g12,
|
|
double &g13,
|
|
double &G2,
|
|
double &g22,
|
|
double &g23,
|
|
double &G3,
|
|
double &g33,
|
|
double &kappa1,
|
|
double &kappa2,
|
|
double &Khat,
|
|
double &rA11,
|
|
double &rA12,
|
|
double &rA13,
|
|
double &rA22,
|
|
double &rA23,
|
|
double &rA33,
|
|
double &rchi,
|
|
double &rG1,
|
|
double &rg11,
|
|
double &rg12,
|
|
double &rg13,
|
|
double &rG2,
|
|
double &rg22,
|
|
double &rg23,
|
|
double &rG3,
|
|
double &rg33,
|
|
double &rKhat,
|
|
double &rTheta,
|
|
double &Theta);
|
|
|
|
static inline void z4c_contract_gamma(
|
|
const double gxx, const double gxy, const double gxz,
|
|
const double gyy, const double gyz, const double gzz,
|
|
const double gxxx, const double gxyx, const double gxzx,
|
|
const double gyyx, const double gyzx, const double gzzx,
|
|
const double gxxy, const double gxyy, const double gxzy,
|
|
const double gyyy, const double gyzy, const double gzzy,
|
|
const double gxxz, const double gxyz, const double gxzz,
|
|
const double gyyz, const double gyzz, const double gzzz,
|
|
double &Gamxa, double &Gamya, double &Gamza)
|
|
{
|
|
double det = gxx * gyy * gzz + gxy * gyz * gxz + gxz * gxy * gyz -
|
|
gxz * gyy * gxz - gxy * gxy * gzz - gxx * gyz * gyz;
|
|
const double gupxx = (gyy * gzz - gyz * gyz) / det;
|
|
const double gupxy = -(gxy * gzz - gyz * gxz) / det;
|
|
const double gupxz = (gxy * gyz - gyy * gxz) / det;
|
|
const double gupyy = (gxx * gzz - gxz * gxz) / det;
|
|
const double gupyz = -(gxx * gyz - gxy * gxz) / det;
|
|
const double gupzz = (gxx * gyy - gxy * gxy) / det;
|
|
|
|
const double Gamxxx = 0.5 * (gupxx * gxxx + gupxy * (2.0 * gxyx - gxxy) + gupxz * (2.0 * gxzx - gxxz));
|
|
const double Gamyxx = 0.5 * (gupxy * gxxx + gupyy * (2.0 * gxyx - gxxy) + gupyz * (2.0 * gxzx - gxxz));
|
|
const double Gamzxx = 0.5 * (gupxz * gxxx + gupyz * (2.0 * gxyx - gxxy) + gupzz * (2.0 * gxzx - gxxz));
|
|
|
|
const double Gamxyy = 0.5 * (gupxx * (2.0 * gxyy - gyyx) + gupxy * gyyy + gupxz * (2.0 * gyzy - gyyz));
|
|
const double Gamyyy = 0.5 * (gupxy * (2.0 * gxyy - gyyx) + gupyy * gyyy + gupyz * (2.0 * gyzy - gyyz));
|
|
const double Gamzyy = 0.5 * (gupxz * (2.0 * gxyy - gyyx) + gupyz * gyyy + gupzz * (2.0 * gyzy - gyyz));
|
|
|
|
const double Gamxzz = 0.5 * (gupxx * (2.0 * gxzz - gzzx) + gupxy * (2.0 * gyzz - gzzy) + gupxz * gzzz);
|
|
const double Gamyzz = 0.5 * (gupxy * (2.0 * gxzz - gzzx) + gupyy * (2.0 * gyzz - gzzy) + gupyz * gzzz);
|
|
const double Gamzzz = 0.5 * (gupxz * (2.0 * gxzz - gzzx) + gupyz * (2.0 * gyzz - gzzy) + gupzz * gzzz);
|
|
|
|
const double Gamxxy = 0.5 * (gupxx * gxxy + gupxy * gyyx + gupxz * (gxzy + gyzx - gxyz));
|
|
const double Gamyxy = 0.5 * (gupxy * gxxy + gupyy * gyyx + gupyz * (gxzy + gyzx - gxyz));
|
|
const double Gamzxy = 0.5 * (gupxz * gxxy + gupyz * gyyx + gupzz * (gxzy + gyzx - gxyz));
|
|
|
|
const double Gamxxz = 0.5 * (gupxx * gxxz + gupxy * (gxyz + gyzx - gxzy) + gupxz * gzzx);
|
|
const double Gamyxz = 0.5 * (gupxy * gxxz + gupyy * (gxyz + gyzx - gxzy) + gupyz * gzzx);
|
|
const double Gamzxz = 0.5 * (gupxz * gxxz + gupyz * (gxyz + gyzx - gxzy) + gupzz * gzzx);
|
|
|
|
const double Gamxyz = 0.5 * (gupxx * (gxyz + gxzy - gyzx) + gupxy * gyyz + gupxz * gzzy);
|
|
const double Gamyyz = 0.5 * (gupxy * (gxyz + gxzy - gyzx) + gupyy * gyyz + gupyz * gzzy);
|
|
const double Gamzyz = 0.5 * (gupxz * (gxyz + gxzy - gyzx) + gupyz * gyyz + gupzz * gzzy);
|
|
|
|
Gamxa = gupxx * Gamxxx + gupyy * Gamxyy + gupzz * Gamxzz +
|
|
2.0 * (gupxy * Gamxxy + gupxz * Gamxxz + gupyz * Gamxyz);
|
|
Gamya = gupxx * Gamyxx + gupyy * Gamyyy + gupzz * Gamyzz +
|
|
2.0 * (gupxy * Gamyxy + gupxz * Gamyxz + gupyz * Gamyyz);
|
|
Gamza = gupxx * Gamzxx + gupyy * Gamzyy + gupzz * Gamzzz +
|
|
2.0 * (gupxy * Gamzxy + gupxz * Gamzxz + gupyz * Gamzyz);
|
|
}
|
|
|
|
static int compute_rhs_z4c_cartesian(
|
|
int *ex, double &T, double *X, double *Y, double *Z,
|
|
double *chi_state, double *chi_constraints, double *trK,
|
|
double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
|
|
double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
|
|
double *Gamx, double *Gamy, double *Gamz,
|
|
double *Lap, double *betax, double *betay, double *betaz,
|
|
double *dtSfx, double *dtSfy, double *dtSfz,
|
|
double *TZ,
|
|
double *chi_rhs, double *trK_rhs,
|
|
double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
|
|
double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
|
|
double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
|
|
double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
|
|
double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
|
|
double *TZ_rhs,
|
|
double *rho, double *Sx, double *Sy, double *Sz,
|
|
double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
|
|
double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
|
|
double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
|
|
double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
|
|
double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
|
|
double *Hcon, double *Mxcon, double *Mycon, double *Mzcon, double *Gmxcon, double *Gmycon, double *Gmzcon,
|
|
int &Symmetry, int &Lev, double &eps, int &co)
|
|
{
|
|
(void)T;
|
|
|
|
const int nx = ex[0];
|
|
const int ny = ex[1];
|
|
const int nz = ex[2];
|
|
const int all = nx * ny * nz;
|
|
|
|
double alpn1[all], chin1[all], gxx[all], gyy[all], gzz[all];
|
|
double chix[all], chiy[all], chiz[all], chixx[all], chixy[all], chixz[all], chiyy[all], chiyz[all], chizz[all];
|
|
double gxxx[all], gxyx[all], gxzx[all], gyyx[all], gyzx[all], gzzx[all];
|
|
double gxxy[all], gxyy[all], gxzy[all], gyyy[all], gyzy[all], gzzy[all];
|
|
double gxxz[all], gxyz[all], gxzz[all], gyyz[all], gyzz[all], gzzz[all];
|
|
double gxxxx[all], gxxxy[all], gxxxz[all], gxxyy[all], gxxyz[all], gxxzz[all];
|
|
double gxyxx[all], gxyxy[all], gxyxz[all], gxyyy[all], gxyyz[all], gxyzz[all];
|
|
double gxzxx[all], gxzxy[all], gxzxz[all], gxzyy[all], gxzyz[all], gxzzz[all];
|
|
double gyyxx[all], gyyxy[all], gyyxz[all], gyyyy[all], gyyyz[all], gyyzz[all];
|
|
double gyzxx[all], gyzxy[all], gyzxz[all], gyzyy[all], gyzyz[all], gyzzz[all];
|
|
double gzzxx[all], gzzxy[all], gzzxz[all], gzzyy[all], gzzyz[all], gzzzz[all];
|
|
double Lapx[all], Lapy[all], Lapz[all], Lapxx[all], Lapxy[all], Lapxz[all], Lapyy[all], Lapyz[all], Lapzz[all];
|
|
double betaxx[all], betaxy[all], betaxz[all], betayx[all], betayy[all], betayz[all], betazx[all], betazy[all], betazz[all];
|
|
double dBxx[all], dBxy[all], dBxz[all], dByx[all], dByy[all], dByz[all], dBzx[all], dBzy[all], dBzz[all];
|
|
double sfxxx[all], sfxxy[all], sfxxz[all], sfxyy[all], sfxyz[all], sfxzz[all];
|
|
double sfyxx[all], sfyxy[all], sfyxz[all], sfyyy[all], sfyyz[all], sfyzz[all];
|
|
double sfzxx[all], sfzxy[all], sfzxz[all], sfzyy[all], sfzyz[all], sfzzz[all];
|
|
double Gamxx[all], Gamxy[all], Gamxz[all], Gamyx[all], Gamyy[all], Gamyz[all], Gamzx[all], Gamzy[all], Gamzz[all];
|
|
double Kx[all], Ky[all], Kz[all], TZx[all], TZy[all], TZz[all];
|
|
double Axxx[all], Axxy[all], Axxz[all], Axyx[all], Axyy[all], Axyz[all];
|
|
double Axzx[all], Axzy[all], Axzz[all], Ayyx[all], Ayyy[all], Ayyz[all];
|
|
double Ayzx[all], Ayzy[all], Ayzz[all], Azzx[all], Azzy[all], Azzz[all];
|
|
#if (GAUGE == 2 || GAUGE == 3 || GAUGE == 4 || GAUGE == 5)
|
|
double reta[all];
|
|
#endif
|
|
|
|
const double SSS[3] = {1.0, 1.0, 1.0};
|
|
const double AAS[3] = {-1.0, -1.0, 1.0};
|
|
const double ASA[3] = {-1.0, 1.0, -1.0};
|
|
const double SAA[3] = {1.0, -1.0, -1.0};
|
|
const double ASS[3] = {-1.0, 1.0, 1.0};
|
|
const double SAS[3] = {1.0, -1.0, 1.0};
|
|
const double SSA[3] = {1.0, 1.0, -1.0};
|
|
|
|
const double ONE = 1.0;
|
|
const double TWO = 2.0;
|
|
const double ZEO = 0.0;
|
|
double chiDivfloor = 1.0e-5;
|
|
|
|
double kappa1 = 2.0e-2;
|
|
double kappa2 = 0.0;
|
|
double FF = 0.75;
|
|
double eta = 2.0;
|
|
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
alpn1[idx] = Lap[idx] + ONE;
|
|
chin1[idx] = chi_state[idx] + ONE;
|
|
gxx[idx] = dxx[idx] + ONE;
|
|
gyy[idx] = dyy[idx] + ONE;
|
|
gzz[idx] = dzz[idx] + ONE;
|
|
}
|
|
|
|
fderivs(ex, betax, betaxx, betaxy, betaxz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, betay, betayx, betayy, betayz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, betaz, betazx, betazy, betazz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);
|
|
fderivs(ex, dtSfx, dBxx, dBxy, dBxz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, dtSfy, dByx, dByy, dByz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, dtSfz, dBzx, dBzy, dBzz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);
|
|
fderivs(ex, chi_state, chix, chiy, chiz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, dxx, gxxx, gxxy, gxxz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, gxy, gxyx, gxyy, gxyz, X, Y, Z, -1.0, -1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, gxz, gxzx, gxzy, gxzz, X, Y, Z, -1.0, 1.0, -1.0, Symmetry, Lev);
|
|
fderivs(ex, dyy, gyyx, gyyy, gyyz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, gyz, gyzx, gyzy, gyzz, X, Y, Z, 1.0, -1.0, -1.0, Symmetry, Lev);
|
|
fderivs(ex, dzz, gzzx, gzzy, gzzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
|
|
fdderivs(ex, dxx, gxxxx, gxxxy, gxxxz, gxxyy, gxxyz, gxxzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fdderivs(ex, dyy, gyyxx, gyyxy, gyyxz, gyyyy, gyyyz, gyyzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fdderivs(ex, dzz, gzzxx, gzzxy, gzzxz, gzzyy, gzzyz, gzzzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fdderivs(ex, gxy, gxyxx, gxyxy, gxyxz, gxyyy, gxyyz, gxyzz, X, Y, Z, -1.0, -1.0, 1.0, Symmetry, Lev);
|
|
fdderivs(ex, gxz, gxzxx, gxzxy, gxzxz, gxzyy, gxzyz, gxzzz, X, Y, Z, -1.0, 1.0, -1.0, Symmetry, Lev);
|
|
fdderivs(ex, gyz, gyzxx, gyzxy, gyzxz, gyzyy, gyzyz, gyzzz, X, Y, Z, 1.0, -1.0, -1.0, Symmetry, Lev);
|
|
|
|
fderivs(ex, Gamx, Gamxx, Gamxy, Gamxz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, Gamy, Gamyx, Gamyy, Gamyz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, Gamz, Gamzx, Gamzy, Gamzz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);
|
|
|
|
fderivs(ex, Lap, Lapx, Lapy, Lapz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, trK, Kx, Ky, Kz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, TZ, TZx, TZy, TZz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
|
|
fdderivs(ex, betax, sfxxx, sfxxy, sfxxz, sfxyy, sfxyz, sfxzz, X, Y, Z, -1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fdderivs(ex, betay, sfyxx, sfyxy, sfyxz, sfyyy, sfyyz, sfyzz, X, Y, Z, 1.0, -1.0, 1.0, Symmetry, Lev);
|
|
fdderivs(ex, betaz, sfzxx, sfzxy, sfzxz, sfzyy, sfzyz, sfzzz, X, Y, Z, 1.0, 1.0, -1.0, Symmetry, Lev);
|
|
|
|
fdderivs(ex, chi_state, chixx, chixy, chixz, chiyy, chiyz, chizz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fdderivs(ex, Lap, Lapxx, Lapxy, Lapxz, Lapyy, Lapyz, Lapzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
|
|
fderivs(ex, Axx, Axxx, Axxy, Axxz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, Axy, Axyx, Axyy, Axyz, X, Y, Z, -1.0, -1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, Axz, Axzx, Axzy, Axzz, X, Y, Z, -1.0, 1.0, -1.0, Symmetry, Lev);
|
|
fderivs(ex, Ayy, Ayyx, Ayyy, Ayyz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
fderivs(ex, Ayz, Ayzx, Ayzy, Ayzz, X, Y, Z, 1.0, -1.0, -1.0, Symmetry, Lev);
|
|
fderivs(ex, Azz, Azzx, Azzy, Azzz, X, Y, Z, 1.0, 1.0, 1.0, Symmetry, Lev);
|
|
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
double point_kappa1 = 0.0;
|
|
f_z4c_rhs_point(
|
|
Axx[idx], Axy[idx], Axz[idx], Ayy[idx], Ayz[idx], Azz[idx],
|
|
alpn1[idx], dtSfx[idx], dtSfy[idx], dtSfz[idx],
|
|
betax[idx], betay[idx], betaz[idx],
|
|
chin1[idx], chiDivfloor,
|
|
Lapx[idx],
|
|
Axxx[idx], Axyx[idx], Axzx[idx], Ayyx[idx], Ayzx[idx], Azzx[idx],
|
|
Lapy[idx],
|
|
Axxy[idx], Axyy[idx], Axzy[idx], Ayyy[idx], Ayzy[idx], Azzy[idx],
|
|
Lapz[idx],
|
|
Axxz[idx], Axyz[idx], Axzz[idx], Ayyz[idx], Ayzz[idx], Azzz[idx],
|
|
betaxx[idx], dBxx[idx], betayx[idx], dByx[idx], betazx[idx], dBzx[idx],
|
|
betaxy[idx], dBxy[idx], betayy[idx], dByy[idx], betazy[idx], dBzy[idx],
|
|
betaxz[idx], dBxz[idx], betayz[idx], dByz[idx], betazz[idx], dBzz[idx],
|
|
chix[idx], chiy[idx], chiz[idx],
|
|
Lapxx[idx], Lapxy[idx], Lapxz[idx], Lapyy[idx], Lapyz[idx], Lapzz[idx],
|
|
sfxxx[idx], sfyxx[idx], sfzxx[idx],
|
|
sfxxy[idx], sfyxy[idx], sfzxy[idx],
|
|
sfxxz[idx], sfyxz[idx], sfzxz[idx],
|
|
sfxyy[idx], sfyyy[idx], sfzyy[idx],
|
|
sfxyz[idx], sfyyz[idx], sfzyz[idx],
|
|
sfxzz[idx], sfyzz[idx], sfzzz[idx],
|
|
chixx[idx], chixy[idx], chixz[idx], chiyy[idx], chiyz[idx], chizz[idx],
|
|
gxxxx[idx], gxyxx[idx], gxzxx[idx], gyyxx[idx], gyzxx[idx], gzzxx[idx],
|
|
gxxxy[idx], gxyxy[idx], gxzxy[idx], gyyxy[idx], gyzxy[idx], gzzxy[idx],
|
|
gxxxz[idx], gxyxz[idx], gxzxz[idx], gyyxz[idx], gyzxz[idx], gzzxz[idx],
|
|
gxxyy[idx], gxyyy[idx], gxzyy[idx], gyyyy[idx], gyzyy[idx], gzzyy[idx],
|
|
gxxyz[idx], gxyyz[idx], gxzyz[idx], gyyyz[idx], gyzyz[idx], gzzyz[idx],
|
|
gxxzz[idx], gxyzz[idx], gxzzz[idx], gyyzz[idx], gyzzz[idx], gzzzz[idx],
|
|
Gamxx[idx], gxxx[idx], gxyx[idx], gxzx[idx],
|
|
Gamyx[idx], gyyx[idx], gyzx[idx],
|
|
Gamzx[idx], gzzx[idx],
|
|
Gamxy[idx], gxxy[idx], gxyy[idx], gxzy[idx],
|
|
Gamyy[idx], gyyy[idx], gyzy[idx],
|
|
Gamzy[idx], gzzy[idx],
|
|
Gamxz[idx], gxxz[idx], gxyz[idx], gxzz[idx],
|
|
Gamyz[idx], gyyz[idx], gyzz[idx],
|
|
Gamzz[idx], gzzz[idx],
|
|
Kx[idx], Ky[idx], Kz[idx],
|
|
TZx[idx], TZy[idx], TZz[idx],
|
|
Gamx[idx], gxx[idx], gxy[idx], gxz[idx],
|
|
Gamy[idx], gyy[idx], gyz[idx],
|
|
Gamz[idx], gzz[idx],
|
|
point_kappa1, kappa2,
|
|
trK[idx],
|
|
Axx_rhs[idx], Axy_rhs[idx], Axz_rhs[idx], Ayy_rhs[idx], Ayz_rhs[idx], Azz_rhs[idx],
|
|
chi_rhs[idx],
|
|
Gamx_rhs[idx], gxx_rhs[idx], gxy_rhs[idx], gxz_rhs[idx],
|
|
Gamy_rhs[idx], gyy_rhs[idx], gyz_rhs[idx],
|
|
Gamz_rhs[idx], gzz_rhs[idx], trK_rhs[idx], TZ_rhs[idx], TZ[idx]);
|
|
}
|
|
|
|
for (int idx = 0; idx < all; ++idx)
|
|
Lap_rhs[idx] = -TWO * alpn1[idx] * trK[idx];
|
|
|
|
#if (GAUGE == 0)
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
betax_rhs[idx] = FF * dtSfx[idx];
|
|
betay_rhs[idx] = FF * dtSfy[idx];
|
|
betaz_rhs[idx] = FF * dtSfz[idx];
|
|
dtSfx_rhs[idx] = Gamx_rhs[idx] - eta * dtSfx[idx];
|
|
dtSfy_rhs[idx] = Gamy_rhs[idx] - eta * dtSfy[idx];
|
|
dtSfz_rhs[idx] = Gamz_rhs[idx] - eta * dtSfz[idx];
|
|
}
|
|
#elif (GAUGE == 1)
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
betax_rhs[idx] = Gamx[idx] - eta * betax[idx];
|
|
betay_rhs[idx] = Gamy[idx] - eta * betay[idx];
|
|
betaz_rhs[idx] = Gamz[idx] - eta * betaz[idx];
|
|
dtSfx_rhs[idx] = ZEO;
|
|
dtSfy_rhs[idx] = ZEO;
|
|
dtSfz_rhs[idx] = ZEO;
|
|
}
|
|
#elif (GAUGE == 2)
|
|
/* Variable-eta gamma-driver, chi-sqrt denominator */
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
const double chin1i = chin1[idx];
|
|
const double det = gxx[idx] * gyy[idx] * gzz[idx]
|
|
+ gxy[idx] * gyz[idx] * gxz[idx] * 2.0
|
|
- gxz[idx] * gyy[idx] * gxz[idx]
|
|
- gxy[idx] * gxy[idx] * gzz[idx]
|
|
- gxx[idx] * gyz[idx] * gyz[idx];
|
|
const double idet = ONE / det;
|
|
const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
|
|
const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
|
|
const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
|
|
const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
|
|
const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
|
|
const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
|
|
const double grdchi2 =
|
|
upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
|
|
+ TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
|
|
const double sqchi = sqrt(chin1i);
|
|
reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - sqchi) * (ONE - sqchi));
|
|
betax_rhs[idx] = FF * dtSfx[idx];
|
|
betay_rhs[idx] = FF * dtSfy[idx];
|
|
betaz_rhs[idx] = FF * dtSfz[idx];
|
|
dtSfx_rhs[idx] = Gamx_rhs[idx] - reta[idx] * dtSfx[idx];
|
|
dtSfy_rhs[idx] = Gamy_rhs[idx] - reta[idx] * dtSfy[idx];
|
|
dtSfz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * dtSfz[idx];
|
|
}
|
|
#elif (GAUGE == 3)
|
|
/* Variable-eta gamma-driver, chi-linear denominator */
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
const double chin1i = chin1[idx];
|
|
const double det = gxx[idx] * gyy[idx] * gzz[idx]
|
|
+ gxy[idx] * gyz[idx] * gxz[idx] * 2.0
|
|
- gxz[idx] * gyy[idx] * gxz[idx]
|
|
- gxy[idx] * gxy[idx] * gzz[idx]
|
|
- gxx[idx] * gyz[idx] * gyz[idx];
|
|
const double idet = ONE / det;
|
|
const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
|
|
const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
|
|
const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
|
|
const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
|
|
const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
|
|
const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
|
|
const double grdchi2 =
|
|
upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
|
|
+ TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
|
|
reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - chin1i) * (ONE - chin1i));
|
|
betax_rhs[idx] = FF * dtSfx[idx];
|
|
betay_rhs[idx] = FF * dtSfy[idx];
|
|
betaz_rhs[idx] = FF * dtSfz[idx];
|
|
dtSfx_rhs[idx] = Gamx_rhs[idx] - reta[idx] * dtSfx[idx];
|
|
dtSfy_rhs[idx] = Gamy_rhs[idx] - reta[idx] * dtSfy[idx];
|
|
dtSfz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * dtSfz[idx];
|
|
}
|
|
#elif (GAUGE == 4)
|
|
/* Variable-eta gamma-driver, first-order, chi-sqrt denominator */
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
const double chin1i = chin1[idx];
|
|
const double det = gxx[idx] * gyy[idx] * gzz[idx]
|
|
+ gxy[idx] * gyz[idx] * gxz[idx] * 2.0
|
|
- gxz[idx] * gyy[idx] * gxz[idx]
|
|
- gxy[idx] * gxy[idx] * gzz[idx]
|
|
- gxx[idx] * gyz[idx] * gyz[idx];
|
|
const double idet = ONE / det;
|
|
const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
|
|
const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
|
|
const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
|
|
const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
|
|
const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
|
|
const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
|
|
const double grdchi2 =
|
|
upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
|
|
+ TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
|
|
const double sqchi = sqrt(chin1i);
|
|
reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - sqchi) * (ONE - sqchi));
|
|
betax_rhs[idx] = Gamx_rhs[idx] - reta[idx] * betax[idx];
|
|
betay_rhs[idx] = Gamy_rhs[idx] - reta[idx] * betay[idx];
|
|
betaz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * betaz[idx];
|
|
dtSfx_rhs[idx] = ZEO;
|
|
dtSfy_rhs[idx] = ZEO;
|
|
dtSfz_rhs[idx] = ZEO;
|
|
}
|
|
#elif (GAUGE == 5)
|
|
/* Variable-eta gamma-driver, first-order, chi-linear denominator */
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
const double chin1i = chin1[idx];
|
|
const double det = gxx[idx] * gyy[idx] * gzz[idx]
|
|
+ gxy[idx] * gyz[idx] * gxz[idx] * 2.0
|
|
- gxz[idx] * gyy[idx] * gxz[idx]
|
|
- gxy[idx] * gxy[idx] * gzz[idx]
|
|
- gxx[idx] * gyz[idx] * gyz[idx];
|
|
const double idet = ONE / det;
|
|
const double upxx = (gyy[idx] * gzz[idx] - gyz[idx] * gyz[idx]) * idet;
|
|
const double upxy = -(gxy[idx] * gzz[idx] - gyz[idx] * gxz[idx]) * idet;
|
|
const double upxz = (gxy[idx] * gyz[idx] - gyy[idx] * gxz[idx]) * idet;
|
|
const double upyy = (gxx[idx] * gzz[idx] - gxz[idx] * gxz[idx]) * idet;
|
|
const double upyz = -(gxx[idx] * gyz[idx] - gxy[idx] * gxz[idx]) * idet;
|
|
const double upzz = (gxx[idx] * gyy[idx] - gxy[idx] * gxy[idx]) * idet;
|
|
const double grdchi2 =
|
|
upxx * chix[idx] * chix[idx] + upyy * chiy[idx] * chiy[idx] + upzz * chiz[idx] * chiz[idx]
|
|
+ TWO * (upxy * chix[idx] * chiy[idx] + upxz * chix[idx] * chiz[idx] + upyz * chiy[idx] * chiz[idx]);
|
|
reta[idx] = 1.31 / TWO * sqrt(grdchi2 / chin1i) / ((ONE - chin1i) * (ONE - chin1i));
|
|
betax_rhs[idx] = Gamx_rhs[idx] - reta[idx] * betax[idx];
|
|
betay_rhs[idx] = Gamy_rhs[idx] - reta[idx] * betay[idx];
|
|
betaz_rhs[idx] = Gamz_rhs[idx] - reta[idx] * betaz[idx];
|
|
dtSfx_rhs[idx] = ZEO;
|
|
dtSfy_rhs[idx] = ZEO;
|
|
dtSfz_rhs[idx] = ZEO;
|
|
}
|
|
#elif (GAUGE == 6 || GAUGE == 7)
|
|
{
|
|
/* Jason's position-dependent damping: rational (6) or exponential (7) */
|
|
int BHN = 0;
|
|
double Porg[9] = {0.0};
|
|
double Mass[3] = {0.0};
|
|
#ifdef fortran1
|
|
extern "C" { void getpbh(int &, double *, double *); }
|
|
#elif defined(fortran2)
|
|
extern "C" { void GETPBH(int &, double *, double *); }
|
|
#else
|
|
extern "C" { void getpbh_(int &, double *, double *); }
|
|
#endif
|
|
{
|
|
#ifdef fortran1
|
|
getpbh(BHN, Porg, Mass);
|
|
#elif defined(fortran2)
|
|
GETPBH(BHN, Porg, Mass);
|
|
#else
|
|
getpbh_(BHN, Porg, Mass);
|
|
#endif
|
|
}
|
|
if (BHN == 2)
|
|
{
|
|
const double M = Mass[0] + Mass[1];
|
|
const double A = 2.0 / M;
|
|
const double w1 = 12.0, w2 = 12.0;
|
|
const double C1 = 1.0 / Mass[0] - A;
|
|
const double C2 = 1.0 / Mass[1] - A;
|
|
const double BH_sep2 = (Porg[3] - Porg[0]) * (Porg[3] - Porg[0])
|
|
+ (Porg[4] - Porg[1]) * (Porg[4] - Porg[1])
|
|
+ (Porg[5] - Porg[2]) * (Porg[5] - Porg[2]);
|
|
const double inv_BH_sep2 = 1.0 / BH_sep2;
|
|
for (int k0 = 0; k0 < nz; ++k0) {
|
|
for (int j0 = 0; j0 < ny; ++j0) {
|
|
for (int i0 = 0; i0 < nx; ++i0) {
|
|
const size_t idx = idx_ex(i0, j0, k0, ex);
|
|
const double xp = X[i0], yp = Y[j0], zp = Z[k0];
|
|
const double r1 = ((Porg[0]-xp)*(Porg[0]-xp) + (Porg[1]-yp)*(Porg[1]-yp) + (Porg[2]-zp)*(Porg[2]-zp)) * inv_BH_sep2;
|
|
const double r2 = ((Porg[3]-xp)*(Porg[3]-xp) + (Porg[4]-yp)*(Porg[4]-yp) + (Porg[5]-zp)*(Porg[5]-zp)) * inv_BH_sep2;
|
|
#if (GAUGE == 6)
|
|
const double reta_val = A + C1 / (1.0 + w1 * r1) + C2 / (1.0 + w2 * r2);
|
|
#else
|
|
const double reta_val = A + C1 * exp(-w1 * r1) + C2 * exp(-w2 * r2);
|
|
#endif
|
|
betax_rhs[idx] = FF * dtSfx[idx];
|
|
betay_rhs[idx] = FF * dtSfy[idx];
|
|
betaz_rhs[idx] = FF * dtSfz[idx];
|
|
dtSfx_rhs[idx] = Gamx_rhs[idx] - reta_val * dtSfx[idx];
|
|
dtSfy_rhs[idx] = Gamy_rhs[idx] - reta_val * dtSfy[idx];
|
|
dtSfz_rhs[idx] = Gamz_rhs[idx] - reta_val * dtSfz[idx];
|
|
}}}
|
|
}
|
|
else
|
|
{
|
|
fprintf(stderr, "z4c_rhs_c: GAUGE %d requires BHN=2, got BHN=%d\n", (int)GAUGE, BHN);
|
|
return 1;
|
|
}
|
|
}
|
|
#else
|
|
#error "z4c_rhs_c.C: unsupported GAUGE value"
|
|
#endif
|
|
|
|
lopsided(ex, X, Y, Z, gxx, gxx_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
lopsided(ex, X, Y, Z, gxy, gxy_rhs, betax, betay, betaz, Symmetry, AAS);
|
|
lopsided(ex, X, Y, Z, gxz, gxz_rhs, betax, betay, betaz, Symmetry, ASA);
|
|
lopsided(ex, X, Y, Z, gyy, gyy_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
lopsided(ex, X, Y, Z, gyz, gyz_rhs, betax, betay, betaz, Symmetry, SAA);
|
|
lopsided(ex, X, Y, Z, gzz, gzz_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
|
|
lopsided(ex, X, Y, Z, Axx, Axx_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
lopsided(ex, X, Y, Z, Axy, Axy_rhs, betax, betay, betaz, Symmetry, AAS);
|
|
lopsided(ex, X, Y, Z, Axz, Axz_rhs, betax, betay, betaz, Symmetry, ASA);
|
|
lopsided(ex, X, Y, Z, Ayy, Ayy_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
lopsided(ex, X, Y, Z, Ayz, Ayz_rhs, betax, betay, betaz, Symmetry, SAA);
|
|
lopsided(ex, X, Y, Z, Azz, Azz_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
|
|
lopsided(ex, X, Y, Z, chi_state, chi_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
lopsided(ex, X, Y, Z, trK, trK_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
|
|
lopsided(ex, X, Y, Z, Gamx, Gamx_rhs, betax, betay, betaz, Symmetry, ASS);
|
|
lopsided(ex, X, Y, Z, Gamy, Gamy_rhs, betax, betay, betaz, Symmetry, SAS);
|
|
lopsided(ex, X, Y, Z, Gamz, Gamz_rhs, betax, betay, betaz, Symmetry, SSA);
|
|
|
|
lopsided(ex, X, Y, Z, Lap, Lap_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
lopsided(ex, X, Y, Z, betax, betax_rhs, betax, betay, betaz, Symmetry, ASS);
|
|
lopsided(ex, X, Y, Z, betay, betay_rhs, betax, betay, betaz, Symmetry, SAS);
|
|
lopsided(ex, X, Y, Z, betaz, betaz_rhs, betax, betay, betaz, Symmetry, SSA);
|
|
#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
|
|
lopsided(ex, X, Y, Z, dtSfx, dtSfx_rhs, betax, betay, betaz, Symmetry, ASS);
|
|
lopsided(ex, X, Y, Z, dtSfy, dtSfy_rhs, betax, betay, betaz, Symmetry, SAS);
|
|
lopsided(ex, X, Y, Z, dtSfz, dtSfz_rhs, betax, betay, betaz, Symmetry, SSA);
|
|
#endif
|
|
lopsided(ex, X, Y, Z, TZ, TZ_rhs, betax, betay, betaz, Symmetry, SSS);
|
|
|
|
for (int idx = 0; idx < all; ++idx)
|
|
{
|
|
double Gamxa = 0.0, Gamya = 0.0, Gamza = 0.0;
|
|
z4c_contract_gamma(
|
|
gxx[idx], gxy[idx], gxz[idx], gyy[idx], gyz[idx], gzz[idx],
|
|
gxxx[idx], gxyx[idx], gxzx[idx], gyyx[idx], gyzx[idx], gzzx[idx],
|
|
gxxy[idx], gxyy[idx], gxzy[idx], gyyy[idx], gyzy[idx], gzzy[idx],
|
|
gxxz[idx], gxyz[idx], gxzz[idx], gyyz[idx], gyzz[idx], gzzz[idx],
|
|
Gamxa, Gamya, Gamza);
|
|
|
|
TZ_rhs[idx] -= alpn1[idx] * (TWO + kappa2) * kappa1 * TZ[idx];
|
|
trK_rhs[idx] += alpn1[idx] * kappa1 * (ONE - kappa2) * TZ[idx];
|
|
Gamx_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamx[idx] - Gamxa);
|
|
Gamy_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamy[idx] - Gamya);
|
|
Gamz_rhs[idx] -= TWO * alpn1[idx] * kappa1 * (Gamz[idx] - Gamza);
|
|
}
|
|
|
|
if (eps > 0.0)
|
|
{
|
|
kodis(ex, X, Y, Z, chi_state, chi_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, trK, trK_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, gxx, gxx_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, gxy, gxy_rhs, AAS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, gxz, gxz_rhs, ASA, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, gyy, gyy_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, gyz, gyz_rhs, SAA, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, gzz, gzz_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Axx, Axx_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Axy, Axy_rhs, AAS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Axz, Axz_rhs, ASA, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Ayy, Ayy_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Ayz, Ayz_rhs, SAA, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Azz, Azz_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Gamx, Gamx_rhs, ASS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Gamy, Gamy_rhs, SAS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Gamz, Gamz_rhs, SSA, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, Lap, Lap_rhs, SSS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, betax, betax_rhs, ASS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, betay, betay_rhs, SAS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, betaz, betaz_rhs, SSA, Symmetry, eps);
|
|
#if (GAUGE == 0 || GAUGE == 2 || GAUGE == 3 || GAUGE == 6 || GAUGE == 7)
|
|
kodis(ex, X, Y, Z, dtSfx, dtSfx_rhs, ASS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, dtSfy, dtSfy_rhs, SAS, Symmetry, eps);
|
|
kodis(ex, X, Y, Z, dtSfz, dtSfz_rhs, SSA, Symmetry, eps);
|
|
#endif
|
|
kodis(ex, X, Y, Z, TZ, TZ_rhs, SSS, Symmetry, eps);
|
|
}
|
|
|
|
if (co == 0)
|
|
{
|
|
#if (ABV == 0)
|
|
f_ricci_gamma(ex, X, Y, Z,
|
|
chi_constraints,
|
|
dxx, gxy, gxz, dyy, gyz, dzz,
|
|
Gamx, Gamy, Gamz,
|
|
Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
|
|
Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
|
|
Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
|
|
Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
|
|
Symmetry);
|
|
#endif
|
|
f_constraint_bssn(ex, X, Y, Z,
|
|
chi_constraints, trK,
|
|
dxx, gxy, gxz, dyy, gyz, dzz,
|
|
Axx, Axy, Axz, Ayy, Ayz, Azz,
|
|
Gamx, Gamy, Gamz,
|
|
Lap, betax, betay, betaz, rho, Sx, Sy, Sz,
|
|
Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
|
|
Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
|
|
Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
|
|
Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
|
|
Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
|
|
Symmetry);
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
extern "C" int f_compute_rhs_Z4c(int *ex, double &T,
|
|
double *X, double *Y, double *Z,
|
|
double *chi, double *trK,
|
|
double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
|
|
double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
|
|
double *Gamx, double *Gamy, double *Gamz,
|
|
double *Lap, double *betax, double *betay, double *betaz,
|
|
double *dtSfx, double *dtSfy, double *dtSfz,
|
|
double *TZ,
|
|
double *chi_rhs, double *trK_rhs,
|
|
double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
|
|
double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
|
|
double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
|
|
double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
|
|
double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
|
|
double *TZ_rhs,
|
|
double *rho, double *Sx, double *Sy, double *Sz,
|
|
double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
|
|
double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
|
|
double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
|
|
double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
|
|
double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
|
|
double *Hcon, double *Mxcon, double *Mycon, double *Mzcon, double *Gmxcon, double *Gmycon, double *Gmzcon,
|
|
int &Symmetry, int &Lev, double &eps, int &co)
|
|
{
|
|
return compute_rhs_z4c_cartesian(
|
|
ex, T, X, Y, Z,
|
|
chi, chi, trK,
|
|
dxx, gxy, gxz, dyy, gyz, dzz,
|
|
Axx, Axy, Axz, Ayy, Ayz, Azz,
|
|
Gamx, Gamy, Gamz,
|
|
Lap, betax, betay, betaz,
|
|
dtSfx, dtSfy, dtSfz,
|
|
TZ,
|
|
chi_rhs, trK_rhs,
|
|
gxx_rhs, gxy_rhs, gxz_rhs, gyy_rhs, gyz_rhs, gzz_rhs,
|
|
Axx_rhs, Axy_rhs, Axz_rhs, Ayy_rhs, Ayz_rhs, Azz_rhs,
|
|
Gamx_rhs, Gamy_rhs, Gamz_rhs,
|
|
Lap_rhs, betax_rhs, betay_rhs, betaz_rhs,
|
|
dtSfx_rhs, dtSfy_rhs, dtSfz_rhs,
|
|
TZ_rhs,
|
|
rho, Sx, Sy, Sz,
|
|
Sxx, Sxy, Sxz, Syy, Syz, Szz,
|
|
Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
|
|
Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
|
|
Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
|
|
Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
|
|
Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
|
|
Symmetry, Lev, eps, co);
|
|
}
|
|
|
|
extern "C" int f_compute_rhs_Z4cnot(int *ex, double &T,
|
|
double *X, double *Y, double *Z,
|
|
double *chi, double *trK,
|
|
double *dxx, double *gxy, double *gxz, double *dyy, double *gyz, double *dzz,
|
|
double *Axx, double *Axy, double *Axz, double *Ayy, double *Ayz, double *Azz,
|
|
double *Gamx, double *Gamy, double *Gamz,
|
|
double *Lap, double *betax, double *betay, double *betaz,
|
|
double *dtSfx, double *dtSfy, double *dtSfz,
|
|
double *TZ,
|
|
double *chi_rhs, double *trK_rhs,
|
|
double *gxx_rhs, double *gxy_rhs, double *gxz_rhs, double *gyy_rhs, double *gyz_rhs, double *gzz_rhs,
|
|
double *Axx_rhs, double *Axy_rhs, double *Axz_rhs, double *Ayy_rhs, double *Ayz_rhs, double *Azz_rhs,
|
|
double *Gamx_rhs, double *Gamy_rhs, double *Gamz_rhs,
|
|
double *Lap_rhs, double *betax_rhs, double *betay_rhs, double *betaz_rhs,
|
|
double *dtSfx_rhs, double *dtSfy_rhs, double *dtSfz_rhs,
|
|
double *TZ_rhs,
|
|
double *rho, double *Sx, double *Sy, double *Sz,
|
|
double *Sxx, double *Sxy, double *Sxz, double *Syy, double *Syz, double *Szz,
|
|
double *Gamxxx, double *Gamxxy, double *Gamxxz, double *Gamxyy, double *Gamxyz, double *Gamxzz,
|
|
double *Gamyxx, double *Gamyxy, double *Gamyxz, double *Gamyyy, double *Gamyyz, double *Gamyzz,
|
|
double *Gamzxx, double *Gamzxy, double *Gamzxz, double *Gamzyy, double *Gamzyz, double *Gamzzz,
|
|
double *Rxx, double *Rxy, double *Rxz, double *Ryy, double *Ryz, double *Rzz,
|
|
double *Hcon, double *Mxcon, double *Mycon, double *Mzcon, double *Gmxcon, double *Gmycon, double *Gmzcon,
|
|
int &Symmetry, int &Lev, double &eps, int &co, double &chitiny)
|
|
{
|
|
const int all = ex[0] * ex[1] * ex[2];
|
|
std::vector<double> chi_clamped(chi, chi + all);
|
|
f_lowerboundset(ex, chi_clamped.data(), chitiny);
|
|
|
|
const int ret = compute_rhs_z4c_cartesian(
|
|
ex, T, X, Y, Z,
|
|
chi_clamped.data(), chi, trK,
|
|
dxx, gxy, gxz, dyy, gyz, dzz,
|
|
Axx, Axy, Axz, Ayy, Ayz, Azz,
|
|
Gamx, Gamy, Gamz,
|
|
Lap, betax, betay, betaz,
|
|
dtSfx, dtSfy, dtSfz,
|
|
TZ,
|
|
chi_rhs, trK_rhs,
|
|
gxx_rhs, gxy_rhs, gxz_rhs, gyy_rhs, gyz_rhs, gzz_rhs,
|
|
Axx_rhs, Axy_rhs, Axz_rhs, Ayy_rhs, Ayz_rhs, Azz_rhs,
|
|
Gamx_rhs, Gamy_rhs, Gamz_rhs,
|
|
Lap_rhs, betax_rhs, betay_rhs, betaz_rhs,
|
|
dtSfx_rhs, dtSfy_rhs, dtSfz_rhs,
|
|
TZ_rhs,
|
|
rho, Sx, Sy, Sz,
|
|
Sxx, Sxy, Sxz, Syy, Syz, Szz,
|
|
Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
|
|
Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
|
|
Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
|
|
Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
|
|
Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
|
|
Symmetry, Lev, eps, co);
|
|
|
|
if (ret != 0 || co != 0)
|
|
return ret;
|
|
|
|
#if (ABV == 0)
|
|
f_ricci_gamma(ex, X, Y, Z,
|
|
chi,
|
|
dxx, gxy, gxz, dyy, gyz, dzz,
|
|
Gamx, Gamy, Gamz,
|
|
Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
|
|
Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
|
|
Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
|
|
Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
|
|
Symmetry);
|
|
#endif
|
|
f_constraint_bssn(ex, X, Y, Z,
|
|
chi, trK,
|
|
dxx, gxy, gxz, dyy, gyz, dzz,
|
|
Axx, Axy, Axz, Ayy, Ayz, Azz,
|
|
Gamx, Gamy, Gamz,
|
|
Lap, betax, betay, betaz, rho, Sx, Sy, Sz,
|
|
Gamxxx, Gamxxy, Gamxxz, Gamxyy, Gamxyz, Gamxzz,
|
|
Gamyxx, Gamyxy, Gamyxz, Gamyyy, Gamyyz, Gamyzz,
|
|
Gamzxx, Gamzxy, Gamzxz, Gamzyy, Gamzyz, Gamzzz,
|
|
Rxx, Rxy, Rxz, Ryy, Ryz, Rzz,
|
|
Hcon, Mxcon, Mycon, Mzcon, Gmxcon, Gmycon, Gmzcon,
|
|
Symmetry);
|
|
return ret;
|
|
}
|