Merge branch '2d_spherical_test' of github.com:zhichen3/Castro into 2…

…d_spherical_test
AMReX-Astro · Nov 13, 2024 · c1723e1 · c1723e1
2 parents 96a07fa + 4f53b47
commit c1723e1
Show file tree

Hide file tree

Showing 16 changed files with 268 additions and 131 deletions.
diff --git a/Docs/source/hse.rst b/Docs/source/hse.rst
@@ -74,6 +74,12 @@ then done on the parabola, again working with :math:`\tilde{p}`.
 Finally, the parabola values are updated to include the hydrostatic
 pressure.
 
+.. index:: castro.ppm_well_balanced
+
+We can do better with PPM, and only use the perturbational pressure,
+$\tilde{p}$, in the characteristic tracing and then add back the
+hydrostatic pressure to the interface afterwards.  This is done via
+``castro.ppm_well_balanced=1``.
 
 Fully fourth-order method
 -------------------------
@@ -165,5 +171,3 @@ test the different HSE approaches.  This sets up a 1-d X-ray burst
 atmosphere (based on the ``flame_wave`` setup).  Richardson
 extrapolation can be used to measure the convergence rate (or just
 look at how the peak velocity changes).
-
-
diff --git a/Exec/gravity_tests/hse_convergence/README.md b/Exec/gravity_tests/hse_convergence/README.md
@@ -31,10 +31,13 @@ To run this problem, use one of the convergence scripts:
 
   * 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, with the following variations:
+    1. SDC-2 + PLM and reflecting BCs
+    2. SDC-2 + PPM and reflecting BCs
+    3. SDC-2 + PLM with HSE BCs
+    4. SDC-2 + PPM with HSE BCs
+    5. 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
+    well-balanced `use_pslope`) and the SDC-4 + reflect give the lowest
     errors and expected (or better) convergence.
diff --git a/Exec/gravity_tests/hse_convergence/convergence_ppm.sh b/Exec/gravity_tests/hse_convergence/convergence_ppm.sh
@@ -102,3 +102,30 @@ ${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}
 
+
+## ppm + reflect + well-balanced
+
+ofile=ppm-reflect-wellbalanced.converge.out
+
+RUNPARAMS="
+castro.lo_bc=3
+castro.hi_bc=3
+castro.ppm_well_balanced=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}
+
diff --git a/Exec/gravity_tests/hse_convergence/convergence_sdc.sh b/Exec/gravity_tests/hse_convergence/convergence_sdc.sh
@@ -67,14 +67,15 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-## sdc-2 + ppm
+## sdc-2 + HSE
 
-ofile=sdc2-ppm.converge.out
+ofile=sdc2.converge.out
 
 RUNPARAMS="
 castro.time_integration_method=2
 castro.sdc_order=2
-castro.ppm_type=1
+castro.ppm_type=0
+castro.use_pslope=1
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 """
@@ -96,15 +97,14 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-## sdc-2 + HSE
+## sdc-2 + ppm
 
-ofile=sdc2.converge.out
+ofile=sdc2-ppm.converge.out
 
 RUNPARAMS="
 castro.time_integration_method=2
 castro.sdc_order=2
-castro.ppm_type=0
-castro.use_pslope=1
+castro.ppm_type=1
 castro.limit_fourth_order=1
 castro.use_reconstructed_gamma1=1
 """
@@ -126,7 +126,6 @@ pfile=`ls -t | grep -i hse_512_plt | head -1`
 fextrema.gnu.ex -v magvel ${pfile} | grep -i magvel >> ${ofile}
 
 
-
 ## sdc-4 + reflect
 
 ofile=sdc4-reflect.converge.out

diff --git a/Exec/science/flame_wave/ci-benchmarks/grid_diag.out b/Exec/science/flame_wave/ci-benchmarks/grid_diag.out
@@ -1,5 +1,5 @@
 #   COLUMN 1                        2                        3                        4                        5                        6                        7                        8                        9                       10                       11                       12                       13                       14                       15                       16                       17                       18                       19                       20                       21                       22                       23
 #   TIMESTEP                     TIME                     MASS                     XMOM                     YMOM                     ZMOM              ANG. MOM. X              ANG. MOM. Y              ANG. MOM. Z              KIN. ENERGY              INT. ENERGY               GAS ENERGY             GRAV. ENERGY             TOTAL ENERGY     CENTER OF MASS X-LOC     CENTER OF MASS Y-LOC     CENTER OF MASS Z-LOC     CENTER OF MASS X-VEL     CENTER OF MASS Y-VEL     CENTER OF MASS Z-VEL      MAXIMUM TEMPERATURE          MAXIMUM DENSITY        MAXIMUM T_S / T_E
            0   0.0000000000000000e+00   1.4700685690736437e+20   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   2.9163021935541997e+37   2.9163021935541997e+37   0.0000000000000000e+00   2.9163021935541997e+37   2.7306641645125415e+04   6.8087525965685256e+02   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   0.0000000000000000e+00   1.3928678114307070e+09   3.0715027784567930e+07   0.0000000000000000e+00
