diff --git a/AMSS_NCKU_source/Z4c_rhs.f90 b/AMSS_NCKU_source/Z4c_rhs.f90 index 3b877ea..3209c4d 100644 --- a/AMSS_NCKU_source/Z4c_rhs.f90 +++ b/AMSS_NCKU_source/Z4c_rhs.f90 @@ -94,29 +94,31 @@ Hcon,Mxcon,Mycon,Mzcon,Gmxcon,Gmycon,Gmzcon, & Symmetry,Lev,eps,co) -#if (ABV == 0) - call 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 - call 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) + if (co == 0) then +#if (ABV == 0) + call 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 + call 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) + endif return @@ -226,11 +228,12 @@ call get_Z4cparameters(kappa1,kappa2,kappa3,FF,eta) -!!! sanity check - dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) & - +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) & - +sum(Gamx)+sum(Gamy)+sum(Gamz) & - +sum(Lap)+sum(betax)+sum(betay)+sum(betaz)+sum(dtSfx)+sum(dtSfy)+sum(dtSfz) & +!!! sanity check +#ifdef DEBUG + dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) & + +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) & + +sum(Gamx)+sum(Gamy)+sum(Gamz) & + +sum(Lap)+sum(betax)+sum(betay)+sum(betaz)+sum(dtSfx)+sum(dtSfy)+sum(dtSfz) & +sum(TZ) if(dX.ne.dX) then if(sum(chi).ne.sum(chi))write(*,*)"Z4c_rhs.f90: find NaN in chi" @@ -257,10 +260,11 @@ if(sum(dtSfx).ne.sum(dtSfx))write(*,*)"Z4c_rhs.f90: find NaN in dtSfx" if(sum(dtSfy).ne.sum(dtSfy))write(*,*)"Z4c_rhs.f90: find NaN in dtSfy" if(sum(dtSfz).ne.sum(dtSfz))write(*,*)"Z4c_rhs.f90: find NaN in dtSfz" - if(sum(TZ).ne.sum(Tz))write(*,*)"Z4c_rhs.f90: find NaN in TZ" - gont = 1 - return - endif + if(sum(TZ).ne.sum(Tz))write(*,*)"Z4c_rhs.f90: find NaN in TZ" + gont = 1 + return + endif +#endif PI = dacos(-ONE) @@ -1263,30 +1267,32 @@ endif -#if (ABV == 0) - call 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 - - call 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) + if (co == 0) then +#if (ABV == 0) + call 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 + + call 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) + endif gont = 0 diff --git a/AMSS_NCKU_source/Z4c_rhs_ss.f90 b/AMSS_NCKU_source/Z4c_rhs_ss.f90 index 173a1e9..115cd11 100644 --- a/AMSS_NCKU_source/Z4c_rhs_ss.f90 +++ b/AMSS_NCKU_source/Z4c_rhs_ss.f90 @@ -121,11 +121,12 @@ call get_Z4cparameters(kappa1,kappa2,kappa3,FF,eta) -!!! sanity check - dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) & - +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) & - +sum(Gamx)+sum(Gamy)+sum(Gamz) & - +sum(Lap)+sum(betax)+sum(betay)+sum(betaz)+sum(dtSfx)+sum(dtSfy)+sum(dtSfz) & +!!! sanity check +#ifdef DEBUG + dX = sum(chi)+sum(trK)+sum(dxx)+sum(gxy)+sum(gxz)+sum(dyy)+sum(gyz)+sum(dzz) & + +sum(Axx)+sum(Axy)+sum(Axz)+sum(Ayy)+sum(Ayz)+sum(Azz) & + +sum(Gamx)+sum(Gamy)+sum(Gamz) & + +sum(Lap)+sum(betax)+sum(betay)+sum(betaz)+sum(dtSfx)+sum(dtSfy)+sum(dtSfz) & +sum(TZ) if(dX.ne.dX) then if(sum(chi).ne.sum(chi))write(*,*)"Z4c_rhs_ss.f90: find NaN in chi" @@ -152,10 +153,11 @@ if(sum(dtSfx).ne.sum(dtSfx))write(*,*)"Z4c_rhs_ss.f90: find NaN in dtSfx" if(sum(dtSfy).ne.sum(dtSfy))write(*,*)"Z4c_rhs_ss.f90: find NaN in dtSfy" if(sum(dtSfz).ne.sum(dtSfz))write(*,*)"Z4c_rhs_ss.f90: find NaN in dtSfz" - if(sum(TZ).ne.sum(Tz))write(*,*)"Z4c_rhs_ss.f90: find NaN in TZ" - gont = 1 - return - endif + if(sum(TZ).ne.sum(Tz))write(*,*)"Z4c_rhs_ss.f90: find NaN in TZ" + gont = 1 + return + endif +#endif PI = dacos(-ONE) @@ -1388,41 +1390,43 @@ call kodis_sh(ex,crho,sigma,R,TZ,TZ_rhs,SSS,Symmetry,eps,sst) endif -#if (ABV == 1) - call ricci_gamma_ss(ex,crho,sigma,R,X, Y, Z, & - drhodx, drhody, drhodz, & - dsigmadx,dsigmady,dsigmadz, & - dRdx,dRdy,dRdz, & - drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz, & - dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz, & - dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz, & - 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,Lev,sst) - call constraint_bssn_ss(ex,crho,sigma,R,X, Y, Z, & - drhodx, drhody, drhodz, & - dsigmadx,dsigmady,dsigmadz, & - dRdx,dRdy,dRdz, & - drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz, & - dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz, & - dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz, & - 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,Lev,sst) -#endif + if (co == 0) then +#if (ABV == 1) + call ricci_gamma_ss(ex,crho,sigma,R,X, Y, Z, & + drhodx, drhody, drhodz, & + dsigmadx,dsigmady,dsigmadz, & + dRdx,dRdy,dRdz, & + drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz, & + dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz, & + dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz, & + 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,Lev,sst) +#endif + call constraint_bssn_ss(ex,crho,sigma,R,X, Y, Z, & + drhodx, drhody, drhodz, & + dsigmadx,dsigmady,dsigmadz, & + dRdx,dRdy,dRdz, & + drhodxx,drhodxy,drhodxz,drhodyy,drhodyz,drhodzz, & + dsigmadxx,dsigmadxy,dsigmadxz,dsigmadyy,dsigmadyz,dsigmadzz, & + dRdxx,dRdxy,dRdxz,dRdyy,dRdyz,dRdzz, & + 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,Lev,sst) + endif gont = 0 diff --git a/AMSS_NCKU_source/makefile b/AMSS_NCKU_source/makefile index 72b9cbd..cbafb76 100644 --- a/AMSS_NCKU_source/makefile +++ b/AMSS_NCKU_source/makefile @@ -64,6 +64,9 @@ lopsided_c.o: lopsided_c.C lopsided_kodis_c.o: lopsided_kodis_c.C ${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@ +z4c_rhs_c.o: z4c_rhs_c.C + ${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@ + #interp_lb_profile.o: interp_lb_profile.C interp_lb_profile.h # ${CXX} $(CXXAPPFLAGS) -c $< $(filein) -o $@ @@ -90,6 +93,13 @@ else CFILES = bssn_rhs_c.o fderivs_c.o fdderivs_c.o kodiss_c.o lopsided_c.o lopsided_kodis_c.o endif +ifeq ($(USE_CXX_Z4C_KERNELS),1) +CFILES += z4c_rhs_c.o +Z4C_F90_OBJ = +else +Z4C_F90_OBJ = Z4c_rhs.o +endif + ## RK4 kernel switch (independent from USE_CXX_KERNELS) ifeq ($(USE_CXX_RK4),1) CFILES += rungekutta4_rout_c.o @@ -119,10 +129,10 @@ F90FILES_BASE = enforce_algebra.o fmisc.o initial_puncture.o prolongrestrict.o\ lopsidediff.o sommerfeld_rout.o getnp4.o diff_new_sh.o\ shellfunctions.o bssn_rhs_ss.o Set_Rho_ADM.o\ getnp4EScalar.o bssnEScalar_rhs.o bssn_constraint.o ricci_gamma.o\ - fadmquantites_bssn.o Z4c_rhs.o Z4c_rhs_ss.o point_diff_new_sh.o\ - cpbc.o getnp4old.o NullEvol.o initial_null.o initial_maxwell.o\ - getnpem2.o empart.o NullNews.o fourdcurvature.o\ - bssn2adm.o adm_constraint.o adm_ricci_gamma.o\ + fadmquantites_bssn.o $(Z4C_F90_OBJ) Z4c_rhs_ss.o point_diff_new_sh.o\ + cpbc.o getnp4old.o NullEvol.o initial_null.o initial_maxwell.o\ + getnpem2.o empart.o NullNews.o fourdcurvature.o\ + bssn2adm.o adm_constraint.o adm_ricci_gamma.o\ scalar_rhs.o initial_scalar.o NullEvol2.o initial_null2.o\ NullNews2.o tool_f.o diff --git a/AMSS_NCKU_source/makefile.inc b/AMSS_NCKU_source/makefile.inc index 331cff1..40a5c39 100755 --- a/AMSS_NCKU_source/makefile.inc +++ b/AMSS_NCKU_source/makefile.inc @@ -48,6 +48,11 @@ endif ## 0 : fall back to original Fortran kernels USE_CXX_KERNELS ?= 1 +## Z4C Cartesian RHS kernel switch +## 1 (default) : use C++ rewrite of Z4c_rhs (main Cartesian path faster) +## 0 : use original Fortran Z4c_rhs.o +USE_CXX_Z4C_KERNELS ?= 1 + ## RK4 kernel implementation switch ## 1 (default) : use C/C++ rewrite of rungekutta4_rout (for optimization experiments) ## 0 : use original Fortran rungekutta4_rout.o diff --git a/AMSS_NCKU_source/z4c_rhs_c.C b/AMSS_NCKU_source/z4c_rhs_c.C new file mode 100644 index 0000000..b56f821 --- /dev/null +++ b/AMSS_NCKU_source/z4c_rhs_c.C @@ -0,0 +1,661 @@ +#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 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) + { + 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], + 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 + + if (eps > 0.0) + { + lopsided_kodis(ex, X, Y, Z, gxx, gxx_rhs, betax, betay, betaz, Symmetry, SSS, eps); + lopsided_kodis(ex, X, Y, Z, gxy, gxy_rhs, betax, betay, betaz, Symmetry, AAS, eps); + lopsided_kodis(ex, X, Y, Z, gxz, gxz_rhs, betax, betay, betaz, Symmetry, ASA, eps); + lopsided_kodis(ex, X, Y, Z, gyy, gyy_rhs, betax, betay, betaz, Symmetry, SSS, eps); + lopsided_kodis(ex, X, Y, Z, gyz, gyz_rhs, betax, betay, betaz, Symmetry, SAA, eps); + lopsided_kodis(ex, X, Y, Z, gzz, gzz_rhs, betax, betay, betaz, Symmetry, SSS, eps); + + lopsided_kodis(ex, X, Y, Z, Axx, Axx_rhs, betax, betay, betaz, Symmetry, SSS, eps); + lopsided_kodis(ex, X, Y, Z, Axy, Axy_rhs, betax, betay, betaz, Symmetry, AAS, eps); + lopsided_kodis(ex, X, Y, Z, Axz, Axz_rhs, betax, betay, betaz, Symmetry, ASA, eps); + lopsided_kodis(ex, X, Y, Z, Ayy, Ayy_rhs, betax, betay, betaz, Symmetry, SSS, eps); + lopsided_kodis(ex, X, Y, Z, Ayz, Ayz_rhs, betax, betay, betaz, Symmetry, SAA, eps); + lopsided_kodis(ex, X, Y, Z, Azz, Azz_rhs, betax, betay, betaz, Symmetry, SSS, eps); + + lopsided_kodis(ex, X, Y, Z, chi_state, chi_rhs, betax, betay, betaz, Symmetry, SSS, eps); + lopsided_kodis(ex, X, Y, Z, trK, trK_rhs, betax, betay, betaz, Symmetry, SSS, eps); + + lopsided_kodis(ex, X, Y, Z, Gamx, Gamx_rhs, betax, betay, betaz, Symmetry, ASS, eps); + lopsided_kodis(ex, X, Y, Z, Gamy, Gamy_rhs, betax, betay, betaz, Symmetry, SAS, eps); + lopsided_kodis(ex, X, Y, Z, Gamz, Gamz_rhs, betax, betay, betaz, Symmetry, SSA, eps); + + lopsided_kodis(ex, X, Y, Z, Lap, Lap_rhs, betax, betay, betaz, Symmetry, SSS, eps); + lopsided_kodis(ex, X, Y, Z, betax, betax_rhs, betax, betay, betaz, Symmetry, ASS, eps); + lopsided_kodis(ex, X, Y, Z, betay, betay_rhs, betax, betay, betaz, Symmetry, SAS, eps); + lopsided_kodis(ex, X, Y, Z, betaz, betaz_rhs, betax, betay, betaz, Symmetry, SSA, eps); +#if (GAUGE == 0) + lopsided_kodis(ex, X, Y, Z, dtSfx, dtSfx_rhs, betax, betay, betaz, Symmetry, ASS, eps); + lopsided_kodis(ex, X, Y, Z, dtSfy, dtSfy_rhs, betax, betay, betaz, Symmetry, SAS, eps); + lopsided_kodis(ex, X, Y, Z, dtSfz, dtSfz_rhs, betax, betay, betaz, Symmetry, SSA, eps); +#endif + lopsided_kodis(ex, X, Y, Z, TZ, TZ_rhs, betax, betay, betaz, Symmetry, SSS, eps); + } + else + { + 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); + } + + 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; +}