Con2PrimFactory: Update C2Ps to use the modified EOS routines

jaykalinani · Jul 23, 2024 · 8534c86 · 8534c86
1 parent 5835c11
commit 8534c86
Show file tree

Hide file tree

Showing 6 changed files with 39 additions and 45 deletions.
diff --git a/Con2PrimFactory/interface.ccl b/Con2PrimFactory/interface.ccl
@@ -6,6 +6,5 @@ INCLUDES HEADER: prims.hxx IN prims.hxx
 INCLUDES HEADER: cons.hxx IN cons.hxx 
 INCLUDES HEADER: c2p_2DNoble.hxx IN c2p_2DNoble.hxx 
 INCLUDES HEADER: c2p_1DPalenzuela.hxx IN c2p_1DPalenzuela.hxx 
-USES INCLUDE HEADER: eos.hxx 
-USES INCLUDE HEADER: eos_idealgas.hxx
+USES INCLUDE HEADER: setup_eos.hxx 
 USES INCLUDE HEADER: roots.hxx
diff --git a/Con2PrimFactory/schedule.ccl b/Con2PrimFactory/schedule.ccl
@@ -1,7 +1,7 @@
 #Schedule definitions for thorn Con2PrimFactory
 
 if (unit_test) {
-  SCHEDULE Con2PrimFactory_Test AT wragh {
+  SCHEDULE Con2PrimFactory_Test AT CCTK_WRAGH AFTER EOSX_Setup_EOS {
   LANG:
     C OPTIONS : meta
   }

diff --git a/Con2PrimFactory/src/c2p.hxx b/Con2PrimFactory/src/c2p.hxx
@@ -13,8 +13,7 @@ c2p is effectively an interface to be used by different c2p implementations.
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
 #include "c2p_utils.hxx"
-#include "eos.hxx"
-#include "eos_idealgas.hxx"
+#include "setup_eos.hxx"
 
 #include <math.h>
 #include "prims.hxx"

diff --git a/Con2PrimFactory/src/c2p_1DPalenzuela.hxx b/Con2PrimFactory/src/c2p_1DPalenzuela.hxx
@@ -13,7 +13,7 @@ public:
   /* Constructor */
   template <typename EOSType>
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DPalenzuela(
-      EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
+      EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
       CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len);
 
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
@@ -34,18 +34,18 @@ public:
   template <typename EOSType>
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
   xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq, CCTK_REAL Bsq,
-                    CCTK_REAL BiSi, EOSType &eos_th, prim_vars &pv,
+                    CCTK_REAL BiSi, EOSType &eos_3p, prim_vars &pv,
                     cons_vars cv, const smat<CCTK_REAL, 3> &gup,
                     const smat<CCTK_REAL, 3> &glo) const;
 
   template <typename EOSType>
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
   funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq, CCTK_REAL BiSi,
-                        CCTK_REAL x, EOSType &eos_th, cons_vars &cv) const;
+                        CCTK_REAL x, EOSType &eos_3p, cons_vars &cv) const;
 
   template <typename EOSType>
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
-  solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv,
+  solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv,
         const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const;
 
   /* Destructor */
