Merge branch 'development' into ppm_one_sided_second

CHANGES.md

@@ -1,3 +1,41 @@
+# 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)
+
+  * update the `hse_convergence` readme to reflect current convergence
+    (#2946)
+
+  * update the `bubble_convergence` plotting script (#2947)
+
+  * new Frontier scaling numbers (#2948)
+
+  * more GPU error printing (#2944)
+
+  * science problem updates: `flame_wave` (#2943)
+
+  * documentation updates (#2939)
+
 # 24.08
 
   * lazy QueueReduction has been enabled for the timing diagnostics

Docs/source/mhd.rst

@@ -17,6 +17,9 @@ supported by MHD.
 
 .. note::
 
+   The MHD solver is still under development and should not be used
+   for science simulations.  See the issues in the github repo.
+
    The MHD solver supports 3-d only.
 
    Currently the MHD solver is single-level only.  AMR support is forthcoming.

Exec/gravity_tests/hse_convergence/README.md

@@ -7,29 +7,34 @@ in the plotfiles.
 
 To run this problem, use one of the convergence scripts:
 
-  * ``convergence_plm.sh`` :
+  * `convergence_plm.sh` :
 
-    this runs CTU + PLM using the default HSE BCs and default
-    use_pslope, then with reflect BCs, then without use_pslope, and
-    finally runs with reflect instead of HSE BCs.
+    this runs CTU + PLM using:
+    1. the default HSE BCs and `use_pslope`
+    2. the HSE BCs with reflection and `use_pslope`
+    3. reflect BCs instead of HSE BCs without `use_pslope`
+    4. reflect BCs with `use_pslope`
 
-    These tests show that the best results come from HSE BCs + reflect vel
+    These tests show that the best results (by far) come from
+    `use_pslope=1` and reflecting BCs
 
   * convergence_ppm.sh :
 
     this runs CTU + PPM in a similar set of configurations as PLM above
-    (with one additional one: grav_source_type = 4)
+    1. the default HSE BCs
+    2. HSE BCs with reflection
+    3. reflecting BCs
+    4. reflecting BCs with `use_pslope`
 
-    These tests show that the best results come from HSE BCs + reflect vel
+    These tests show that the best results (by far) come from
+    reflecting BCs with `use_pslope=1`, just like the PLM case.
 
   * convergence_sdc.sh :
 
-    this uses the TRUE_SDC integration, first with SDC-2 + PLM  and reflecting BCs,
-    the SDC-2 + PPM and reflecting BCs, then the same but HSE BCs, and finally
-    SDC-4 + reflect
+    this uses the TRUE_SDC integration, first with SDC-2 + PLM and
+    reflecting BCs, the SDC-2 + PPM and reflecting BCs, then the same
+    but HSE BCs, and finally SDC-4 + reflect
 
     These tests show that the PLM + reflect (which uses the
     well-balanced use_pslope) and the SDC-4 + reflect give the lowest
-    errors and expected (or better) convergence:
-
-
+    errors and expected (or better) convergence.

Exec/gravity_tests/hse_convergence/convergence_plm.sh

@@ -58,43 +58,15 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-## plm + hse reflect + no pslope
-
-ofile=plm-hsereflect-nopslope.converge.out
-
-RUNPARAMS="
-castro.ppm_type=0
-castro.use_pslope=0
-castro.hse_interp_temp=1
-castro.hse_reflect_vels=1
-"""
-
-${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out
-pfile=`ls -t | grep -i hse_64_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel > ${ofile}
-
-${EXEC} inputs.ppm.128 ${RUNPARAMS} >& 128.out
-pfile=`ls -t | grep -i hse_128_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
-
-${EXEC} inputs.ppm.256 ${RUNPARAMS} >& 256.out
-pfile=`ls -t | grep -i hse_256_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
-
-${EXEC} inputs.ppm.512 ${RUNPARAMS} >& 512.out
-pfile=`ls -t | grep -i hse_512_plt | head -1`
-fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
-
-
-## plm + reflect
+## plm + reflect + nopslope
 
-ofile=plm-reflect.converge.out
+ofile=plm-reflect-nopslope.converge.out
 
 RUNPARAMS="
 castro.ppm_type=0
-castro.use_pslope=1
 castro.lo_bc=3
 castro.hi_bc=3
+castro.use_pslope=0
 """
 
 ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out
@@ -114,16 +86,15 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
+## plm + reflect + pslope
 
-## plm + reflect + nopslope
-
-ofile=plm-reflect-nopslope.converge.out
+ofile=plm-reflect-pslope.converge.out
 
 RUNPARAMS="
 castro.ppm_type=0
 castro.lo_bc=3
 castro.hi_bc=3
-castro.use_pslope=0
+castro.use_pslope=1
 """
 
 ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out

Exec/gravity_tests/hse_convergence/convergence_ppm.sh

@@ -50,12 +50,13 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-## ppm + grav_source_type = 4
+## ppm + reflect
 
-ofile=ppm-grav4.converge.out
+ofile=ppm-reflect.converge.out
 
 RUNPARAMS="
-castro.grav_source_type=4
+castro.lo_bc=3
+castro.hi_bc=3
 """
 
 ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out
@@ -75,13 +76,14 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-## ppm + reflect
+## ppm + reflect + pslope
 
-ofile=ppm-reflect.converge.out
+ofile=ppm-reflect-pslope.converge.out
 
 RUNPARAMS="
 castro.lo_bc=3
 castro.hi_bc=3
+castro.use_pslope=1
 """
 
 ${EXEC} inputs.ppm.64 ${RUNPARAMS} >& 64.out

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/(gamma_const - 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/radiation_tests/Rad2Tshock/problem_initialize.H

@@ -14,8 +14,8 @@ void problem_initialize ()
 
     eos_state.rho = problem::rho0;
     eos_state.T = problem::T0;
-    for (int n = 0; n < NumSpec; n++) {
-        eos_state.xn[n] = 0.0_rt;
+    for (auto & X : eos_state.xn) {
+        X = 0.0_rt;
     }
     eos_state.xn[0]  = 1.0_rt;
 

Exec/radiation_tests/Rad2Tshock/problem_initialize_rad_data.H

@@ -12,6 +12,9 @@ void problem_initialize_rad_data (int i, int j, int k,
                                   const GeometryData& geomdata)
 {
 
+    amrex::ignore_unused(nugroup);
+    amrex::ignore_unused(dnugroup);
+
     const Real* dx = geomdata.CellSize();
     const Real* problo = geomdata.ProbLo();
 

Exec/radiation_tests/Rad2Tshock/problem_initialize_state_data.H

@@ -16,7 +16,7 @@ void problem_initialize_state_data (int i, int j, int k,
     // Provides the simulation to be run in the x,y,or z direction
     // where length direction is the length side in a square prism
 
-    Real length_cell;
+    Real length_cell{};
     if (problem::idir == 1) {
         length_cell = problo[0] + dx[0] * (static_cast<Real>(i) + 0.5_rt);
     } else if (problem::idir == 2) {