-           1   5.9806047036494849e-08   1.4700700110495903e+20   3.3564221217849212e+23   7.2876252590226133e+23   1.2612378556916341e+20   3.1318391966955516e+23  -2.7141686347841742e+24   1.9066097493127676e+28   7.8605324800324878e+29   2.9163081863790706e+37   2.9163082649843738e+37   0.0000000000000000e+00   2.9163082649843738e+37   2.7306641717925955e+04   6.8087482283157317e+02   0.0000000000000000e+00   2.2831716153358752e+03   4.9573321027203638e+03   8.5794407491595881e-01   1.3929561431310942e+09   3.0715326643212341e+07   3.3873763041849706e-04
-           2   1.2559269877663918e-07   1.4700717750115462e+20   7.0484895445039834e+23   1.5304335700823335e+24   5.5619061962888426e+20   1.3810900989697918e+24  -1.1969287484132617e+25   4.0039912979259806e+28   3.0853275723280981e+30   2.9163145904135894e+37   2.9163148989463411e+37   0.0000000000000000e+00   2.9163148989463411e+37   2.7306641954244522e+04   6.8087461718998338e+02   0.0000000000000000e+00   4.7946567401095926e+03   1.0410604407871944e+04   3.7834249257966728e+00   1.3933944939675827e+09   3.0715562157563787e+07   3.3857586854903492e-04
+           1   5.9806047036494849e-08   1.4700700110495903e+20   3.3564221215590610e+23   7.2876252589700228e+23   1.2612378556594561e+20   3.1318391966786556e+23  -2.7141686347838510e+24   1.9066097493139366e+28   7.8605324800584047e+29   2.9163081863790706e+37   2.9163082649843738e+37   0.0000000000000000e+00   2.9163082649843738e+37   2.7306641717925955e+04   6.8087482283158556e+02   0.0000000000000000e+00   2.2831716151822361e+03   4.9573321026845897e+03   8.5794407489407010e-01   1.3929529135508094e+09   3.0715326539046615e+07   3.3873763041849706e-04
+           2   1.2559269877663918e-07   1.4700717750115462e+20   7.0484895413073359e+23   1.5304335700221490e+24   5.5619061951241526e+20   1.3810900989051759e+24  -1.1969287484120623e+25   4.0039912979436722e+28   3.0853275724059276e+30   2.9163145904135894e+37   2.9163148989463416e+37   0.0000000000000000e+00   2.9163148989463416e+37   2.7306641954244522e+04   6.8087461719003727e+02   0.0000000000000000e+00   4.7946567379351090e+03   1.0410604407462544e+04   3.7834249250044056e+00   1.3933796339641171e+09   3.0715561813912578e+07   3.3857586854903492e-04
diff --git a/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt b/Exec/science/flame_wave/ci-benchmarks/job_info_params.txt
@@ -64,6 +64,7 @@
     castro.dual_energy_eta1 = 1
     castro.dual_energy_eta2 = 0.0001
     castro.use_pslope = 0
+    castro.ppm_well_balanced = 0
     castro.pslope_cutoff_density = -1e+20
     castro.limit_fluxes_on_small_dens = 0
     castro.speed_limit = 0

diff --git a/Exec/science/flame_wave/ci-benchmarks/species_diag.out b/Exec/science/flame_wave/ci-benchmarks/species_diag.out
@@ -1,5 +1,5 @@
 #   COLUMN 1                        2                        3                        4                        5                        6                        7                        8                        9                       10                       11                       12                       13                       14                       15
 #   TIMESTEP                     TIME                 Mass He4                 Mass C12                 Mass O16                Mass Ne20                Mass Mg24                Mass Si28                 Mass S32                Mass Ar36                Mass Ca40                Mass Ti44                Mass Cr48                Mass Fe52                Mass Ni56
            0   0.0000000000000000e+00   2.8142378112400895e-15   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.3932235330395157e-24   7.1117997526547418e-14