@@ -56,10 +56,10 @@ public:
 template <typename EOSType>
 CCTK_HOST CCTK_DEVICE
     CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_1DPalenzuela::c2p_1DPalenzuela(
-        EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
+        EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
         CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len) {
 
-  GammaIdealFluid = eos_th.gamma;
+  GammaIdealFluid = eos_3p.gamma;
   maxIterations = maxIter;
   tolerance = tol;
   rho_strict = rho_str;
@@ -131,7 +131,7 @@ template <typename EOSType>
 CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
 c2p_1DPalenzuela::xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq,
                                     CCTK_REAL Bsq, CCTK_REAL BiSi,
-                                    EOSType &eos_th, prim_vars &pv,
+                                    EOSType &eos_3p, prim_vars &pv,
                                     cons_vars cv, const smat<CCTK_REAL, 3> &gup,
                                     const smat<CCTK_REAL, 3> &glo) const {
   const CCTK_REAL qPalenzuela = cv.tau / cv.dens;
@@ -184,7 +184,7 @@ c2p_1DPalenzuela::xPalenzuelaToPrim(CCTK_REAL xPalenzuela_Sol, CCTK_REAL Ssq,
 
   pv.Ye = cv.dYe / cv.dens;
 
-  pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+  pv.press = eos_3p.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
 
   pv.Bvec = cv.dBvec;
 
@@ -196,7 +196,7 @@ template <typename EOSType>
 CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
 c2p_1DPalenzuela::funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq,
                                         CCTK_REAL BiSi, CCTK_REAL x,
-                                        EOSType &eos_th, cons_vars &cv) const {
+                                        EOSType &eos_3p, cons_vars &cv) const {
   // computes f(x) from x and q,r,s,t
   const CCTK_REAL qPalenzuela = cv.tau / cv.dens;
   const CCTK_REAL rPalenzuela = Ssq / pow(cv.dens, 2);
@@ -223,14 +223,14 @@ c2p_1DPalenzuela::funcRoot_1DPalenzuela(CCTK_REAL Ssq, CCTK_REAL Bsq,
 
   // (iv)
   CCTK_REAL P_loc =
-      eos_th.press_from_valid_rho_eps_ye(rho_loc, eps_loc, Ye_loc);
+      eos_3p.press_from_valid_rho_eps_ye(rho_loc, eps_loc, Ye_loc);
 
   return (x - (1.0 + eps_loc + P_loc / rho_loc) * W_loc);
 }
 
 template <typename EOSType>
 CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
-c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
+c2p_1DPalenzuela::solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds,
                         cons_vars cv, const smat<CCTK_REAL, 3> &glo,
                         c2p_report &rep) const {
 
@@ -311,7 +311,7 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
   CCTK_REAL a = xPalenzuela_lowerBound;
   CCTK_REAL b = xPalenzuela_upperBound;
   auto fn = [&](auto x) {
-    return funcRoot_1DPalenzuela(Ssq, Bsq, BiSi, x, eos_th, cv);
+    return funcRoot_1DPalenzuela(Ssq, Bsq, BiSi, x, eos_3p, cv);
   };
   auto result = Algo::brent(fn, a, b, minbits, maxiters, rep.iters);
 
@@ -334,7 +334,7 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
     status = ROOTSTAT::NOT_CONVERGED;
   }
 
-  xPalenzuelaToPrim(xPalenzuela_Sol, Ssq, Bsq, BiSi, eos_th, pv, cv, gup, glo);
+  xPalenzuelaToPrim(xPalenzuela_Sol, Ssq, Bsq, BiSi, eos_3p, pv, cv, gup, glo);
 
   // set to atmo if computed rho is below floor density
   if (pv.rho < atmo.rho_cut) {
@@ -344,7 +344,7 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
   }
 
   // check the validity of the computed eps
-  auto rgeps = eos_th.range_eps_from_valid_rho_ye(pv.rho, pv.Ye);
+  auto rgeps = eos_3p.range_eps_from_valid_rho_ye(pv.rho, pv.Ye);
   if (pv.eps > rgeps.max) {
     //printf("(pv.eps > rgeps.max) is true, adjusting cons.. \n");
     rep.adjust_cons = true;
@@ -383,8 +383,8 @@ c2p_1DPalenzuela::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
     }
     pv.vel *= v_lim / sol_v;
     pv.w_lor = w_lim;
-    pv.eps = std::min(std::max(eos_th.rgeps.min, pv.eps), eos_th.rgeps.max);
-    pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+    pv.eps = std::min(std::max(eos_3p.rgeps.min, pv.eps), eos_3p.rgeps.max);
+    pv.press = eos_3p.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
 
     rep.adjust_cons = true;
   }

