Merge branch 'main' into amr

AsterX/interface.ccl

@@ -104,3 +104,35 @@ CCTK_REAL analysis_quantities TYPE=gf CENTERING={ccc}
   b2small
   total_energy
 } "Analysis quantities"
+
+#grid functions required for upwindCT computation of electric field
+
+CCTK_REAL vtilde_xface TYPE=gf CENTERING={vcc} TAGS='checkpoint="no"'
+{
+  vtilde_y_xface, vtilde_z_xface
+} "staggered vtilde components on x-face where vtilde^i = \alpha*v^i-\Beta^i"
+
+CCTK_REAL vtilde_yface TYPE=gf CENTERING={cvc} TAGS='checkpoint="no"'
+{
+  vtilde_x_yface, vtilde_z_yface
+} "staggered vtilde components on y-face where vtilde^i = \alpha*v^i-\Beta^i"
+
+CCTK_REAL vtilde_zface TYPE=gf CENTERING={ccv} TAGS='checkpoint="no"'
+{
+  vtilde_x_zface, vtilde_y_zface
+} "staggered vtilde components on z-face where vtilde^i = \alpha*v^i-\Beta^i"
+
+CCTK_REAL a_xface TYPE=gf CENTERING={vcc} TAGS='checkpoint="no"'
+{
+  amax_xface, amin_xface
+} "staggered vtilde components on x-face where vtilde^i = \alpha*v^i-\Beta^i"
+
+CCTK_REAL a_yface TYPE=gf CENTERING={cvc} TAGS='checkpoint="no"'
+{
+  amax_yface, amin_yface
+} "staggered vtilde components on y-face where vtilde^i = \alpha*v^i-\Beta^i"
+
+CCTK_REAL a_zface TYPE=gf CENTERING={ccv} TAGS='checkpoint="no"'
+{
+  amax_zface, amin_zface
+} "staggered vtilde components on z-face where vtilde^i = \alpha*v^i-\Beta^i"

AsterX/param.ccl

@@ -11,6 +11,10 @@ CCTK_REAL lorenz_damp_fac "Damping factor used in the Lorenz Gauge (see Farris e
  0:* :: "should be positive"
 } 0.0
 
+BOOLEAN use_uct "Shall we use the Upwind-CT method to compute the electric field (instead of the fluxCT approach)?" STEERABLE=always
+{
+} "no"
+
 KEYWORD flux_type "Flux solver" STEERABLE=always
 {
   "LxF" :: ""
@@ -102,6 +106,10 @@ USES CCTK_REAL rho_strict
 USES BOOLEAN Ye_lenient
 USES CCTK_INT max_iter
 USES CCTK_REAL c2p_tol
+USES CCTK_REAL r_atmo
+USES CCTK_REAL n_rho_atmo
+USES CCTK_REAL n_press_atmo
+USES BOOLEAN thermal_eos_atmo
 
 SHARES: ReconX
 USES KEYWORD reconstruction_method

AsterX/schedule.ccl

@@ -167,11 +167,15 @@ SCHEDULE AsterX_Fluxes IN AsterX_RHSGroup
   READS: dens(everywhere) tau(everywhere) mom(everywhere)
   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: 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)
+  WRITES: vtilde_xface(interior) vtilde_yface(interior) vtilde_zface(interior)
+  WRITES: a_xface(interior) a_yface(interior) a_zface(interior)
   SYNC: Aux_in_RHSof_A_Psi
   SYNC: flux_x flux_y flux_z
+  SYNC: vtilde_xface vtilde_yface vtilde_zface a_xface a_yface a_zface
 } "Calculate the hydro fluxes"
 
 SCHEDULE AsterX_SourceTerms IN AsterX_RHSGroup AFTER AsterX_Fluxes