-           1   5.9806047036494849e-08   2.8142323856670446e-15   5.4329394398418374e-21   7.4194443455948811e-24   7.3934514829832856e-24   7.3935541949306949e-24   7.3933935989353088e-24   7.3932321468740865e-24   7.3932312019180345e-24   7.3932308328287333e-24   7.3932307869392985e-24   7.3932307850741519e-24   7.3932307849933044e-24   7.1118070045957103e-14
-           2   1.2559269877663918e-07   2.8142264165381050e-15   1.1401993450709496e-20   7.4929528148974785e-24   7.3944694290388374e-24   7.3939305042586412e-24   7.3935838894001731e-24   7.3932425308122493e-24   7.3932405326301655e-24   7.3932397568503098e-24   7.3932396603620498e-24   7.3932396564405566e-24   7.3932396562705492e-24   7.1118158758588142e-14
+           1   5.9806047036494849e-08   2.8142323856670446e-15   5.4329394398523593e-21   7.4194443428333686e-24   7.3934514828234477e-24   7.3935541949253214e-24   7.3933935989348298e-24   7.3932321468740835e-24   7.3932312019180345e-24   7.3932308328287333e-24   7.3932307869392985e-24   7.3932307850741519e-24   7.3932307849933044e-24   7.1118070045957103e-14
+           2   1.2559269877663918e-07   2.8142264165381050e-15   1.1401993450700379e-20   7.4929526498089992e-24   7.3944694218302123e-24   7.3939305039811099e-24   7.3935838893792728e-24   7.3932425308121303e-24   7.3932405326301640e-24   7.3932397568503098e-24   7.3932396603620498e-24   7.3932396564405580e-24   7.3932396562705492e-24   7.1118158758588142e-14
diff --git a/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out b/Exec/science/wdmerger/ci-benchmarks/wdmerger_collision_2D.out
@@ -1,29 +1,29 @@
  plotfile = plt00086
  time = 1.25
               variables          minimum value          maximum value
- density                      8.6936468335e-05           19568007.309
- xmom                        -5.4959005475e+14       1.3559838981e+14
- ymom                        -2.5540206053e+15       2.5540214496e+15
+ density                      8.6940338039e-05           19441641.375
+ xmom                        -5.4953770416e+14       1.3594264808e+14
+ ymom                        -2.4933243003e+15       2.4933250525e+15
  zmom                                        0                      0
- rho_E                        7.4987846833e+11       5.0676543192e+24
- rho_e                        7.1083104678e+11       5.0648015049e+24
- Temp                             242292.44287           1409581841.1
- rho_He4                      8.6936468335e-17           3.5973411848
- rho_C12                      3.4774587333e-05            7826946.398
- rho_O16                         5.2161881e-05           11740633.889
- rho_Ne20                     8.6936468335e-17           181819.49413
- rho_Mg24                     8.6936468335e-17           1190.7712386
- rho_Si28                     8.6936468335e-17           6.6817951193
- rho_S32                      8.6936468335e-17       0.00019451843176
- rho_Ar36                     8.6936468335e-17       1.9568007618e-05
- rho_Ca40                     8.6936468335e-17       1.9568007341e-05
- rho_Ti44                     8.6936468335e-17       1.9568007318e-05
- rho_Cr48                     8.6936468335e-17       1.9568007318e-05
- rho_Fe52                     8.6936468335e-17       1.9568007318e-05
- rho_Ni56                     8.6936468335e-17       1.9568007318e-05
- phiGrav                     -5.8709462562e+17      -2.3375498549e+16
- grav_x                          -685026429.13          -51428.265677
- grav_y                          -739654246.49           739654206.24
+ rho_E                        7.4973602186e+11       5.0768248379e+24
+ rho_e                        7.1068648973e+11       5.0744783673e+24
+ Temp                             242282.60874           1404450633.1
+ rho_He4                      8.6940338039e-17            3.398107373
+ rho_C12                      3.4776135215e-05           7775850.9371
+ rho_O16                      5.2164202823e-05           11664450.012
+ rho_Ne20                     8.6940338039e-17             172485.537
+ rho_Mg24                     8.6940338039e-17            1043.054267
+ rho_Si28                     8.6940338039e-17           5.9869391361
+ rho_S32                      8.6940338039e-17       0.00016459247232
+ rho_Ar36                     8.6940338039e-17       1.9441643669e-05
+ rho_Ca40                     8.6940338039e-17       1.9441641397e-05
+ rho_Ti44                     8.6940338039e-17       1.9441641384e-05
+ rho_Cr48                     8.6940338039e-17       1.9441641384e-05
+ rho_Fe52                     8.6940338039e-17       1.9441641384e-05
+ rho_Ni56                     8.6940338039e-17       1.9441641384e-05
+ phiGrav                      -5.870743119e+17       -2.337549858e+16
+ grav_x                           -685044085.4          -51428.268861
+ grav_y                          -739591083.78           739591039.26
  grav_z                                      0                      0