diff --git a/Con2PrimFactory/src/c2p_2DNoble.hxx b/Con2PrimFactory/src/c2p_2DNoble.hxx
@@ -15,7 +15,7 @@ public:
   /* Constructor */
   template <typename EOSType>
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_2DNoble(
-      EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
+      EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
       CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len);
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline CCTK_REAL
       get_Ssq_Exact(vec<CCTK_REAL, 3> &mom,
@@ -45,11 +45,11 @@ public:
   template <typename EOSType>
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
   WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq, CCTK_REAL BiSi,
-          EOSType &eos_th, prim_vars &pv, cons_vars cv,
+          EOSType &eos_3p, prim_vars &pv, cons_vars cv,
           const smat<CCTK_REAL, 3> &gup) const;
   template <typename EOSType>
   CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
-  solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv,
+  solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds, cons_vars cv,
         const smat<CCTK_REAL, 3> &glo, c2p_report &rep) const;
 
   /* Destructor */
@@ -60,10 +60,10 @@ public:
 template <typename EOSType>
 CCTK_HOST
     CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline c2p_2DNoble::c2p_2DNoble(
-        EOSType &eos_th, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
+        EOSType &eos_3p, atmosphere &atm, CCTK_INT maxIter, CCTK_REAL tol,
         CCTK_REAL rho_str, CCTK_REAL vwlim, CCTK_REAL B_lim, bool ye_len) {
 
-  GammaIdealFluid = eos_th.gamma;
+  GammaIdealFluid = eos_3p.gamma;
   maxIterations = maxIter;
   tolerance = tol;
   rho_strict = rho_str;
@@ -151,7 +151,7 @@ c2p_2DNoble::get_dPdVsq_funcZVsq(CCTK_REAL Z, CCTK_REAL Vsq,
 template <typename EOSType>
 CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
 c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq,
-                     CCTK_REAL BiSi, EOSType &eos_th, prim_vars &pv,
+                     CCTK_REAL BiSi, EOSType &eos_3p, prim_vars &pv,
                      cons_vars cv, const smat<CCTK_REAL, 3> &gup) const {
   CCTK_REAL W_Sol = 1.0 / sqrt(1.0 - vsq_Sol);
 
@@ -178,7 +178,7 @@ c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq,
 
   pv.Ye = cv.dYe / cv.dens;
 
-  pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+  pv.press = eos_3p.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
 
   pv.Bvec = cv.dBvec;
 
@@ -188,7 +188,7 @@ c2p_2DNoble::WZ2Prim(CCTK_REAL Z_Sol, CCTK_REAL vsq_Sol, CCTK_REAL Bsq,
 
 template <typename EOSType>
 CCTK_HOST CCTK_DEVICE CCTK_ATTRIBUTE_ALWAYS_INLINE inline void
-c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
+c2p_2DNoble::solve(EOSType &eos_3p, prim_vars &pv, prim_vars &pv_seeds,
                    cons_vars cv, const smat<CCTK_REAL, 3> &glo,
                    c2p_report &rep) const {
 
@@ -260,12 +260,12 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
   pv_seeds.rho = cv.dens / pv_seeds.w_lor;
 
   CCTK_REAL eps_min =
-      eos_th.range_eps_from_valid_rho_ye(pv_seeds.rho, pv_seeds.Ye).min;
+      eos_3p.range_eps_from_valid_rho_ye(pv_seeds.rho, pv_seeds.Ye).min;
   CCTK_REAL eps_last = max({pv_seeds.eps, eps_min});
 
   /* get pressure seed from updated pv_seeds.rho */
   pv_seeds.press =
-      eos_th.press_from_valid_rho_eps_ye(pv_seeds.rho, eps_last, pv_seeds.Ye);
+      eos_3p.press_from_valid_rho_eps_ye(pv_seeds.rho, eps_last, pv_seeds.Ye);
 
   /* get Z seed */
   CCTK_REAL Z_Seed =
@@ -387,7 +387,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
   CCTK_REAL vsq_Sol = x[1];
 
   /* Write prims if C2P succeeded */
-  WZ2Prim(Z_Sol, vsq_Sol, Bsq, BiSi, eos_th, pv, cv, gup);
+  WZ2Prim(Z_Sol, vsq_Sol, Bsq, BiSi, eos_3p, pv, cv, gup);
 
   // set to atmo if computed rho is below floor density
   if (pv.rho < atmo.rho_cut) {
@@ -397,7 +397,7 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
   }
 
   // check the validity of the computed eps
-  auto rgeps = eos_th.range_eps_from_valid_rho_ye(pv.rho, pv.Ye);
+  auto rgeps = eos_3p.range_eps_from_valid_rho_ye(pv.rho, pv.Ye);
   if (pv.eps > rgeps.max) {
     //printf("(pv.eps > rgeps.max) is true, adjusting cons.. \n");
     rep.adjust_cons = true;
@@ -434,8 +434,8 @@ c2p_2DNoble::solve(EOSType &eos_th, prim_vars &pv, prim_vars &pv_seeds,
     }
     pv.vel *= v_lim / sol_v;
     pv.w_lor = w_lim;
-    pv.eps = std::min(std::max(eos_th.rgeps.min, pv.eps), eos_th.rgeps.max);
-    pv.press = eos_th.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
+    pv.eps = std::min(std::max(eos_3p.rgeps.min, pv.eps), eos_3p.rgeps.max);
+    pv.press = eos_3p.press_from_valid_rho_eps_ye(pv.rho, pv.eps, pv.Ye);
 
     rep.adjust_cons = true;
   }

diff --git a/Con2PrimFactory/src/test.cxx b/Con2PrimFactory/src/test.cxx
@@ -7,8 +7,8 @@
 #include "c2p.hxx"
 #include "c2p_1DPalenzuela.hxx"
 #include "c2p_2DNoble.hxx"
-#include <eos.hxx>
-#include <eos_idealgas.hxx>
+
+#include "setup_eos.hxx"
 
 namespace Con2PrimFactory {
 
@@ -18,10 +18,6 @@ using namespace EOSX;
 extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTS;
 
-  // Initializing the IG EOS
-  eos::range rgeps(0, 100), rgrho(1e-13, 20), rgye(0.5, 0.5);
-  const eos_idealgas eos_th(2.0, 938.985, rgeps, rgrho, rgye);
-
   // Setting up atmosphere
   atmosphere atmo(1e-10, 1e-8, 0.5, 1e-8, 0.001);
 
@@ -30,8 +26,8 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) {
                              1.0, 0.0, 1.0}; // xx, xy, xz, yy, yz, zz
 
   // Con2Prim objects
-  c2p_2DNoble c2p_Noble(eos_th, atmo, 100, 1e-8, 1e8, 1, 1, true);
-  c2p_1DPalenzuela c2p_Pal(eos_th, atmo, 100, 1e-8, 1e8, 1, 1, true);
+  c2p_2DNoble c2p_Noble(*eos_3p_ig, atmo, 100, 1e-8, 1e8, 1, 1, true);
+  c2p_1DPalenzuela c2p_Pal(*eos_3p_ig, atmo, 100, 1e-8, 1e8, 1, 1, true);
 
   // Construct error report object:
   c2p_report rep_Noble;
@@ -49,7 +45,7 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) {
 
   // Testing C2P Noble
   CCTK_VINFO("Testing C2P Noble...");
-  c2p_Noble.solve(eos_th, pv, pv_seeds, cv, g, rep_Noble);
+  c2p_Noble.solve(*eos_3p_ig, pv, pv_seeds, cv, g, rep_Noble);
 
   printf("pv_seeds, pv: \n"
          "rho: %f, %f \n"
@@ -90,7 +86,7 @@ extern "C" void Con2PrimFactory_Test(CCTK_ARGUMENTS) {
 
   // Testing C2P Palenzuela
   CCTK_VINFO("Testing C2P Palenzuela...");
-  c2p_Pal.solve(eos_th, pv, pv_seeds, cv, g, rep_Pal);
+  c2p_Pal.solve(*eos_3p_ig, pv, pv_seeds, cv, g, rep_Pal);
 
   printf("pv_seeds, pv: \n"
          "rho: %f, %f \n"