@@ -191,10 +195,14 @@ SCHEDULE AsterX_RHS IN AsterX_RHSGroup AFTER AsterX_SourceTerms
 {
   LANG: C
   READS: ADMBaseX::metric(everywhere) ADMBaseX::lapse(everywhere) ADMBaseX::shift(everywhere)
+  READS: HydroBaseX::vel(everywhere) HydroBaseX::press(everywhere)
   READS: flux_x(everywhere) flux_y(everywhere) flux_z(everywhere)
   READS: densrhs(everywhere) taurhs(everywhere) momrhs(everywhere)
   READS: Psi(everywhere)
   READS: Aux_in_RHSof_A_Psi(everywhere)
+  READS: dBx_stag(everywhere) dBy_stag(everywhere) dBz_stag(everywhere)
+  READS: vtilde_xface(everywhere) vtilde_yface(everywhere) vtilde_zface(everywhere)
+  READS: a_xface(everywhere) a_yface(everywhere) a_zface(everywhere)
   WRITES: densrhs(interior) taurhs(interior) momrhs(interior)
   WRITES: Avec_x_rhs(interior) Avec_y_rhs(interior) Avec_z_rhs(interior) Psi_rhs(interior)
   SYNC: densrhs taurhs momrhs

AsterX/src/con2prim.cxx

@@ -55,28 +55,45 @@ void AsterX_Con2Prim_typeEoS(CCTK_ARGUMENTS, EOSIDType &eos_cold,
 
   const smat<GF3D2<const CCTK_REAL>, 3> gf_g{gxx, gxy, gxz, gyy, gyz, gzz};
 
-  // Setting up atmosphere
-  const CCTK_REAL rho_atmo_cut = rho_abs_min * (1 + atmo_tol);
-  const CCTK_REAL gm1 = eos_cold.gm1_from_valid_rmd(rho_abs_min);
-  CCTK_REAL eps_atm = eos_cold.sed_from_valid_gm1(gm1);
-  eps_atm = std::min(std::max(eos_th.rgeps.min, eps_atm), eos_th.rgeps.max);
-  const CCTK_REAL p_atm =
-      eos_th.press_from_valid_rho_eps_ye(rho_abs_min, eps_atm, Ye_atmo);
-  atmosphere atmo(rho_abs_min, eps_atm, Ye_atmo, p_atm, rho_atmo_cut);
-
-  // Construct Noble c2p object:
-  c2p_2DNoble c2p_Noble(eos_th, atmo, max_iter, c2p_tol, rho_strict, vw_lim,
-                        B_lim, Ye_lenient);
-
-  // Construct Palenzuela c2p object:
-  c2p_1DPalenzuela c2p_Pal(eos_th, atmo, max_iter, c2p_tol, rho_strict, vw_lim,
-                           B_lim, Ye_lenient);
-
   // Loop over the interior of the grid
   cctk_grid.loop_int_device<
       1, 1, 1>(grid.nghostzones, [=] CCTK_DEVICE(
                                      const PointDesc
                                          &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+    // Setting up atmosphere
+    CCTK_REAL rho_atm = 0.0;   // dummy initialization
+    CCTK_REAL press_atm = 0.0; // dummy initialization
+    CCTK_REAL eps_atm = 0.0;   // dummy initialization
+    CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
+
+    // Grading rho
+    rho_atm = (radial_distance > r_atmo)
+                  ? (rho_abs_min * pow((r_atmo / radial_distance), n_rho_atmo))
+                  : rho_abs_min;
+    const CCTK_REAL rho_atmo_cut = rho_atm * (1 + atmo_tol);
+
+    // Grading pressure based on either cold or thermal EOS
+    if (thermal_eos_atmo) {
+      press_atm = (radial_distance > r_atmo)
+                      ? (p_atmo * pow(r_atmo / radial_distance, n_press_atmo))
+                      : p_atmo;
+      eps_atm = eos_th.eps_from_valid_rho_press_ye(rho_atm, press_atm, Ye_atmo);
+    } else {
+      const CCTK_REAL gm1 = eos_cold.gm1_from_valid_rmd(rho_atm);
+      eps_atm = eos_cold.sed_from_valid_gm1(gm1);
+      eps_atm = std::min(std::max(eos_th.rgeps.min, eps_atm), eos_th.rgeps.max);
+      press_atm = eos_th.press_from_valid_rho_eps_ye(rho_atm, eps_atm, Ye_atmo);
+    }
+    atmosphere atmo(rho_atm, eps_atm, Ye_atmo, press_atm, rho_atmo_cut);
+
+    // Construct Noble c2p object:
+    c2p_2DNoble c2p_Noble(eos_th, atmo, max_iter, c2p_tol, rho_strict, vw_lim,
+                          B_lim, Ye_lenient);
+
+    // Construct Palenzuela c2p object:
+    c2p_1DPalenzuela c2p_Pal(eos_th, atmo, max_iter, c2p_tol, rho_strict,
+                             vw_lim, B_lim, Ye_lenient);
+
     /* Get covariant metric */
     const smat<CCTK_REAL, 3> glo(
         [&](int i, int j) ARITH_INLINE { return calc_avg_v2c(gf_g(i, j), p); });

AsterX/src/con2prim_rpa.cxx

@@ -25,37 +25,57 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
 
   const smat<GF3D2<const CCTK_REAL>, 3> gf_g{gxx, gxy, gxz, gyy, gyz, gzz};
 
-  // Setting up initial data EOS
-  const CCTK_REAL n = 1 / (poly_gamma - 1); // Polytropic index
-  const CCTK_REAL adiab_ind_id = 1.0 / (poly_gamma - 1);
-  const CCTK_REAL rmd_p = pow(poly_k, -n); //Polytropic density scale
-  const auto eos_id = make_eos_barotr_poly(adiab_ind_id, rmd_p, rho_max);
-
-  // Setting up evolution EOS
-  const CCTK_REAL adiab_ind_evol = 1.0 / (gl_gamma - 1);
-  const auto eos = make_eos_idealgas(adiab_ind_evol, eps_max, rho_max);
-
-  // Setting up atmosphere
-  const CCTK_REAL rho_atmo_cut = rho_abs_min * (1 + atmo_tol);
-  CCTK_REAL eps_atm = eos_id.at_rho(rho_abs_min).eps();
-  eps_atm  = eos.range_eps(rho_abs_min, Ye_atmo).limit_to(eps_atm);
-  CCTK_REAL p_atm  = eos.at_rho_eps_ye(rho_abs_min, eps_atm, Ye_atmo).press();
-  const atmosphere atmo(rho_abs_min, eps_atm, Ye_atmo, p_atm, rho_atmo_cut);
-
-  CCTK_REAL dummy_Ye = 0.5;
-  CCTK_REAL dummy_dYe = 0.5;
-
-  // Get a recovery function
-  con2prim_mhd cv2pv(eos, 1e-5, 1, 100, 100, atmo, 1e-8, 100);
-  //  con2prim_mhd cv2pv(eos, rho_strict, ye_lenient, max_z, max_b, atmo,
-  //  c2p_acc,
-  //                  max_iter);
-
   // Loop over the interior of the grid
   cctk_grid.loop_int_device<
       1, 1, 1>(grid.nghostzones, [=] CCTK_DEVICE(
                                      const PointDesc
                                          &p) CCTK_ATTRIBUTE_ALWAYS_INLINE {
+    // Setting up initial data EOS
+    const CCTK_REAL n = 1 / (poly_gamma - 1); // Polytropic index
+    const CCTK_REAL adiab_ind_id = 1.0 / (poly_gamma - 1);
+    const CCTK_REAL rmd_p = pow(poly_k, -n); // Polytropic density scale
+    const auto eos_id = make_eos_barotr_poly(adiab_ind_id, rmd_p, rho_max);
+
+    // Setting up evolution EOS
+    const CCTK_REAL adiab_ind_evol = 1.0 / (gl_gamma - 1);
+    const auto eos = make_eos_idealgas(adiab_ind_evol, eps_max, rho_max);
+
+    // Setting up atmosphere
+
+    CCTK_REAL rho_atm = 0.0;   // dummy initialization
+    CCTK_REAL press_atm = 0.0; // dummy initialization
+    CCTK_REAL eps_atm = 0.0;   // dummy initialization
+    CCTK_REAL radial_distance = sqrt(p.x * p.x + p.y * p.y + p.z * p.z);
+
+    // Grading rho
+    rho_atm = (radial_distance > r_atmo)
+                  ? (rho_abs_min * pow((r_atmo / radial_distance), n_rho_atmo))
+                  : rho_abs_min;
+    const CCTK_REAL rho_atmo_cut = rho_atm * (1 + atmo_tol);
+
+    // Grading pressure and eps
+    if (thermal_eos_atmo) {
+      press_atm = (radial_distance > r_atmo)
+                      ? (p_atmo * pow(r_atmo / radial_distance, n_press_atmo))
+                      : p_atmo;
+      //TODO: eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps() does not exist in RePrimAnd
+      //Currently computing eps from ideal gas EOS
+      //eps_atm = eos.at_rho_press_ye(rho_atm, press_atm, Ye_atmo).eps();
+      eps_atm = press_atm/(rho_atm*(gl_gamma - 1.));
+    } else {
+      eps_atm = eos_id.at_rho(rho_atm).eps();
+      eps_atm = eos.range_eps(rho_atm, Ye_atmo).limit_to(eps_atm);
+      press_atm = eos.at_rho_eps_ye(rho_atm, eps_atm, Ye_atmo).press();
+    }
+    const atmosphere atmo(rho_atm, eps_atm, Ye_atmo, press_atm, rho_atmo_cut);
+
+    CCTK_REAL dummy_Ye = 0.5;
+    CCTK_REAL dummy_dYe = 0.5;
+
+    // Get a recovery function
+    con2prim_mhd cv2pv(eos, rho_strict, Ye_lenient, vw_lim, B_lim, atmo, c2p_tol,
+                       max_iter);
+
     /* Get covariant metric */
     const smat<CCTK_REAL, 3> glo(
         [&](int i, int j) ARITH_INLINE { return calc_avg_v2c(gf_g(i, j), p); });
@@ -117,7 +137,7 @@ extern "C" void AsterX_Con2Prim(CCTK_ARGUMENTS) {
             // velz(p.I), Bvecx(p.I), Bvecy(p.I), Bvecz(p.I),
             Avec_x(p.I), Avec_y(p.I), Avec_z(p.I));
       }
-      //set to atmo
+      // set to atmo
       cv.bcons(0) = dBx(p.I);
       cv.bcons(1) = dBy(p.I);
       cv.bcons(2) = dBz(p.I);

AsterX/src/fluxes.cxx

@@ -46,8 +46,16 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
   /* grid functions */
   const vec<GF3D2<const CCTK_REAL>, dim> gf_vels{velx, vely, velz};
   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};
   const smat<GF3D2<const CCTK_REAL>, dim> gf_g{gxx, gxy, gxz, gyy, gyz, gzz};
+  /* grid functions for Upwind CT */
+  const vec<GF3D2<CCTK_REAL>, dim> vtildes_one{vtilde_y_xface, vtilde_z_yface,
+                                               vtilde_x_zface};
+  const vec<GF3D2<CCTK_REAL>, dim> vtildes_two{vtilde_z_xface, vtilde_x_yface,
+                                               vtilde_y_zface};
+  const vec<GF3D2<CCTK_REAL>, dim> amax{amax_xface, amax_yface, amax_zface};
+  const vec<GF3D2<CCTK_REAL>, dim> amin{amin_xface, amin_yface, amin_zface};
 
   static_assert(dir >= 0 && dir < 3, "");
 
@@ -156,6 +164,10 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
   constexpr array<int, dim> face_centred = {!(dir == 0), !(dir == 1),
                                             !(dir == 2)};
 
+  constexpr array<int, dim> dir_arr = {(dir==0) ? 2 : ( (dir==1) ? 0 : 1 ), 
+                                       dir,
+                                       (dir==0) ? 1 : ( (dir==1) ? 2 : 0 )};
+
   grid.loop_int_device<
       face_centred[0], face_centred[1],
       face_centred
@@ -194,10 +206,6 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
           eos_th.press_from_valid_rho_eps_ye(rho_rc(1), eps_min, ye_rc(1));
     }
 
-    const vec<vec<CCTK_REAL, 2>, 3> Bs_rc([&](int i) ARITH_INLINE {
-      return vec<CCTK_REAL, 2>{reconstruct_pt(gf_Bvecs(i), p, false, false)};
-    });
-
     /* 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(
@@ -222,6 +230,26 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
       return 1 / sqrt(1 - calc_contraction(vlows_rc, vels_rc)(f));
     });
 
+    // Introduce reconstructed Bs
+    // Use staggered dB for i == dir
+
+    vec<vec<CCTK_REAL, 2>, 3> Bs_rc;
+    array<CCTK_REAL,2> Bs_rc_dummy; // note: can't copy array<,2> to vec<,2>, only construct
+
+    Bs_rc(dir)(0) = gf_dBstags(dir)(p.I)/sqrtg;
+    Bs_rc(dir)(1) = Bs_rc(dir)(0);
+
+    Bs_rc_dummy = reconstruct_pt(gf_Bvecs(dir_arr[0]), p, false, false);
+    Bs_rc(dir_arr[0])(0) = Bs_rc_dummy[0];
+    Bs_rc(dir_arr[0])(1) = Bs_rc_dummy[1];
+
+    Bs_rc_dummy = reconstruct_pt(gf_Bvecs(dir_arr[2]), p, false, false);
+    Bs_rc(dir_arr[2])(0) = Bs_rc_dummy[0];
+    Bs_rc(dir_arr[2])(1) = Bs_rc_dummy[1];
+
+    // End of setting Bs
+
+
     /* alpha * b0 = W * B^i * v_i */
     const vec<CCTK_REAL, 2> alp_b0_rc([&](int f) ARITH_INLINE {
       return w_lorentz_rc(f) * calc_contraction(Bs_rc, vlows_rc)(f);
@@ -449,6 +477,31 @@ void CalcFlux(CCTK_ARGUMENTS, EOSType &eos_th) {
       printf("  vtilde_rc = %16.8e, %16.8e.\n", vtilde_rc(0), vtilde_rc(1));
       assert(0);
     }
+
+    /* Begin code for upwindCT */
+        // if dir==0: dir1=1, dir2=2 | dir==1: dir1=2, dir2=0 | dir==2; dir1=0,
+        // dir2=1
+
+        const int dir1 = (dir == 0) ? 1 : ((dir == 1) ? 2 : 0);
+        const int dir2 = (dir == 0) ? 2 : ((dir == 1) ? 0 : 1);
+
+        amax(dir)(p.I) = max({CCTK_REAL(0), lambda(0)(0), lambda(0)(1),
+                              lambda(0)(2), lambda(0)(3), lambda(1)(0),
+                              lambda(1)(1), lambda(1)(2), lambda(1)(3)});
+
+        amin(dir)(p.I) = -1 * (min({CCTK_REAL(0), lambda(0)(0), lambda(0)(1),
+                                    lambda(0)(2), lambda(0)(3), lambda(1)(0),
+                                    lambda(1)(1), lambda(1)(2), lambda(1)(3)}));
+
+        vtildes_one(dir)(p.I) = (amax(dir)(p.I) * vtildes_rc(dir1)(0) +
+                                 amin(dir)(p.I) * vtildes_rc(dir1)(1)) /
+                                (amax(dir)(p.I) + amin(dir)(p.I));
+        vtildes_two(dir)(p.I) = (amax(dir)(p.I) * vtildes_rc(dir2)(0) +
+                                 amin(dir)(p.I) * vtildes_rc(dir2)(1)) /
+                                (amax(dir)(p.I) + amin(dir)(p.I));
+
+        /* End code for upwindCT */
+
   });
 }