Z4co: implement position depenent eta as in LazEv

EinsteinToolkit · Sep 10, 2024 · 4eff1c2 · 4eff1c2
1 parent 09eee68
commit 4eff1c2
Show file tree

Hide file tree

Showing 4 changed files with 30 additions and 1 deletion.
diff --git a/Z4co/param.ccl b/Z4co/param.ccl
@@ -49,6 +49,21 @@ CCTK_REAL eta "eta" STEERABLE=always
   *:* :: ""
 } 2.0
 
+CCTK_REAL veta_width "" STEERABLE=always
+{
+  (0:* :: ""
+} 50.0
+
+CCTK_REAL veta_central "" STEERABLE=always
+{
+  (0:* :: ""
+} 2.0
+
+CCTK_REAL veta_outer "" STEERABLE=always
+{
+  (0:* :: ""
+} 0.25
+
 CCTK_REAL chi_floor "Floor for chi" STEERABLE=always
 {
   (0:* :: ""

diff --git a/Z4co/src/rhs.cxx b/Z4co/src/rhs.cxx
@@ -226,7 +226,14 @@ extern "C" void Z4co_RHS(CCTK_ARGUMENTS) {
   const CCTK_REAL ckappa2 = kappa2;
   const CCTK_REAL cmuL = f_mu_L;
   const CCTK_REAL cmuS = f_mu_S;
-  const CCTK_REAL ceta = eta;
+  // const CCTK_REAL ceta = eta;
+  const auto calceta = [=](const vreal x, const CCTK_REAL y,
+                           const CCTK_REAL z) {
+    const vreal r = sqrt(x * x + y * y + z * z);
+    const CCTK_REAL is4 =
+        1.0 / (veta_width * veta_width * veta_width * veta_width);
+    return (veta_central - veta_outer) * exp(-r * r * r * r * is4) + veta_outer;
+  };
 
   // Loop
 #ifdef __CUDACC__

diff --git a/Z4co/wolfram/Z4co_set_rhs_with_adv.hxx b/Z4co/wolfram/Z4co_set_rhs_with_adv.hxx
@@ -34,6 +34,9 @@ noinline([&]() __attribute__((__flatten__, __hot__)) {
     const GF3D2index index2(layout2, p.I);
     const GF3D5index index5(layout5, p.I);
 
+const vreal x0 = p.x + Arith::iota<vreal>() * p.dx;
+const vreal ceta = calceta(x0, p.y, p.z);
+
 const auto &tmp_eTtt = gf_eTtt(mask, index2);
 const auto &tmp_eTt = gf_eTt(mask, index2);
 const auto &tmp_eT = gf_eT(mask, index2);

diff --git a/Z4co/wolfram/Z4co_set_rhs_with_adv.wl b/Z4co/wolfram/Z4co_set_rhs_with_adv.wl
@@ -47,6 +47,10 @@ $MainPrint[] :=
     pr["    const GF3D5index index5(layout5, p.I);"];
     pr[];
 
+    pr["const vreal x0 = p.x + Arith::iota<vreal>() * p.dx;"]
+    pr["const vreal ceta = calceta(x0, p.y, p.z);"];
+    pr[];
+
     PrintListInitializations[TmunuVarlist, "gf_", "index2"];
     PrintListInitializations[EvolVarlist, "tl_", "index5"];
     PrintListInitializations[dEvolVarlist, "tl_", "index5"];