- rho_enuc                    -4.7815621457e+12       7.6360058391e+23
+ rho_enuc                    -7.5385126121e+12       7.1503781318e+23
 
diff --git a/Exec/science/xrb_spherical/analysis/r_profile.py b/Exec/science/xrb_spherical/analysis/r_profile.py
@@ -0,0 +1,56 @@
+#!/usr/bin/env python3
+
+# Spherical R profile at different theta
+
+import os
+import sys
+import yt
+import matplotlib.pyplot as plt
+import numpy as np
+from functools import reduce
+import itertools
+
+import matplotlib.ticker as ptick
+from yt.frontends.boxlib.api import CastroDataset
+from yt.units import cm
+
+
+plotfile = sys.argv[1]
+ds = CastroDataset(plotfile)
+
+rmin = ds.domain_left_edge[0]
+rmax = rmin + 5000.0*cm
+#rmax = ds.domain_right_edge[0]
+print(ds.domain_left_edge[1])
+fig, _ax = plt.subplots(2,2)
+
+axes = list(itertools.chain(*_ax))
+
+fig.set_size_inches(7.0, 8.0)
+
+fields = ["Temp", "density", "x_velocity", "y_velocity"]
+nice_names = [r"$T$ (K)", r"$\rho$ (g/${cm}^3$)", r"$u$ (cm/s)", r"$v$ (cm/s)"]
+
+# 4 rays at different theta values
+thetal = ds.domain_left_edge[1]
+thetar = ds.domain_right_edge[1]
+thetas = [thetal, 0.25*thetar, 0.5*thetar, 0.75*thetar]
+
+for i, f in enumerate(fields):
+
+    for theta in thetas:
+        # simply go from (rmin, theta) -> (rmax, theta). Doesn't need to convert to physical R-Z
+        ray = ds.ray((rmin, theta, 0*cm), (rmax, theta, 0*cm))
+        isrt = np.argsort(ray["t"])
+        axes[i].plot(ray['r'][isrt], ray[f][isrt], label=r"$\theta$ = {:.4f}".format(float(theta)))
+
+    axes[i].set_xlabel(r"$r$ (cm)")
+    axes[i].set_ylabel(nice_names[i])
+    axes[i].set_yscale("symlog")
+
+    if i == 0:
+        axes[0].legend(frameon=False, loc="lower left")
+
+#fig.set_size_inches(10.0, 9.0)
+plt.tight_layout()
+plt.savefig("{}_profiles.png".format(os.path.basename(plotfile)))
diff --git a/Exec/science/xrb_spherical/analysis/slice.py b/Exec/science/xrb_spherical/analysis/slice.py
@@ -7,6 +7,8 @@
 import matplotlib.pyplot as plt
 from yt.frontends.boxlib.api import CastroDataset
 
