Merge branch 'main' into add_small_e_prot

CHANGES.md

@@ -1,3 +1,52 @@
+# 4.5.0
+
+  * add protection against negative density/energy in the 4th order
+    compressible method when converting from averages to centers
+    (#309)
+
+  * implement a "clean state" routine for the compressible solvers
+    that can enforce a density floor (#329)
+
+  * updated the citation information (#324) and the zenodo file (#325)
+
+  * code cleaning (#323)
+
+  * added python 3.13 support (#321)
+
+  * CI doc build fixes (#322) and added the ability to skip executing
+    notebooks when building docs (#319)
+
+  * fixed issues with how derived variables are recognized when
+    reading in an output file and add vorticity and machnumber
+    to compressible derives (#294, #320)
+
+  * expand the documentation on the compressible solvers (#317)
+
+  * remove unused APIs from docs (#318) and fix autodoc issues (#316)
+
+  * implement a sponge term in the compressible solvers (#313, #336)
+
+  * add some asserts to the compressible solver to catch unphysical
+    mass and energy (#312)
+
+  * add the compressible convection problem (#293, #310, #315, #337,
+    #338, #339)
+
+  * add a multimode Rayleigh-Taylor problem to the compressible
+    solvers (#307)
+
+  * improve the documentation on working with data (#304, #305)
+
+  * spherical compressible solver fixes (#295, #296, #298, #299, #300)
+
+  * move the 4th order reconstruction into mesh (#302)
+
+  * use a safer divide in compressible when density is zero (#297)
+
+  * add the ability to have a problem-dependent external source and
+    clean-up how the compressible sources are treated (#289, #290,
+    #292)
+
 # 4.4.0
 
   * the gresho problem (compressible) now takes Mach number as input

pyro/compressible/interface.py

@@ -3,7 +3,7 @@
 
 
 @njit(cache=True)
-def states(idir, ng, dx, dt,
+def states(idir, ng, dx, dloga, dt,
            irho, iu, iv, ip, ix, nspec,
            gamma, qv, dqv):
     r"""
@@ -212,6 +212,27 @@ def states(idir, ng, dx, dt,
                     q_l[i, j + 1, m] = q_l[i, j + 1, m] + sum_l
                     q_r[i, j,  m] = q_r[i, j,  m] + sum_r
 
+            # Geometric Source term from converting conserved-variable to primitive
+            # It's only there for non Cartesian coord.
+
+            if idir == 1:
+                rho_source = -0.5 * dt * dloga[i, j] * q[irho] * q[iu]
+
+                q_l[i + 1, j, irho] += rho_source
+                q_r[i, j, irho] += rho_source
+
+                q_l[i + 1, j, ip] += rho_source * cs * cs
+                q_r[i, j, ip] += rho_source * cs * cs
+
+            else:
+                rho_source = -0.5 * dt * dloga[i, j] * q[irho] * q[iv]
+
+                q_l[i, j + 1, irho] += rho_source
+                q_r[i, j, irho] += rho_source
+
+                q_l[i, j + 1, ip] += rho_source * cs * cs
+                q_r[i, j, ip] += rho_source * cs * cs
+
     return q_l, q_r
 
 

pyro/compressible/problems/inputs.convection

@@ -1,7 +1,7 @@
 # simple inputs files for the four-corner problem.
 
 [driver]
-max_steps = 10000
+max_steps = 100000
 tmax = 25.0
 
 

pyro/compressible/problems/inputs.convection.big

@@ -1,7 +1,7 @@
 # simple inputs files for the four-corner problem.
 
 [driver]
-max_steps = 10000
+max_steps = 100000
 tmax = 25.0
 
 

pyro/compressible/unsplit_fluxes.py

@@ -206,11 +206,13 @@ def interface_states(my_data, rp, ivars, tc, dt):
     tm_states = tc.timer("interfaceStates")
     tm_states.begin()
 
-    V_l, V_r = ifc.states(1, myg.ng, myg.Lx, dt,
+    _V_l, _V_r = ifc.states(1, myg.ng, myg.Lx, myg.dlogAx, dt,
                           ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix,
                           ivars.naux,
                           gamma,
                           q, ldx)
+    V_l = ai.ArrayIndexer(d=_V_l, grid=myg)
+    V_r = ai.ArrayIndexer(d=_V_r, grid=myg)
 
     tm_states.end()
 
@@ -225,7 +227,7 @@ def interface_states(my_data, rp, ivars, tc, dt):
     # left and right primitive variable states
     tm_states.begin()
 
-    _V_l, _V_r = ifc.states(2, myg.ng, myg.Ly, dt,
+    _V_l, _V_r = ifc.states(2, myg.ng, myg.Ly, myg.dlogAy, dt,
                             ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix,
                             ivars.naux,
                             gamma,

pyro/compressible_fv4/simulation.py

@@ -1,3 +1,5 @@
+import numpy as np
+
 import pyro.compressible_fv4.fluxes as flx
 from pyro import compressible_rk
 from pyro.compressible import cons_to_prim, get_external_sources, get_sponge_factor
@@ -70,6 +72,9 @@ def preevolve(self):
         this solver assumes that the initialization was done to
         cell-centers and converts it to cell-averages."""
 
+        # this should only work for dx == dy
+        assert np.abs(self.cc_data.grid.dx - self.cc_data.grid.dy) < 1.e-12 * self.cc_data.grid.dx, "grid cells need to be square"
+
         # we just initialized cell-centers, but we need to store averages
         for var in self.cc_data.names:
             self.cc_data.from_centers(var)