Merge pull request #61 from jaykalinani/z_vec

Merge into main

AsterX/interface.ccl

@@ -99,6 +99,16 @@ CCTK_REAL saved_prims TYPE=gf CENTERING={ccc}
   saved_eps
 } "Saved primitive variables as initial guesses for con2prim"
 
+CCTK_REAL zvec TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"'
+{
+  zvec_x, zvec_y, zvec_z
+} ""
+
+CCTK_REAL svec TYPE=gf CENTERING={ccc} TAGS='checkpoint="no"'
+{
+  svec_x, svec_y, svec_z
+} ""
+
 #grid functions required for upwindCT computation of electric field
 
 CCTK_REAL vtilde_xface TYPE=gf CENTERING={vcc} TAGS='checkpoint="no"'

AsterX/param.ccl

@@ -21,6 +21,13 @@ KEYWORD flux_type "Flux solver" STEERABLE=always
   "HLLE" :: ""
 } "LxF"
 
+KEYWORD recon_type "Vector field used in reconstruction" STEERABLE=always
+{
+  "v_vec" :: "Eulerian 3-velocity"
+  "z_vec" :: "Lorentz factor x Eulerian 3-velocity"
+  "s_vec" :: "Purely hydrodynamic part of conserved momentum without volume factor"
+} "v_vec"
+
 BOOLEAN local_estimate_error "Use error estimation criteria of this thorn" STEERABLE=always
 {
 } yes

AsterX/schedule.ccl

@@ -72,8 +72,12 @@ SCHEDULE AsterX_Prim2Con_Initial IN AsterX_InitialGroup AFTER AsterX_ComputeBFro
   WRITES: dens(interior) tau(interior) mom(interior) dB(interior)
   WRITES: Psi(everywhere)
   WRITES: saved_prims
+  WRITES: zvec
+  WRITES: svec
   SYNC: dens tau mom dB
   SYNC: saved_prims
+  SYNC: zvec
+  SYNC: svec
 } "Compute conserved variables from primitive variables at initial"
 
 
@@ -125,10 +129,14 @@ SCHEDULE AsterX_Con2Prim IN AsterX_Con2PrimGroup AFTER AsterX_ComputedBFromdBsta
   WRITES: con2prim_flag(interior)
   WRITES: HydroBaseX::rho(interior) HydroBaseX::vel(interior) HydroBaseX::eps(interior) HydroBaseX::press(interior) HydroBaseX::Bvec(interior)
   WRITES: saved_prims(interior)
+  WRITES: zvec(interior)
+  WRITES: svec(interior)
   WRITES: dens(interior) tau(interior) mom(interior) dB(interior)
   SYNC: con2prim_flag
   SYNC: HydroBaseX::rho HydroBaseX::vel HydroBaseX::eps HydroBaseX::press HydroBaseX::Bvec
   SYNC: saved_prims
+  SYNC: zvec
+  SYNC: svec
   SYNC: dens tau mom dB
 } "Calculate primitive variables from conservative variables"
 
@@ -148,6 +156,8 @@ SCHEDULE AsterX_Fluxes IN AsterX_RHSGroup
   READS: HydroBaseX::rho(everywhere) HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere)
   READS: HydroBaseX::Bvec(everywhere)
   READS: dBx_stag(everywhere) dBy_stag(everywhere) dBz_stag(everywhere)
+  READS: zvec_x(everywhere) zvec_y(everywhere) zvec_z(everywhere)
+  READS: svec_x(everywhere) svec_y(everywhere) svec_z(everywhere)
   READS: Avec_x(everywhere) Avec_y(everywhere) Avec_z(everywhere) Psi(everywhere)
   WRITES: flux_x(interior) flux_y(interior) flux_z(interior)
   WRITES: Aux_in_RHSof_A_Psi(interior)
