diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 6ba371706..45aafae5d 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -68,6 +68,13 @@ def cons_to_prim(U, gamma, ivars, myg): e_min = e.v().min() rho_min = q.v(n=ivars.irho).min() + if e_min < 0: + eidx = np.asarray(e < 0).nonzero() + i_idx = eidx[0] + j_jdx = eidx[1] + for i, j in zip(i_idx, j_jdx): + print(f" e < 0: {i}, {j}, {e[i, j]}") + assert e_min > 0.0 and rho_min > 0.0, f"invalid state, min(rho) = {rho_min}, min(e) = {e_min}" q[:, :, ivars.ip] = eos.pres(gamma, q[:, :, ivars.irho], e) @@ -452,8 +459,19 @@ def evolve(self): def clean_state(self, U): """enforce minimum density and eint on the conserved state U""" - U.v(n=self.ivars.idens)[:, :] = np.maximum(U.v(n=self.ivars.idens), - self.rp.get_param("compressible.small_dens")) + U[..., self.ivars.idens] = np.maximum(U[..., self.ivars.idens], + self.rp.get_param("compressible.small_dens")) + + if self.small_eint > 0: + + ekin = 0.5 * (U[..., self.ivars.ixmom]**2 + + U[..., self.ivars.iymom]**2) / U[..., self.ivars.idens] + + rhoe = U[..., self.ivars.iener] - ekin + + U[..., self.ivars.iener] = np.where(rhoe < self.small_eint, + U[..., self.ivars.idens] * self.small_eint + ekin, + U[..., self.ivars.iener]) def dovis(self): """ diff --git a/pyro/compressible_fv4/_defaults b/pyro/compressible_fv4/_defaults index ae1c2979c..7cd1eb710 100644 --- a/pyro/compressible_fv4/_defaults +++ b/pyro/compressible_fv4/_defaults @@ -24,6 +24,7 @@ grav = 0.0 ; gravitational acceleration (in y-direction) riemann = CGF small_dens = -1.e200 ; minimum allowed density +small_eint_factor = -1.e200 ; multiplicative factor on the initial minimum eint to use for limiting e [sponge] do_sponge = 0 ; do we include a sponge source term diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index bca1522f3..8fd5685b1 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -1,6 +1,6 @@ import pyro.compressible_fv4.fluxes as flx from pyro import compressible_rk -from pyro.compressible import get_external_sources, get_sponge_factor +from pyro.compressible import cons_to_prim, get_external_sources, get_sponge_factor from pyro.mesh import fv @@ -73,3 +73,12 @@ def preevolve(self): # we just initialized cell-centers, but we need to store averages for var in self.cc_data.names: self.cc_data.from_centers(var) + + efac = self.rp.get_param("compressible.small_eint_factor") + + print(type(self.cc_data)) + + q = cons_to_prim(self.cc_data.data, self.rp.get_param("eos.gamma"), self.ivars, self.cc_data.grid) + + self.small_eint = efac * q.v(n=self.ivars.ip).min() + print("small_eint = ", self.small_eint) diff --git a/pyro/compressible_sdc/_defaults b/pyro/compressible_sdc/_defaults index 11568dd35..837e73068 100644 --- a/pyro/compressible_sdc/_defaults +++ b/pyro/compressible_sdc/_defaults @@ -21,6 +21,8 @@ riemann = CGF small_dens = -1.e200 ; minimum allowed density +small_eint_factor = -1.e200 ; multiplicative factor on the initial minimum eint to use for limiting e + [sponge] do_sponge = 0 ; do we include a sponge source term