Merge branch 'development' into 2d_spherical_rotation

CHANGES.md

@@ -1,3 +1,24 @@
+# 24.11
+
+  * a new well-balanced method was added to the CTU PPM solver.  This
+    does the characteristic projection only on the perturbed pressure
+    and then adds back in the hydrostatic pressure. It can be enabled
+    via `castro.ppm_well_balanced` (#2945)
+
+  * fixed a bug in the div{U} calculation for artificial viscosity on
+    symmetry boundaries (#2983)
+
+  * more development on 2D spherical geometry (#2973, #2975, #2981)
+
+  * updates to the massive star plotting scripts (#2979)
+
+  * add some new checks to prevent running unsupported combinations
+    of solvers (#2978)
+
+  * documentation updates (#2977)
+
+  * `flame_wave` can now be run in 1D (#2976)
+
 # 24.10
 
   * update initial model for `subchandra` when doing ASE NSE (#2970)

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).
-
-

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.

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}
+

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

Exec/radiation_tests/Rad2Tshock/ci-benchmarks/Rad2TShock-1d.out

@@ -1,40 +1,40 @@
  plotfile = plt_00025
- time = 0.00026083835470476599
+ time = 0.00026084061812361401
               variables          minimum value          maximum value
- density                      5.4596904436e-13       1.2720098871e-12
- xmom                         1.2327065163e-07       1.4587103381e-07
+ density                      5.4596904436e-13       1.2717694091e-12
+ xmom                         1.2328047524e-07       1.4591226172e-07
  ymom                                        0                      0
  zmom                                        0                      0
- rho_E                          0.021940622286         0.039199658391
- rho_e                         0.0068091599527         0.032333529333
- Temp                             100.00003116            217.3564598
- rho_X                        5.4596904436e-13       1.2720098871e-12
- rad                          7.5662679486e-07       1.4083425328e-05
- pressure                      0.0045394399687         0.021555686223
- kineng                        0.0064467244856         0.015131462333
- soundspeed                        117717.6283           173551.23637
+ rho_E                          0.021940622285         0.039194728059
+ rho_e                         0.0068091599514         0.032333527949
+ Temp                             100.00003115           217.29660609
+ rho_X                        5.4596904436e-13       1.2717694091e-12
+ rad                          7.5662675907e-07       1.4083418723e-05
+ pressure                      0.0045394399678           0.0215556853
+ kineng                        0.0064473434191         0.015131462334
+ soundspeed                       117717.62829           173527.33922
  Gamma_1                          1.6666666667           1.6666666667
- MachNumber                      0.60689740202           1.9999997205
- uplusc                           269483.97988           362284.15407
- uminusc                         -66810.704217            117717.5954
- entropy                          2458725989.8             2507243428
+ MachNumber                       0.6068973935           1.9999997207
+ uplusc                           269518.17809           362247.33487
+ uminusc                         -66804.446805           117717.59541
+ entropy                          2458725989.8           2507184972.2
  magvort                                     0                      0
- divu                            -20286.061564           245.50876324
- eint_E                            12471696011            27108028479
- eint_e                            12471696011            27108028479
- logden                           -12.26283198          -11.895509513
- StateErr_0                   5.4596904436e-13       1.2720098871e-12
- StateErr_1                       100.00003116            217.3564598
+ divu                            -20283.836079           245.94831837
+ eint_E                            12471696009            27100563708
+ eint_e                            12471696009            27100563708
+ logden                           -12.26283198          -11.895591626
+ StateErr_0                   5.4596904436e-13       1.2717694091e-12
+ StateErr_1                       100.00003115           217.29660609
  StateErr_2                                  1                      1
  X(X)                                        1                      1
  abar                                        1                      1
- x_velocity                       102258.95554            235435.2237
+ x_velocity                       102299.57948           235435.22371
  y_velocity                                  0                      0
  z_velocity                                  0                      0
- magvel                           102258.95554            235435.2237
- radvel                           -235435.2237           114038.34336
- circvel                                     0         0.005524271728
- magmom                       1.2327065163e-07       1.4587103381e-07
+ magvel                           102299.57948           235435.22371
+ radvel                          -235435.22371           114042.93363
+ circvel                                     0        0.0047841596539
+ magmom                       1.2328047524e-07       1.4591226172e-07
  angular_momentum_x                          0                      0
  angular_momentum_y                          0                      0
  angular_momentum_z                         -0                     -0

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

Exec/science/xrb_spherical/GNUmakefile

@@ -0,0 +1,40 @@
+PRECISION  = DOUBLE
+PROFILE    = FALSE
+
+DEBUG      = FALSE
+
+DIM        = 2
+
+COMP       = gnu
+
+USE_MPI    = TRUE
+
+USE_GRAV   = TRUE
+USE_REACT  = FALSE
+
+USE_ROTATION  = FALSE
+USE_DIFFUSION = FALSE
+
+# define the location of the CASTRO top directory
+CASTRO_HOME  ?= ../../..
+
+USE_JACOBIAN_CACHING = TRUE
+USE_MODEL_PARSER = TRUE
+NUM_MODELS := 2
+
+# This sets the EOS directory in $(MICROPHYSICS_HOME)/eos
+EOS_DIR     := helmholtz
+
+# This sets the network directory in $(MICROPHYSICS_HOME)/networks
+NETWORK_DIR := subch_base
+
+INTEGRATOR_DIR := VODE
+
+CONDUCTIVITY_DIR := stellar
+
+PROBLEM_DIR ?= ./
+
+Bpack   := $(PROBLEM_DIR)/Make.package
+Blocs   := $(PROBLEM_DIR)
+
+include $(CASTRO_HOME)/Exec/Make.Castro

Exec/science/xrb_spherical/Make.package

@@ -0,0 +1,2 @@
+CEXE_headers += initial_model.H
+

Exec/science/xrb_spherical/README.md

@@ -0,0 +1,6 @@
+# xrb_spherical
+
+This is the full-star XRB flame setup based on flame_wave.
+This setup uses a spherical 2D geometry to model XRB flame
+on a spherical shell with initial temperature perturbation
+on the north pole.

Exec/science/xrb_spherical/_prob_params

@@ -0,0 +1,68 @@
+
+dtemp                   real              3.81e8_rt                 y
+
+theta_half_max          real              1.745e-2_rt               y
+
+theta_half_width        real              4.9e-3_rt                 y
+
+# cutoff mass fraction of the first species for refinement
+X_min                   real              1.e-4_rt                  y
+
+# do we dynamically refine based on density? or based on height?
+tag_by_density          integer           1                         y
+
+# used for tagging if tag_by_density = 1
+cutoff_density          real              500.e0_rt                 y
+
+# used if we are refining based on height rather than density
+refine_height           real              3600                      y
+
+T_hi                    real              5.e8_rt                   y
+
+T_star                  real              1.e8_rt                   y
+
+T_lo                    real              5.e7_rt                   y
+
+dens_base               real              2.e6_rt                   y
+
+H_star                  real              500.e0_rt                 y
+
+atm_delta               real              25.e0_rt                  y
+
+fuel1_name              string         "helium-4"                y
+
+fuel2_name              string         ""                        y
+
+fuel3_name              string         ""                        y
+
+fuel4_name              string         ""                        y
+
+ash1_name               string         "iron-56"                 y
+
+ash2_name               string         ""                        y
+
+ash3_name               string         ""                        y
+
+fuel1_frac              real              1.0_rt                    y
+
+fuel2_frac              real              0.0_rt                    y
+
+fuel3_frac              real              0.0_rt                    y
+
+fuel4_frac              real              0.0_rt                    y
+
+ash1_frac               real              1.0_rt                    y
+
+ash2_frac               real              0.0_rt                    y
+
+ash3_frac               real              0.0_rt                    y
+
+low_density_cutoff      real              1.e-4_rt                  y
+
+smallx                  real              1.e-10_rt                 y
+
+r_refine_distance       real              1.e30_rt                  y
+
+max_hse_tagging_level   integer           2                         y
+
+max_base_tagging_level  integer           1                         y

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)))