From f248336b2b7679b2e6f57697abf3096b63b98306 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Wed, 8 Jan 2025 15:27:27 -0500 Subject: [PATCH 1/4] 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/4] 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/4] 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/4] 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")