@@ -165,7 +175,10 @@ SCHEDULE AsterX_SourceTerms IN AsterX_RHSGroup AFTER AsterX_Fluxes
   READS: ADMBaseX::lapse(everywhere)
   READS: ADMBaseX::shift(everywhere)
   READS: ADMBaseX::curv(everywhere)
-  READS: HydroBaseX::rho(everywhere) HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere)
+  READS: HydroBaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere)
+  READS: HydroBaseX::vel(everywhere)
+  READS: zvec_x(everywhere), zvec_y(everywhere), zvec_z(everywhere)
+  READS: svec_x(everywhere), svec_y(everywhere), svec_z(everywhere)
   READS: HydroBaseX::Bvec(everywhere)
   WRITES: densrhs(interior) taurhs(interior) momrhs(interior)
   SYNC: densrhs taurhs momrhs
@@ -217,7 +230,10 @@ if(update_tmunu){
   {
     LANG: C
     READS: ADMBaseX::metric(everywhere) ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere)
-    READS: HydroBaseX::rho(everywhere) HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere)
+    READS: HydroBaseX::rho(everywhere) HydroBaseX::press(everywhere) HydroBaseX::eps(everywhere)
+    READS: HydroBaseX::vel(everywhere)
+    READS: zvec_x(everywhere), zvec_y(everywhere), zvec_z(everywhere)
+    READS: svec_x(everywhere), svec_y(everywhere), svec_z(everywhere)
     READS: HydroBaseX::Bvec(everywhere)
     READS: TmunuBaseX::eTtt(interior) TmunuBaseX::eTti(interior) TmunuBaseX::eTij(interior)
     WRITES: TmunuBaseX::eTtt(interior) TmunuBaseX::eTti(interior) TmunuBaseX::eTij(interior)

AsterX/src/con2prim.cxx

@@ -219,6 +219,14 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold,
     pv.scatter(rho(p.I), eps(p.I), dummy_Ye, press(p.I), velx(p.I), vely(p.I),
                velz(p.I), wlor, Bvecx(p.I), Bvecy(p.I), Bvecz(p.I), Ex, Ey, Ez);
 