+from yt.units import cm
+
 """
 Given a plot file and field name, it gives slice plots at the top,
 middle, and bottom of the domain (shell).
@@ -54,8 +56,10 @@ def slice(fname:str, field:str,
                "mid":(r_center*np.sin(theta_center), r_center*np.cos(theta_center)),
                "bot":(r_center*np.sin(thetar)+0.5*width, r_center*np.cos(thetar))}
 
-    sp = yt.SlicePlot(ds, 'phi', field, width=box_widths)
+    # Note we can also set center during SlicePlot, however then we would enter in [r_center, theta_center, 0]
+    # rather than the physical R-Z coordinate if we do it via sp.set_center
 
+    sp = yt.SlicePlot(ds, 'phi', field, width=box_widths)
     sp.set_center(centers[loc])
 
     sp.set_cmap(field, "viridis")

diff --git a/Exec/science/xrb_spherical/inputs.He.1000Hz b/Exec/science/xrb_spherical/inputs.He.1000Hz
@@ -7,25 +7,19 @@ geometry.is_periodic = 0       0
 geometry.coord_sys   = 2                  # 0 => cart, 1 => RZ  2=>spherical
 geometry.prob_lo     = 1.1e6           0
 geometry.prob_hi     = 1.13072e6       3.1415926
-amr.n_cell           = 192             1152
+amr.n_cell           = 768            2304  #192             1152
 
 # >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
 # 0 = Interior           3 = Symmetry
 # 1 = Inflow             4 = SlipWall
 # 2 = Outflow            5 = NoSlipWall
 # >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
-castro.lo_bc       =  1   3   # Inflow in lower R and Symmetry across Theta
+castro.lo_bc       =  3   3   # Inflow in lower R and Symmetry across Theta
 castro.hi_bc       =  2   3   # Outflow in upper boundaries
 
 # Allow non-square zones
 castro.allow_non_unit_aspect_zones = 1
 
-# Hydrostatic condition along R-direction
-castro.xl_ext_bc_type = 1
-castro.hse_interp_temp = 0
-castro.hse_fixed_temp = 1.e8   # this should match problem.T_star
-castro.hse_reflect_vels = 1
-
 # Fill ambient states with outflow velocity in R-direction
 castro.fill_ambient_bc = 1
 castro.ambient_fill_dir = 0
@@ -44,6 +38,7 @@ castro.small_dens = 1.e-5
 castro.ppm_type = 1
 castro.grav_source_type = 2
 castro.use_pslope = 1
+castro.ppm_well_balanced = 1
 castro.pslope_cutoff_density = 1.e4
 
 gravity.gravity_type = ConstantGrav
@@ -144,7 +139,6 @@ problem.X_min = 1.e-2
 
 problem.r_refine_distance = 9.216e4
 
-
 # Microphysics
 
 integrator.rtol_spec = 1.e-6

diff --git a/Source/driver/Castro.cpp b/Source/driver/Castro.cpp
@@ -306,14 +306,26 @@ Castro::read_params ()
 #elif (AMREX_SPACEDIM == 2)
     if ( dgeom.IsSPHERICAL() )
       {
+        if ( (lo_bc[1] != amrex::PhysBCType::symmetry) && (dgeom.ProbLo(1) == 0.0) )
+        {
+          std::cerr << "ERROR:Castro::read_params: must set theta=0 boundary condition to Symmetry for spherical\n";
+          amrex::Error();
+        }
+
+        if ( (hi_bc[1] != amrex::PhysBCType::symmetry) && (std::abs(dgeom.ProbHi(1) - M_PI) <= 1.e-4_rt) )
+        {
+          std::cerr << "ERROR:Castro::read_params: must set theta=pi boundary condition to Symmetry for spherical\n";
+          amrex::Error();
+        }
+
         if ( (dgeom.ProbLo(1) < 0.0_rt) && (dgeom.ProbHi(1) > M_PI) )
         {
-          amrex::Abort("Theta must be within [0, Pi] for spherical coordinate system in 2D");
+          amrex::Abort("ERROR:Castro::read_params: Theta must be within [0, Pi] for spherical coordinate system in 2D");
         }
 
         if ( dgeom.ProbLo(0) < static_cast<Real>(NUM_GROW) * dgeom.CellSize(0) )
         {
-          amrex::Abort("R-min must be large enough so ghost cells doesn't extend to negative R");
+          amrex::Abort("ERROR:Castro::read_params: R-min must be large enough so ghost cells doesn't extend to negative R");
         }
       }
 #elif (AMREX_SPACEDIM == 3)

diff --git a/Source/driver/_cpp_parameters b/Source/driver/_cpp_parameters
@@ -161,6 +161,10 @@ dual_energy_eta2             Real          1.0e-4
 # does well with HSE
 use_pslope                   bool           0
 
+# for PPM, do we only use the perturbational pressure in the characteristic
+# tracing?  This is more indepth than the simple `use_pslope` approach.
+ppm_well_balanced            bool           0
+
 # if we are using pslope, below what density to we turn off the well-balanced
 # reconstruction?
 pslope_cutoff_density        Real          -1.e20