From f248336b2b7679b2e6f57697abf3096b63b98306 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 8 Jan 2025 15:27:27 -0500 Subject: [PATCH 1/5] add some protection against negatives in 4th order --- pyro/compressible_fv4/fluxes.py | 13 +++++++++++-- pyro/mesh/fv.py | 8 +++++++- pyro/pyro_sim.py | 2 ++ 3 files changed, 20 insertions(+), 3 deletions(-) diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index 6b3d67528..ff29c1775 100644 --- a/pyro/compressible_fv4/fluxes.py +++ b/pyro/compressible_fv4/fluxes.py @@ -52,10 +52,10 @@ def fluxes(myd, rp, ivars): # convert U from cell-averages to cell-centers U_cc = np.zeros_like(U_avg) - U_cc[:, :, ivars.idens] = myd.to_centers("density") + U_cc[:, :, ivars.idens] = myd.to_centers("density", is_positive=True) U_cc[:, :, ivars.ixmom] = myd.to_centers("x-momentum") U_cc[:, :, ivars.iymom] = myd.to_centers("y-momentum") - U_cc[:, :, ivars.iener] = myd.to_centers("energy") + U_cc[:, :, ivars.iener] = myd.to_centers("energy", is_positive=True) # compute the primitive variables of both the cell-center and averages q_bar = comp.cons_to_prim(U_avg, gamma, ivars, myd.grid) @@ -66,6 +66,15 @@ def fluxes(myd, rp, ivars): for n in range(ivars.nq): q_avg.v(n=n, buf=3)[:, :] = q_cc.v(n=n, buf=3) + myg.dx**2/24.0 * q_bar.lap(n=n, buf=3) + # enforce positivity + q_avg.v(n=ivars.irho, buf=3)[:, :] = np.where(q_avg.v(n=ivars.irho, buf=3) > 0, + q_avg.v(n=ivars.irho, buf=3), + q_cc.v(n=ivars.irho, buf=3)) + + q_avg.v(n=ivars.ip, buf=3)[:, :] = np.where(q_avg.v(n=ivars.ip, buf=3) > 0, + q_avg.v(n=ivars.ip, buf=3), + q_cc.v(n=ivars.ip, buf=3)) + # flattening -- there is a single flattening coefficient (xi) for all directions use_flattening = rp.get_param("compressible.use_flattening") diff --git a/pyro/mesh/fv.py b/pyro/mesh/fv.py index b5d0fafbd..816c86709 100644 --- a/pyro/mesh/fv.py +++ b/pyro/mesh/fv.py @@ -3,6 +3,8 @@ """ +import numpy as np + from pyro.mesh.patch import CellCenterData2d @@ -13,7 +15,7 @@ class FV2d(CellCenterData2d): """ - def to_centers(self, name): + def to_centers(self, name, is_positive=False): """ convert variable name from an average to cell-centers """ a = self.get_var(name) @@ -21,6 +23,10 @@ def to_centers(self, name): ng = self.grid.ng c[:, :] = a[:, :] c.v(buf=ng-1)[:, :] = a.v(buf=ng-1) - self.grid.dx**2*a.lap(buf=ng-1)/24.0 + + if is_positive: + c.v(buf=ng-1)[:, :] = np.where(c.v(buf=ng-1) >= 0.0, c.v(buf=ng-1), a.v(buf=ng-1)) + return c def from_centers(self, name): diff --git a/pyro/pyro_sim.py b/pyro/pyro_sim.py index bdaf6efa9..60e21b989 100755 --- a/pyro/pyro_sim.py +++ b/pyro/pyro_sim.py @@ -30,6 +30,8 @@ "lm_atm", "swe"] +import warnings +warnings.filterwarnings("error") class Pyro: """ From 504f5dd833ad3e9289747c3cb24d52c734372e18 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 9 Jan 2025 11:28:36 -0500 Subject: [PATCH 2/5] revert --- pyro/pyro_sim.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyro/pyro_sim.py b/pyro/pyro_sim.py index 60e21b989..bdaf6efa9 100755 --- a/pyro/pyro_sim.py +++ b/pyro/pyro_sim.py @@ -30,8 +30,6 @@ "lm_atm", "swe"] -import warnings -warnings.filterwarnings("error") class Pyro: """ From b60ed1631c695fc07a4c53609153bee82053b579 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 25 Jan 2025 09:45:26 -0500 Subject: [PATCH 3/5] this now uses a mask --- pyro/compressible_fv4/fluxes.py | 20 ++++++++++++++++++-- pyro/mesh/patch.py | 6 +++--- 2 files changed, 21 insertions(+), 5 deletions(-) diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index ff29c1775..00bbfd50e 100644 --- a/pyro/compressible_fv4/fluxes.py +++ b/pyro/compressible_fv4/fluxes.py @@ -52,10 +52,26 @@ def fluxes(myd, rp, ivars): # convert U from cell-averages to cell-centers U_cc = np.zeros_like(U_avg) - U_cc[:, :, ivars.idens] = myd.to_centers("density", is_positive=True) + + U_cc[:, :, ivars.idens] = myd.to_centers("density") U_cc[:, :, ivars.ixmom] = myd.to_centers("x-momentum") U_cc[:, :, ivars.iymom] = myd.to_centers("y-momentum") - U_cc[:, :, ivars.iener] = myd.to_centers("energy", is_positive=True) + U_cc[:, :, ivars.iener] = myd.to_centers("energy") + + # the mask will be 1 in any zone where the density or energy + # is unphysical do to the conversion from averages to centers + + rhoe = U_cc[..., ivars.iener] - 0.5 * (U_cc[..., ivars.ixmom]**2 + + U_cc[..., ivars.iymom]**2) / U_cc[..., ivars.idens] + + mask = myg.scratch_array(dtype=np.uint8) + mask[:, :] = np.where(np.logical_or(U_cc[:, :, ivars.idens] < 0, rhoe < 0), + 1, 0) + + U_cc[..., ivars.idens] = np.where(mask == 1, U_avg[..., ivars.idens], U_cc[..., ivars.idens]) + U_cc[..., ivars.ixmom] = np.where(mask == 1, U_avg[..., ivars.ixmom], U_cc[..., ivars.ixmom]) + U_cc[..., ivars.iymom] = np.where(mask == 1, U_avg[..., ivars.iymom], U_cc[..., ivars.iymom]) + U_cc[..., ivars.iener] = np.where(mask == 1, U_avg[..., ivars.iener], U_cc[..., ivars.iener]) # compute the primitive variables of both the cell-center and averages q_bar = comp.cons_to_prim(U_avg, gamma, ivars, myd.grid) diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py index bdf4d6687..ab5159f3e 100644 --- a/pyro/mesh/patch.py +++ b/pyro/mesh/patch.py @@ -146,15 +146,15 @@ def __init__(self, nx, ny, *, ng=1, self.xr2d = ArrayIndexer(d=xr2d, grid=self) self.yr2d = ArrayIndexer(d=yr2d, grid=self) - def scratch_array(self, *, nvar=1): + def scratch_array(self, *, nvar=1, dtype=np.float64): """ return a standard numpy array dimensioned to have the size and number of ghostcells as the parent grid """ if nvar == 1: - _tmp = np.zeros((self.qx, self.qy), dtype=np.float64) + _tmp = np.zeros((self.qx, self.qy), dtype=dtype) else: - _tmp = np.zeros((self.qx, self.qy, nvar), dtype=np.float64) + _tmp = np.zeros((self.qx, self.qy, nvar), dtype=dtype) return ArrayIndexer(d=_tmp, grid=self) def coarse_like(self, N): From 5ee8ca9539bffb484c815391ee5cdde89a678cef Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 25 Jan 2025 09:50:00 -0500 Subject: [PATCH 4/5] fix flake9 --- pyro/compressible_fv4/fluxes.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index 00bbfd50e..d59352544 100644 --- a/pyro/compressible_fv4/fluxes.py +++ b/pyro/compressible_fv4/fluxes.py @@ -52,7 +52,6 @@ def fluxes(myd, rp, ivars): # convert U from cell-averages to cell-centers U_cc = np.zeros_like(U_avg) - U_cc[:, :, ivars.idens] = myd.to_centers("density") U_cc[:, :, ivars.ixmom] = myd.to_centers("x-momentum") U_cc[:, :, ivars.iymom] = myd.to_centers("y-momentum") From 0c9b5a5bf026d1ae4ae65c30c4ed230384949ff9 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Sat, 25 Jan 2025 13:47:40 -0500 Subject: [PATCH 5/5] add small e protection --- pyro/compressible/simulation.py | 22 ++++++++++++++++++++-- pyro/compressible_fv4/_defaults | 1 + pyro/compressible_fv4/simulation.py | 11 ++++++++++- pyro/compressible_sdc/_defaults | 2 ++ 4 files changed, 33 insertions(+), 3 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 5f03f68e4..bd5f39b04 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) @@ -442,8 +449,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 2fff3541d..9d5d3bb4a 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 @@ -68,3 +68,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