+    zvec_x(p.I) = wlor * pv.vel(0); 
+    zvec_y(p.I) = wlor * pv.vel(1);
+    zvec_z(p.I) = wlor * pv.vel(2);
+
+    svec_x(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(0); 
+    svec_y(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(1);
+    svec_z(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(2);
+
     // Write back cv
     cv.scatter(dens(p.I), momx(p.I), momy(p.I), momz(p.I), tau(p.I), dummy_Ye,
                dBx(p.I), dBy(p.I), dBz(p.I));

AsterX/src/con2prim_rpa.cxx

@@ -154,6 +154,14 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
     pv.scatter(rho(p.I), eps(p.I), dumye, press(p.I), velx(p.I), vely(p.I),
                velz(p.I), wlor, Ex, Ey, Ez, Bvecx(p.I), Bvecy(p.I), Bvecz(p.I));
 
+    zvec_x(p.I) = wlor * pv.vel(0); 
+    zvec_y(p.I) = wlor * pv.vel(1);
+    zvec_z(p.I) = wlor * pv.vel(2);
+
+    svec_x(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(0); 
+    svec_y(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(1);
+    svec_z(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(2);
+
     // Write back cv
     if (rep.adjust_cons) {
       cv.scatter(dens(p.I), tau(p.I), dumye, momx(p.I), momy(p.I), momz(p.I),

AsterX/src/fluxes.cxx

@@ -25,6 +25,7 @@ using namespace ReconX;
 
 enum class flux_t { LxF, HLLE };
 enum class eos_t { IdealGas, Hybrid, Tabulated };
+enum class rec_var_t { v_vec, z_vec, s_vec };
 
 // Calculate the fluxes in direction `dir`. This function is more
 // complex because it has to handle any direction, but as reward,
@@ -45,6 +46,8 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
   const vec<GF3D2<CCTK_REAL>, dim> fluxBzs{fxBz, fyBz, fzBz};
   /* grid functions */
   const vec<GF3D2<const CCTK_REAL>, dim> gf_vels{velx, vely, velz};
+  const vec<GF3D2<const CCTK_REAL>, dim> gf_zvec{zvec_x, zvec_y, zvec_z};
+  const vec<GF3D2<const CCTK_REAL>, dim> gf_svec{svec_x, svec_y, svec_z};
   const vec<GF3D2<const CCTK_REAL>, dim> gf_Bvecs{Bvecx, Bvecy, Bvecz};
   const vec<GF3D2<const CCTK_REAL>, dim> gf_dBstags{dBx_stag, dBy_stag, dBz_stag};
   const vec<GF3D2<const CCTK_REAL>, dim> gf_beta{betax, betay, betaz};
@@ -59,6 +62,17 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
 
   static_assert(dir >= 0 && dir < 3, "");
 
+  rec_var_t rec_var;
+  if (CCTK_EQUALS(recon_type, "v_vec")) {
+    rec_var = rec_var_t::v_vec;
+  } else if (CCTK_EQUALS(recon_type, "z_vec")) {
+    rec_var = rec_var_t::z_vec;
+  } else if (CCTK_EQUALS(recon_type, "s_vec")) {
+    rec_var = rec_var_t::s_vec;
+  } else {
+    CCTK_ERROR("Unknown value for parameter \"recon_type\"");
+  }
+
   reconstruction_t reconstruction;
   if (CCTK_EQUALS(reconstruction_method, "Godunov"))
     reconstruction = reconstruction_t::Godunov;
@@ -174,10 +188,23 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
           [2]>(grid.nghostzones, [=] CCTK_DEVICE(
                                      const PointDesc
                                          &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+
     /* Reconstruct primitives from the cells on left (indice 0) and right
      * (indice 1) side of this face rc = reconstructed variables or
      * computed from reconstructed variables */
 
+    /* Interpolate metric components from vertices to faces */
+    const CCTK_REAL alp_avg = calc_avg_v2f(alp, p, dir);
+    const vec<CCTK_REAL, 3> betas_avg(
+        [&](int i) ARITH_INLINE { return calc_avg_v2f(gf_beta(i), p, dir); });
+    const smat<CCTK_REAL, 3> g_avg([&](int i, int j) ARITH_INLINE {
+      return calc_avg_v2f(gf_g(i, j), p, dir);
+    });
+
+    /* determinant of spatial metric */
+    const CCTK_REAL detg_avg = calc_det(g_avg);
+    const CCTK_REAL sqrtg = sqrt(detg_avg);
+
     vec<CCTK_REAL, 2> rho_rc{reconstruct_pt(rho, p, true, true)};
 
     // set to atmo if reconstructed rho is less than atmo or is negative
@@ -188,9 +215,6 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
       rho_rc(1) = rho_abs_min;
     }
 
-    const vec<vec<CCTK_REAL, 2>, 3> vels_rc([&](int i) ARITH_INLINE {
-      return vec<CCTK_REAL, 2>{reconstruct_pt(gf_vels(i), p, false, false)};
-    });
     vec<CCTK_REAL, 2> press_rc{reconstruct_pt(press, p, false, true)};
     // TODO: Correctly reconstruct Ye
     const vec<CCTK_REAL, 2> ye_rc{ye_min, ye_max};
@@ -206,28 +230,13 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
           eos_th.press_from_valid_rho_eps_ye(rho_rc(1), eps_min, ye_rc(1));
     }
 
-    /* Interpolate metric components from vertices to faces */
-    const CCTK_REAL alp_avg = calc_avg_v2f(alp, p, dir);
-    const vec<CCTK_REAL, 3> betas_avg(
-        [&](int i) ARITH_INLINE { return calc_avg_v2f(gf_beta(i), p, dir); });
-    const smat<CCTK_REAL, 3> g_avg([&](int i, int j) ARITH_INLINE {
-      return calc_avg_v2f(gf_g(i, j), p, dir);
+    const vec<CCTK_REAL, 2> eps_rc([&](int f) ARITH_INLINE {
+      return eos_th.eps_from_valid_rho_press_ye(rho_rc(f), press_rc(f),
+                                                ye_rc(f));
     });
 
-    /* determinant of spatial metric */
-    const CCTK_REAL detg_avg = calc_det(g_avg);
-    const CCTK_REAL sqrtg = sqrt(detg_avg);
-    /* co-velocity measured by Eulerian observer: v_j */
-    const vec<vec<CCTK_REAL, 2>, 3> vlows_rc = calc_contraction(g_avg, vels_rc);
-    /* vtilde^i = alpha * v^i - beta^i */
-    const vec<vec<CCTK_REAL, 2>, 3> vtildes_rc([&](int i) ARITH_INLINE {
-      return vec<CCTK_REAL, 2>([&](int f) ARITH_INLINE {
-        return alp_avg * vels_rc(i)(f) - betas_avg(i);
-      });
-    });
-    /* Lorentz factor: W = 1 / sqrt(1 - v^2) */
-    const vec<CCTK_REAL, 2> w_lorentz_rc([&](int f) ARITH_INLINE {
-      return 1 / sqrt(1 - calc_contraction(vlows_rc, vels_rc)(f));
+    const vec<CCTK_REAL, 2> rhoh_rc([&](int f) ARITH_INLINE {
+      return rho_rc(f) + rho_rc(f)*eps_rc(f) + press_rc(f);
     });
 
     // Introduce reconstructed Bs
@@ -249,6 +258,78 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
 
     // End of setting Bs
 
+    vec<vec<CCTK_REAL, 2>, 3> vels_rc;
+    vec<vec<CCTK_REAL, 2>, 3> vlows_rc;
+    vec<CCTK_REAL, 2> w_lorentz_rc;
+    array<CCTK_REAL,2> vels_rc_dummy; // note: can't copy array<,2> to vec<,2>, only construct
+    switch (rec_var) {
+    case rec_var_t::v_vec : {
+
+      for(int i = 0; i <= 2; ++i) { // loop over components
+	      vels_rc_dummy = reconstruct_pt(gf_vels(i), p, false, false);
+	      vels_rc(i)(0) = vels_rc_dummy[0];
+	      vels_rc(i)(1) = vels_rc_dummy[1];
+      }
+
+      /* co-velocity measured by Eulerian observer: v_j */
+      vlows_rc = calc_contraction(g_avg, vels_rc);
+
+      /* Lorentz factor: W = 1 / sqrt(1 - v^2) */
+      w_lorentz_rc(0) = 1 / sqrt(1 - calc_contraction(vlows_rc, vels_rc)(0));
+      w_lorentz_rc(1) = 1 / sqrt(1 - calc_contraction(vlows_rc, vels_rc)(1));
+      break;
+
+    };
+    case rec_var_t::z_vec : {
+
+      const vec<vec<CCTK_REAL, 2>, 3> zvec_rc([&](int i) ARITH_INLINE {
+         return vec<CCTK_REAL, 2>{reconstruct_pt(gf_zvec(i), p, false, false)};
+      });
+
+      const vec<vec<CCTK_REAL, 2>, 3> zveclow_rc = calc_contraction(g_avg, zvec_rc);
+
+      w_lorentz_rc(0) = sqrt(1 + calc_contraction(zveclow_rc, zvec_rc)(0));
+      w_lorentz_rc(1) = sqrt(1 + calc_contraction(zveclow_rc, zvec_rc)(1));
+
+      for(int i = 0; i <= 2; ++i) { // loop over components
+	      for(int j = 0; j <= 1; ++j) { // loop over left and right state
+		      vels_rc(i)(j) = zvec_rc(i)(j)/w_lorentz_rc(j);
+		      vlows_rc(i)(j) = zveclow_rc(i)(j)/w_lorentz_rc(j);
+	       }
+      }
+      break;
+
+    };
+    case rec_var_t::s_vec : {
+
+      const vec<vec<CCTK_REAL, 2>, 3> svec_rc([&](int i) ARITH_INLINE {
+         return vec<CCTK_REAL, 2>{reconstruct_pt(gf_svec(i), p, false, false)};
+      });
+
+      const vec<vec<CCTK_REAL, 2>, 3> sveclow_rc = calc_contraction(g_avg, svec_rc);
+
+      w_lorentz_rc(0) = sqrt(0.5+sqrt(0.25+calc_contraction(sveclow_rc, svec_rc)(0)/rhoh_rc(0)/rhoh_rc(0)));
+      w_lorentz_rc(1) = sqrt(0.5+sqrt(0.25+calc_contraction(sveclow_rc, svec_rc)(1)/rhoh_rc(1)/rhoh_rc(1)));
+
+      //printf("  wlor = %16.8e, %16.8e\n", w_lorentz_rc(0), w_lorentz_rc(1));
+
+      for(int i = 0; i <= 2; ++i) { // loop over components
+	      for(int j = 0; j <= 1; ++j) { // loop over left and right state
+		      vels_rc(i)(j) = svec_rc(i)(j)/w_lorentz_rc(j)/w_lorentz_rc(j)/rhoh_rc(j);
+		      vlows_rc(i)(j) = sveclow_rc(i)(j)/w_lorentz_rc(j)/w_lorentz_rc(j)/rhoh_rc(j);
+	       }
+      }
+      break;
+
+    };
+    }
+
+    /* vtilde^i = alpha * v^i - beta^i */
+    const vec<vec<CCTK_REAL, 2>, 3> vtildes_rc([&](int i) ARITH_INLINE {
+      return vec<CCTK_REAL, 2>([&](int f) ARITH_INLINE {
+        return alp_avg * vels_rc(i)(f) - betas_avg(i);
+      });
+    });
 
     /* alpha * b0 = W * B^i * v_i */
     const vec<CCTK_REAL, 2> alp_b0_rc([&](int f) ARITH_INLINE {
@@ -281,11 +362,6 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
     // vars
 
     // Ideal gas case {
-    /* eps for ideal gas EOS */
-    const vec<CCTK_REAL, 2> eps_rc([&](int f) ARITH_INLINE {
-      return eos_th.eps_from_valid_rho_press_ye(rho_rc(f), press_rc(f),
-                                                ye_rc(f));
-    });
     /* cs2 for ideal gas EOS */
     const vec<CCTK_REAL, 2> cs2_rc([&](int f) ARITH_INLINE {
       return eos_th.csnd_from_valid_rho_eps_ye(rho_rc(f), eps_rc(f), ye_rc(f)) *

AsterX/src/prim2con.cxx

@@ -50,6 +50,19 @@ extern "C" void AsterX_Prim2Con_Initial(CCTK_ARGUMENTS) {
         saved_vely(p.I) = pv.vel(1);
         saved_velz(p.I) = pv.vel(2);
         saved_eps(p.I) = pv.eps;
+
+	const vec<CCTK_REAL, 3> v_up{pv.vel(0),pv.vel(1),pv.vel(2)};
+        const vec<CCTK_REAL, 3> v_low = calc_contraction(g, v_up);
+        CCTK_REAL wlor = calc_wlorentz(v_low, v_up);
+
+	zvec_x(p.I) = wlor * pv.vel(0);
+	zvec_y(p.I) = wlor * pv.vel(1);
+	zvec_z(p.I) = wlor * pv.vel(2);
+
+	svec_x(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(0);
+	svec_y(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(1);
+	svec_z(p.I) = (pv.rho+pv.rho*pv.eps+pv.press)*wlor*wlor*pv.vel(2);
+
       });
 
   /* Initilaize Psi to 0.0 */

AsterX/src/rhs.cxx

@@ -196,6 +196,7 @@ extern "C" void AsterX_RHS(CCTK_ARGUMENTS) {
         momyrhs(p.I) += calcupdate_hydro(gf_fmomy, p);
         momzrhs(p.I) += calcupdate_hydro(gf_fmomz, p);
         taurhs(p.I) += calcupdate_hydro(gf_ftau, p);
+
         if (isnan(densrhs(p.I))) {
           printf("calcupdate = %f, ", calcupdate_hydro(gf_fdens, p));
           printf("densrhs = %f, gf_fdens = %f, %f, %f, %f, %f, %f \n",
@@ -204,6 +205,7 @@ extern "C" void AsterX_RHS(CCTK_ARGUMENTS) {
                  gf_fdens(1)(p.I + p.DI[1]), gf_fdens(2)(p.I + p.DI[2]));
         }
         assert(!isnan(densrhs(p.I)));
+
       });
 
   grid.loop_int_device<1, 0, 0>(grid.nghostzones,

AsterX/src/source.cxx

@@ -19,6 +19,8 @@ template <int FDORDER> void SourceTerms(CCTK_ARGUMENTS) {
   DECLARE_CCTK_ARGUMENTSX_AsterX_SourceTerms;
   DECLARE_CCTK_PARAMETERS;
 
+  const bool use_v_vec = CCTK_EQUALS(recon_type, "v_vec");
+
   /* grid functions */
   const vec<GF3D2<const CCTK_REAL>, 3> gf_beta{betax, betay, betaz};
   const smat<GF3D2<const CCTK_REAL>, 3> gf_g{gxx, gxy, gxz, gyy, gyz, gzz};
@@ -66,9 +68,27 @@ template <int FDORDER> void SourceTerms(CCTK_ARGUMENTS) {
         });
 
         /* Computing v_j */
-        const vec<CCTK_REAL, 3> v_up{velx(p.I), vely(p.I), velz(p.I)};
-        const vec<CCTK_REAL, 3> v_low = calc_contraction(g_avg, v_up);
-        const CCTK_REAL w_lorentz = calc_wlorentz(v_low, v_up);
+        vec<CCTK_REAL, 3> v_up;
+        vec<CCTK_REAL, 3> v_low;
+        CCTK_REAL w_lorentz;
+	if (use_v_vec) {
+
+            v_up(0) = velx(p.I);
+	    v_up(1) = vely(p.I);
+	    v_up(2) = velz(p.I);
+            v_low = calc_contraction(g_avg, v_up);
+            w_lorentz = calc_wlorentz(v_low, v_up);
+
+	} else {
+
+            const vec<CCTK_REAL, 3> z_up{zvec_x(p.I), zvec_y(p.I), zvec_z(p.I)};
+            const vec<CCTK_REAL, 3> z_low = calc_contraction(g_avg, z_up);
+            w_lorentz = calc_wlorentz_zvec(z_low, z_up);
+            v_up  = z_up/w_lorentz;
+            v_low = z_low/w_lorentz;
+
+	}
+
         /* Computing [ \rho(1+\epsilon) + Pgas ]*W^2 = \rho * h * W^2 */
         const CCTK_REAL rhoenthalpyW2 =
             (rho(p.I) * (1.0 + eps(p.I)) + press(p.I)) * w_lorentz * w_lorentz;
@@ -141,6 +161,7 @@ template <int FDORDER> void SourceTerms(CCTK_ARGUMENTS) {
         momyrhs(p.I) = alp_avg * sqrt_detg * mom_source(1);
         momzrhs(p.I) = alp_avg * sqrt_detg * mom_source(2);
         taurhs(p.I) = alp_avg * sqrt_detg * tau_source;
+
       }); // end of loop over grid
 }