#include "macrodef.h" #include "bssn_rhs.h" #include "fmisc.h" #include "ricci_gamma.h" #include "share_func.h" #include "tool.h" #include #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]; 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; } #else #error "z4c_rhs_c.C currently supports GAUGE == 0 or GAUGE == 1 for Z4C" #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) 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) 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 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; }