Merge branch 'development' into ppm_full_well_balancing

CHANGES.md

@@ -1,3 +1,24 @@
+# 24.10
+
+  * update initial model for `subchandra` when doing ASE NSE (#2970)
+
+  * `massive_star` plot tweaks (#2968)
+
+  * start of work on 2D spherical geometry (#2953, #2954, #2955,
+    #2957, #2958, #2959, #2961, #2962, #2964, #2965)
+
+  * the gresho problem now takes Mach number instead of p0 as input
+    (#2951m #2963)
+
+  * the PPM geometric source terms in the normal predictor are now
+    traced to the interfaces (#2473)
+
+  * `subch_planar` now works in 1D (#2952)
+
+  * remove old `get_const_grav()` function (#2956)
+
+  * clang-tidy fixes to radiation (#2950)
+
 # 24.09
 
   * Code clean-ups / clang-tidy (#2942, #2949)
@@ -9,7 +30,7 @@
 
   * new Frontier scaling numbers (#2948)
 
-  * more GPU error printing (@3944)
+  * more GPU error printing (#2944)
 
   * science problem updates: `flame_wave` (#2943)
 

Exec/hydro_tests/gresho_vortex/_prob_params

@@ -1,11 +1,15 @@
-p0      real     1.0_rt       y
+# initial Mach number
+M0      real     0.1_rt       y
 
+# ambient density
 rho0    real     1.0_rt       y
 
+# reference pressure -- this is set automatically
+p0      real     -1.0_rt
+
+# reference timescale
 t_r     real     1.0_rt       y
 
 x_r     real     0.0_rt
 
 q_r     real     0.0_rt
-
-nsub    integer   4           y

Exec/hydro_tests/gresho_vortex/inputs.2d

@@ -79,16 +79,8 @@ amr.derive_plot_vars=ALL
 # PROBLEM PARAMETERS
 problem.rho0 = 1.0
 
-# p0 = rho0 u_phi**2/(gamma * M**2) - 1/2  -- specify M to get p0
-
-# M ~ 0.001
-problem.p0 = 599999.5
-
-# M ~ 0.01
-# problem.p0 = 5999.5
-
-# M ~ 0.1
-# problem.p0 = 59.5
+# Mach number of the flow
+problem.M0 = 0.1
 
 # EOS
 eos.eos_assume_neutral = 1

Exec/hydro_tests/gresho_vortex/problem_initialize.H

@@ -27,6 +27,17 @@ void problem_initialize ()
     problem::x_r = probhi[0] - problo[0];
     problem::q_r = 0.4_rt * M_PI * problem::x_r / problem::t_r;
 
+    // pressure peaks at r = 1/5 from the center
+    // where p = p0 + 25/2 r^2, so peak pressure is p = p0 + 1/2
+    //
+    // Mach number is |u_phi| / sqrt(gamma p / rho), so we get solve
+    // for p0.  Note that the peak velocity is q,
+    //
+    // p0 = rho q^2 / (gamma M^2) - 1/2
+
+    problem::p0 = problem::rho0 * std::pow(problem::q_r, 2.0) /
+        (eos_rp::eos_gamma * std::pow(problem::M0, 2.0)) - 0.5_rt;
+
 }
 
 #endif

Exec/hydro_tests/gresho_vortex/problem_initialize_state_data.H

@@ -22,71 +22,46 @@ void problem_initialize_state_data (int i, int j, int k,
     const Real* dx = geomdata.CellSize();
 
     Real x = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
-    Real xl = problo[0] + dx[0] * (static_cast<Real>(i));
 
     Real y = 0.0;
-    Real yl = 0.0;
-
 #if AMREX_SPACEDIM >= 2
     y = problo[1] + dx[1] * (static_cast<Real>(j) + 0.5_rt);
-    yl = problo[1] + dx[1] * (static_cast<Real>(j));
 #endif
 
     Real z = 0.0;
-    Real zl = 0.0;
-
 #if AMREX_SPACEDIM == 3
     z = problo[2] + dx[2] * (static_cast<Real>(k) + 0.5_rt);
-    zl = problo[2] + dx[2] * (static_cast<Real>(k));
 #endif
 
     Real reint = 0.0;
     Real u_tot = 0.0;
 
-    for (int kk = 0; kk < problem::nsub; kk++) {
-        Real zz = zl + dx[2] * (static_cast<Real>(kk) + 0.5_rt) / problem::nsub;
-
-        for (int jj = 0; jj < problem::nsub; jj++) {
-            Real yy = yl + dx[1] * (static_cast<Real>(jj) + 0.5_rt) / problem::nsub;
-
-            for (int ii = 0; ii < problem::nsub; ii++) {
-                Real xx = xl + dx[0] * (static_cast<Real>(ii) + 0.5_rt) / problem::nsub;
-
-                Real r = std::sqrt((xx - problem::center[0]) * (xx - problem::center[0]) +
-                                   (yy - problem::center[1]) * (yy - problem::center[1]));
-
-                Real u_phi;
-                Real p;
+    Real r = std::sqrt(amrex::Math::powi<2>(x - problem::center[0]) +
+                       amrex::Math::powi<2>(y - problem::center[1]));
 
-                if (r < 0.2_rt) {
-                    u_phi = 5.0_rt * r;
-                    p = problem::p0 + 12.5_rt * r * r;
+    Real u_phi;
+    Real p;
 
-                } else if (r < 0.4_rt) {
-                    u_phi = 2.0_rt - 5.0_rt * r;
-                    p = problem::p0 + 12.5_rt * r * r + 4.0_rt *
-                        (1.0_rt - 5.0_rt * r - std::log(0.2_rt) + std::log(r));
+    if (r < 0.2_rt) {
+        u_phi = 5.0_rt * r;
+        p = problem::p0 + 12.5_rt * r * r;
 
-                } else {
-                    u_phi = 0.0_rt;
-                    p = problem::p0 - 2.0_rt + 4.0_rt * std::log(2.0_rt);
-                }
+    } else if (r < 0.4_rt) {
+        u_phi = 2.0_rt - 5.0_rt * r;
+        p = problem::p0 + 12.5_rt * r * r + 4.0_rt *
+            (1.0_rt - 5.0_rt * r - std::log(0.2_rt) + std::log(r));
 
-                u_tot += u_phi;
-                reint += p/(eos_rp::eos_gamma - 1.0_rt);
-            }
-        }
+    } else {
+        u_phi = 0.0_rt;
+        p = problem::p0 - 2.0_rt + 4.0_rt * std::log(2.0_rt);
     }
 
-    Real u_phi = u_tot / (problem::nsub * problem::nsub * problem::nsub);
-    reint = reint / (problem::nsub * problem::nsub * problem::nsub);
+    reint = p / (eos_rp::eos_gamma - 1.0_rt);
 
     state(i,j,k,URHO) = problem::rho0;
 
     // phi unit vector: \hat{\phi} = -sin(phi) \hat{x} + cos(phi) \hat{y}
     // with cos(phi) = x/r; sin(phi) = y/r
-    Real r = std::sqrt((x - problem::center[0]) * (x - problem::center[0]) +
-                       (y - problem::center[1]) * (y - problem::center[1]));
 
     //  -sin(phi) = y/r
     state(i,j,k,UMX) = -problem::rho0 * problem::q_r * u_phi * ((y - problem::center[1]) / r);

Exec/science/massive_star/analysis/initial_model.py

@@ -45,8 +45,8 @@ def find_r_for_rho(r, rho, rho_want):
 fig = plt.figure()
 ax = fig.add_subplot(211)
 
-l1 = ax.plot(data[:,0], data[:,idens], label="density")
-l2 = ax.plot(data[:,0], data[:,itemp], label="temperature")
+l1 = ax.plot(data[:,0], data[:,idens], label=r"$\rho$")
+l2 = ax.plot(data[:,0], data[:,itemp], label="$T$")
 
 # show where the refinement kicks in
 rho_refine = 2.e4
@@ -72,12 +72,12 @@ def find_r_for_rho(r, rho, rho_want):
 
 ax2.set_ylabel(r"$Y_e$")
 
-l3 = ax2.plot(data[:,0], data[:,iye], color="C2", label="Ye")
+l3 = ax2.plot(data[:,0], data[:,iye], color="C2", label="$Y_e$")
 
 lns = l1 + l2 + l3
 labs = [l.get_label() for l in lns]
 
-ax.legend(lns, labs, frameon=False, loc=6)
+ax.legend(lns, labs, frameon=False, loc=6, fontsize="small")
 
 # species plot
 
@@ -119,11 +119,10 @@ def find_r_for_rho(r, rho, rho_want):
 
 ax.grid(color="b", alpha=0.5, ls=":")
 
-ax.legend(frameon=True, edgecolor="w", ncol=1, framealpha=0.5)
+ax.legend(frameon=True, edgecolor="w", ncol=1, framealpha=0.5, fontsize="small")
 
-fig.set_size_inches((8, 12))
+fig.set_size_inches((6, 9))
 
 fig.tight_layout()
 
 fig.savefig("initial_model.pdf")
-

Exec/science/subchandra/GNUmakefile.nse_net

@@ -7,22 +7,23 @@ DEBUG      = FALSE
 
 DIM        = 2
 
-COMP	   = gnu
+COMP       = gnu
 
 USE_MPI    = TRUE
 
 USE_GRAV   = TRUE
 USE_REACT  = TRUE
 
 USE_NSE_NET = TRUE
+USE_SIMPLIFIED_SDC = TRUE
 
 USE_MODEL_PARSER = TRUE
 
 # This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
 EOS_DIR     := helmholtz
 
 # This sets the network directory in $(MICROPHYSICS_HOME)/networks
-NETWORK_DIR := ase
+NETWORK_DIR := subch_base
 SCREEN_METHOD = chabrier1998
 
 INTEGRATOR_DIR := VODE