From 1ca7ccb710211d4ea5c3cb245065b135a448a238 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 08:23:19 -0400 Subject: [PATCH 01/27] compressible: add the ability for a problem-dependent external source --- pyro/compressible/problems/heating.py | 68 +++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 pyro/compressible/problems/heating.py diff --git a/pyro/compressible/problems/heating.py b/pyro/compressible/problems/heating.py new file mode 100644 index 000000000..bb2d568f1 --- /dev/null +++ b/pyro/compressible/problems/heating.py @@ -0,0 +1,68 @@ +"""A test of the energy sources. This uses a uniform domain and +slowly adds heat to the center over time.""" + +import math + +import numpy as np + +from pyro.util import msg + +DEFAULT_INPUTS = "inputs.sedov" + +PROBLEM_PARAMS = {"heating.rho_ambient": 1.0, # ambient density + "heating.p_ambient": 10.0, # ambient pressure + "heating.r_src": 0.1, # physical size of the heating src + "heating.e_rate": 0.1} # energy generation rate (energy / mass / time) + + +def init_data(my_data, rp): + """ initialize the heating problem """ + + if rp.get_param("driver.verbose"): + msg.bold("initializing the heating problem...") + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + gamma = rp.get_param("eos.gamma") + + # initialize the components, remember, that ener here is rho*eint + # + 0.5*rho*v**2, where eint is the specific internal energy + # (erg/g) + dens[:, :] = rp.get_param("heating.rho_ambient") + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + ener[:, :] = rp.get_param("heating.p_ambient") / (gamma - 1.0) + + +def source_terms(my_data, rp): + """source terms to be added to the evolution""" + + e_src = my_data.get_var("E_src") + + grid = my_data.grid + + xctr = 0.5 * (grid.xmin + grid.xmax) + yctr = 0.5 * (grid.ymin + grid.ymax) + + dist = np.sqrt((my_data.grid.x2d - xctr)**2 + + (my_data.grid.y2d - yctr)**2) + + e_rate = rp.get_param("heating.e_rate") + r_src = rp.get_param("heating.r_src") + + e_src[:, :] = e_rate * np.exp(-(dist / r_src)**2) + + +def finalize(): + """ print out any information to the user at the end of the run """ + + print(""" + The script analysis/sedov_compare.py can be used to analyze these + results. That will perform an average at constant radius and + compare the radial profiles to the exact solution. Sample exact + data is provided as analysis/cylindrical-sedov.out + """) From 6e05e62f6f75ed0106cd0e22d3ff2ddaab70da69 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 08:45:51 -0400 Subject: [PATCH 02/27] compressible: add a get_external_sources method this will make it easier to do external sources problem-by-problem --- pyro/compressible/__init__.py | 2 +- pyro/compressible/simulation.py | 20 ++++++++++++++++++++ pyro/compressible/unsplit_fluxes.py | 28 +++++++++++----------------- 3 files changed, 32 insertions(+), 18 deletions(-) diff --git a/pyro/compressible/__init__.py b/pyro/compressible/__init__.py index d0381535f..91b3411c0 100644 --- a/pyro/compressible/__init__.py +++ b/pyro/compressible/__init__.py @@ -6,4 +6,4 @@ __all__ = ["simulation"] -from .simulation import Simulation, Variables, cons_to_prim, prim_to_cons +from .simulation import Simulation, Variables, cons_to_prim, prim_to_cons, get_external_sources diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 843285611..8b7cbb1e4 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -91,6 +91,26 @@ def prim_to_cons(q, gamma, ivars, myg): return U +def get_external_sources(t, U, ivars, rp, myg): + """compute the external sources, including gravity""" + + U_src = myg.scratch_array(nvar=ivars.nvar) + + grav = rp.get_param("compressible.grav") + + if myg.coord_type == 1: + # gravity points in the radial direction for spherical + U_src[:, :, ivars.ixmom] = U[:, :, ivars.idens] * grav + U_src[:, :, ivars.iener] = U[:, :, ivars.ixmom] * grav + + else: + # gravity points in the vertical (y) direction for Cartesian + U_src[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav + U_src[:, :, ivars.iener] = U[:, :, ivars.iymom] * grav + + return U_src + + class Simulation(NullSimulation): """The main simulation class for the corner transport upwind compressible hydrodynamics solver diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index 5c1fd5462..483b2d50a 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -293,25 +293,19 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr, ymom_src = my_aux.get_var("ymom_src") E_src = my_aux.get_var("E_src") - grav = rp.get_param("compressible.grav") + U_src = comp.get_external_sources(my_data.t, my_data.data, ivars, rp, myg) - # Calculate external source (gravity), geometric, and pressure terms - if myg.coord_type == 1: - # assume gravity points in r-direction in spherical. - dens_src.v()[:, :] = 0.0 - xmom_src.v()[:, :] = dens.v()*grav + \ - ymom.v()**2 / (dens.v()*myg.x2d.v()) - \ - (pres.ip(1) - pres.v()) / myg.Lx.v() - ymom_src.v()[:, :] = -(pres.jp(1) - pres.v()) / myg.Ly.v() - \ - xmom.v()*ymom.v() / (dens.v()*myg.x2d.v()) - E_src.v()[:, :] = xmom.v()*grav + dens_src[:, :] = U_src[:, :, ivars.idens] + xmom_src[:, :] = U_src[:, :, ivars.ixmom] + ymom_src[:, :] = U_src[:, :, ivars.iymom] + E_src[:, :] = U_src[:, :, ivars.iener] - else: - # assume gravity points in y-direction in cartesian - dens_src.v()[:, :] = 0.0 - xmom_src.v()[:, :] = 0.0 - ymom_src.v()[:, :] = dens.v()*grav - E_src.v()[:, :] = ymom.v()*grav + # add geometric and pressure terms + if myg.coord_type == 1: + xmom_src.v()[:, :] += ymom.v()**2 / (dens.v()*myg.x2d.v()) - \ + (pres.ip(1) - pres.v()) / myg.Lx.v() + ymom_src.v()[:, :] += -(pres.jp(1) - pres.v()) / myg.Ly.v() - \ + xmom.v()*ymom.v() / (dens.v()*myg.x2d.v()) my_aux.fill_BC("dens_src") my_aux.fill_BC("xmom_src") From 47b5eb269941dbd69602ebe7181530ec6e08dbf4 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 08:58:07 -0400 Subject: [PATCH 03/27] fix linting --- pyro/compressible/__init__.py | 3 ++- pyro/compressible/simulation.py | 2 ++ 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/pyro/compressible/__init__.py b/pyro/compressible/__init__.py index 91b3411c0..d78b69b9b 100644 --- a/pyro/compressible/__init__.py +++ b/pyro/compressible/__init__.py @@ -6,4 +6,5 @@ __all__ = ["simulation"] -from .simulation import Simulation, Variables, cons_to_prim, prim_to_cons, get_external_sources +from .simulation import (Simulation, Variables, cons_to_prim, + get_external_sources, prim_to_cons) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 8b7cbb1e4..f7ea73e58 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -94,6 +94,8 @@ def prim_to_cons(q, gamma, ivars, myg): def get_external_sources(t, U, ivars, rp, myg): """compute the external sources, including gravity""" + _ = t # maybe unused + U_src = myg.scratch_array(nvar=ivars.nvar) grav = rp.get_param("compressible.grav") From bb483b00312027804e622b750c43d1a7272a9de2 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 09:10:45 -0400 Subject: [PATCH 04/27] create a U_old --- pyro/compressible/simulation.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index f7ea73e58..78e3ec8ee 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -290,9 +290,11 @@ def evolve(self): self.cc_data, self.rp, self.ivars) - old_dens = dens.copy() - old_xmom = xmom.copy() - old_ymom = ymom.copy() + # save the old state (without ghost cells) + U_old = np.stack([dens.v().copy(), + xmom.v().copy(), + ymom.v().copy(), + ener.v().copy()], axis=2) # Conservative update @@ -308,32 +310,30 @@ def evolve(self): # Now apply external sources - # For SphericalPolar (coord_type == 1): - # There are gravity (external) sources, - # geometric terms due to local unit vectors, and pressure gradient - # since we don't include pressure in xmom and ymom fluxes - # due to incompatible divergence and gradient in non-Cartesian geometry + # For SphericalPolar (coord_type == 1) there are geometric + # terms and the pressure gradient since we don't include + # pressure in xmom and ymom fluxes # For Cartesian2d (coord_type == 0): # There is only gravity sources. if myg.coord_type == 1: xmom.v()[:, :] += 0.5*self.dt * \ - ((dens.v() + old_dens.v())*grav + + ((dens.v() + U_old[:, :, self.ivars.ixmom])*grav + (ymom.v()**2 / dens.v() + - old_ymom.v()**2 / old_dens.v()) / myg.x2d.v()) - \ + U_old[:, :, self.ivars.iymom]**2 / U_old[:. :, self.ivars.idens]) / myg.x2d.v()) - \ self.dt * (qx.ip(1, n=self.ivars.ip) - qx.v(n=self.ivars.ip)) / myg.Lx.v() ymom.v()[:, :] += 0.5*self.dt * \ (-xmom.v()*ymom.v() / dens.v() - - old_xmom.v()*old_ymom.v() / old_dens.v()) / myg.x2d.v() - \ + U_old[:, :, self.ivars.ixmom] * U_old[:, :, self.ivars.iymom] / U_old[:, :, self.ivars.idens]) / myg.x2d.v() - \ self.dt * (qy.jp(1, n=self.ivars.ip) - qy.v(n=self.ivars.ip)) / myg.Ly.v() - ener.v()[:, :] += 0.5*self.dt*(xmom.v() + old_xmom.v())*grav + ener.v()[:, :] += 0.5 * self.dt * (xmom.v() + U_old[:, :, self.ivars.ixmom]) * grav else: - ymom.v()[:, :] += 0.5*self.dt*(dens.v() + old_dens.v())*grav - ener.v()[:, :] += 0.5*self.dt*(ymom.v() + old_ymom.v())*grav + ymom.v()[:, :] += 0.5 * self.dt * (dens.v() + U_old[:, :, self.ivars.idens]) * grav + ener.v()[:, :] += 0.5 * self.dt * (ymom.v() + U_old[:, :, self.ivars.iymom]) * grav if self.particles is not None: self.particles.update_particles(self.dt) From b8ab28b8f9e3146d4bbefda51075fa11cbfbecdb Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 09:12:14 -0400 Subject: [PATCH 05/27] fix typo --- pyro/compressible/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 78e3ec8ee..be57f3282 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -321,7 +321,7 @@ def evolve(self): xmom.v()[:, :] += 0.5*self.dt * \ ((dens.v() + U_old[:, :, self.ivars.ixmom])*grav + (ymom.v()**2 / dens.v() + - U_old[:, :, self.ivars.iymom]**2 / U_old[:. :, self.ivars.idens]) / myg.x2d.v()) - \ + U_old[:, :, self.ivars.iymom]**2 / U_old[:, :, self.ivars.idens]) / myg.x2d.v()) - \ self.dt * (qx.ip(1, n=self.ivars.ip) - qx.v(n=self.ivars.ip)) / myg.Lx.v() ymom.v()[:, :] += 0.5*self.dt * \ From 59a1921783990a3d570f6b5cd5eac467aba46b1b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 09:22:58 -0400 Subject: [PATCH 06/27] fix --- pyro/compressible/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index be57f3282..c518f6e5b 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -319,7 +319,7 @@ def evolve(self): if myg.coord_type == 1: xmom.v()[:, :] += 0.5*self.dt * \ - ((dens.v() + U_old[:, :, self.ivars.ixmom])*grav + + ((dens.v() + U_old[:, :, self.ivars.idens])*grav + (ymom.v()**2 / dens.v() + U_old[:, :, self.ivars.iymom]**2 / U_old[:, :, self.ivars.idens]) / myg.x2d.v()) - \ self.dt * (qx.ip(1, n=self.ivars.ip) - qx.v(n=self.ivars.ip)) / myg.Lx.v() From b476e6fc1d356ac10ebc2106d79be0130a5d8f7d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 10:04:24 -0400 Subject: [PATCH 07/27] fix stack --- pyro/compressible/simulation.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index c518f6e5b..48bfbb9d6 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -291,10 +291,11 @@ def evolve(self): self.ivars) # save the old state (without ghost cells) - U_old = np.stack([dens.v().copy(), - xmom.v().copy(), - ymom.v().copy(), - ener.v().copy()], axis=2) + U_old = myg.scratch_array(nvar=self.ivars.nvar) + U_old[:, :, self.ivars.idens] = dens[:, :] + U_old[:, :, self.ivars.ixmom] = xmom[:, :] + U_old[:, :, self.ivars.iymom] = ymom[:, :] + U_old[:, :, self.ivars.iener] = ener[:, :] # Conservative update @@ -319,21 +320,21 @@ def evolve(self): if myg.coord_type == 1: xmom.v()[:, :] += 0.5*self.dt * \ - ((dens.v() + U_old[:, :, self.ivars.idens])*grav + + ((dens.v() + U_old.v(n=self.ivars.idens))*grav + (ymom.v()**2 / dens.v() + - U_old[:, :, self.ivars.iymom]**2 / U_old[:, :, self.ivars.idens]) / myg.x2d.v()) - \ + U_old.v(n=self.ivars.iymom)**2 / U_old(n=self.ivars.idens)) / myg.x2d.v()) - \ self.dt * (qx.ip(1, n=self.ivars.ip) - qx.v(n=self.ivars.ip)) / myg.Lx.v() ymom.v()[:, :] += 0.5*self.dt * \ (-xmom.v()*ymom.v() / dens.v() - - U_old[:, :, self.ivars.ixmom] * U_old[:, :, self.ivars.iymom] / U_old[:, :, self.ivars.idens]) / myg.x2d.v() - \ + U_old.v(n=self.ivars.ixmom) * U_old.v(n=self.ivars.iymom) / U_old.v(n=self.ivars.idens)) / myg.x2d.v() - \ self.dt * (qy.jp(1, n=self.ivars.ip) - qy.v(n=self.ivars.ip)) / myg.Ly.v() - ener.v()[:, :] += 0.5 * self.dt * (xmom.v() + U_old[:, :, self.ivars.ixmom]) * grav + ener.v()[:, :] += 0.5 * self.dt * (xmom.v() + U_old.v(n=self.ivars.ixmom)) * grav else: - ymom.v()[:, :] += 0.5 * self.dt * (dens.v() + U_old[:, :, self.ivars.idens]) * grav - ener.v()[:, :] += 0.5 * self.dt * (ymom.v() + U_old[:, :, self.ivars.iymom]) * grav + ymom.v()[:, :] += 0.5 * self.dt * (dens.v() + U_old.v(n=self.ivars.idens)) * grav + ener.v()[:, :] += 0.5 * self.dt * (ymom.v() + U_old.v(n=self.ivars.iymom)) * grav if self.particles is not None: self.particles.update_particles(self.dt) From 110d7dc334a04598954a1fc8e8d57259771d42f5 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 10:48:30 -0400 Subject: [PATCH 08/27] do the corrector for the source --- pyro/compressible/simulation.py | 69 +++++++++++++++++++++-------- pyro/compressible/unsplit_fluxes.py | 2 +- 2 files changed, 51 insertions(+), 20 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 48bfbb9d6..6d7d0ab90 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -91,26 +91,42 @@ def prim_to_cons(q, gamma, ivars, myg): return U -def get_external_sources(t, U, ivars, rp, myg): +def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): """compute the external sources, including gravity""" _ = t # maybe unused - U_src = myg.scratch_array(nvar=ivars.nvar) + S = myg.scratch_array(nvar=ivars.nvar) grav = rp.get_param("compressible.grav") - if myg.coord_type == 1: - # gravity points in the radial direction for spherical - U_src[:, :, ivars.ixmom] = U[:, :, ivars.idens] * grav - U_src[:, :, ivars.iener] = U[:, :, ivars.ixmom] * grav + if U_old is None: + # we are just computing the sources from the current state U + + if myg.coord_type == 1: + # gravity points in the radial direction for spherical + S[:, :, ivars.ixmom] = U[:, :, ivars.idens] * grav + S[:, :, ivars.iener] = U[:, :, ivars.ixmom] * grav + + else: + # gravity points in the vertical (y) direction for Cartesian + S[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav + S[:, :, ivars.iener] = U[:, :, ivars.iymom] * grav else: - # gravity points in the vertical (y) direction for Cartesian - U_src[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav - U_src[:, :, ivars.iener] = U[:, :, ivars.iymom] * grav + # we want to compute gravity using the time-updated momentum + # we assume that U is an approximation to U^{n+1}, which includes + # a full dt * S_old + + S[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav + S_old_ymom = U_old[:, :, ivars.idens] * grav - return U_src + # we want the corrected ymom that has a time-centered source + ymom_new = U[:, :, ivars.iymom] + 0.5 * dt * (S[:, :, ivars.iymom] - S_old_ymom) + + S[:, :, ivars.iener] = ymom_new * grav + + return S class Simulation(NullSimulation): @@ -315,13 +331,9 @@ def evolve(self): # terms and the pressure gradient since we don't include # pressure in xmom and ymom fluxes - # For Cartesian2d (coord_type == 0): - # There is only gravity sources. - if myg.coord_type == 1: xmom.v()[:, :] += 0.5*self.dt * \ - ((dens.v() + U_old.v(n=self.ivars.idens))*grav + - (ymom.v()**2 / dens.v() + + ((ymom.v()**2 / dens.v() + U_old.v(n=self.ivars.iymom)**2 / U_old(n=self.ivars.idens)) / myg.x2d.v()) - \ self.dt * (qx.ip(1, n=self.ivars.ip) - qx.v(n=self.ivars.ip)) / myg.Lx.v() @@ -330,11 +342,30 @@ def evolve(self): U_old.v(n=self.ivars.ixmom) * U_old.v(n=self.ivars.iymom) / U_old.v(n=self.ivars.idens)) / myg.x2d.v() - \ self.dt * (qy.jp(1, n=self.ivars.ip) - qy.v(n=self.ivars.ip)) / myg.Ly.v() - ener.v()[:, :] += 0.5 * self.dt * (xmom.v() + U_old.v(n=self.ivars.ixmom)) * grav + # now the external sources (including gravity). We are going + # to do a predictor-corrector here: + # + # * compute old sources using old state: S^n = S(U^n) + # * update state full dt using old sources: U^{n+1,*} += dt * S^n + # * compute new sources using this updated state: S^{n+1) = S(U^{n+1,*}) + # * correct: U^{n+1} = U^{n+1,*} + dt/2 (S^{n+1} - S^n) - else: - ymom.v()[:, :] += 0.5 * self.dt * (dens.v() + U_old.v(n=self.ivars.idens)) * grav - ener.v()[:, :] += 0.5 * self.dt * (ymom.v() + U_old.v(n=self.ivars.iymom)) * grav + S_old = get_external_sources(self.cc_data.t, self.dt, U_old, + self.ivars, self.rp, myg) + + for n in range(self.ivars.nvar): + var = self.cc_data.get_var_by_index(n) + var.v()[:, :] += self.dt * S_old.v(n=n) + + # now get the new time source + + S_new = get_external_sources(self.cc_data.t, self.dt, self.cc_data.data, + self.ivars, self.rp, myg, U_old=U_old) + + # and correct + for n in range(self.ivars.nvar): + var = self.cc_data.get_var_by_index(n) + var.v()[:, :] += 0.5 * self.dt * (S_new.v(n=n) - S_old.v(n=n)) if self.particles is not None: self.particles.update_particles(self.dt) diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index 483b2d50a..4123e97b6 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -293,7 +293,7 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr, ymom_src = my_aux.get_var("ymom_src") E_src = my_aux.get_var("E_src") - U_src = comp.get_external_sources(my_data.t, my_data.data, ivars, rp, myg) + U_src = comp.get_external_sources(my_data.t, dt, my_data.data, ivars, rp, myg) dens_src[:, :] = U_src[:, :, ivars.idens] xmom_src[:, :] = U_src[:, :, ivars.ixmom] From 275bb467f28d674de86fd72827f85c843e9d03cc Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 11:01:01 -0400 Subject: [PATCH 09/27] fix pylint --- pyro/compressible/simulation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 6d7d0ab90..8d210900f 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -245,7 +245,6 @@ def evolve(self): ymom = self.cc_data.get_var("y-momentum") ener = self.cc_data.get_var("energy") - grav = self.rp.get_param("compressible.grav") gamma = self.rp.get_param("eos.gamma") myg = self.cc_data.grid From f5d6a8478e60931c646d2794eaf92e1ccd6b15bd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 11:38:58 -0400 Subject: [PATCH 10/27] this works now --- pyro/compressible/problems/heating.py | 19 +++++------ pyro/compressible/problems/inputs.heating | 41 +++++++++++++++++++++++ pyro/compressible/simulation.py | 13 +++++-- pyro/pyro_sim.py | 10 +++++- pyro/simulation_null.py | 6 +++- 5 files changed, 74 insertions(+), 15 deletions(-) create mode 100644 pyro/compressible/problems/inputs.heating diff --git a/pyro/compressible/problems/heating.py b/pyro/compressible/problems/heating.py index bb2d568f1..c70519669 100644 --- a/pyro/compressible/problems/heating.py +++ b/pyro/compressible/problems/heating.py @@ -7,7 +7,7 @@ from pyro.util import msg -DEFAULT_INPUTS = "inputs.sedov" +DEFAULT_INPUTS = "inputs.heating" PROBLEM_PARAMS = {"heating.rho_ambient": 1.0, # ambient density "heating.p_ambient": 10.0, # ambient pressure @@ -38,23 +38,22 @@ def init_data(my_data, rp): ener[:, :] = rp.get_param("heating.p_ambient") / (gamma - 1.0) -def source_terms(my_data, rp): +def source_terms(myg, U, ivars, rp): """source terms to be added to the evolution""" - e_src = my_data.get_var("E_src") + S = myg.scratch_array(nvar=ivars.nvar) - grid = my_data.grid + xctr = 0.5 * (myg.xmin + myg.xmax) + yctr = 0.5 * (myg.ymin + myg.ymax) - xctr = 0.5 * (grid.xmin + grid.xmax) - yctr = 0.5 * (grid.ymin + grid.ymax) - - dist = np.sqrt((my_data.grid.x2d - xctr)**2 + - (my_data.grid.y2d - yctr)**2) + dist = np.sqrt((myg.x2d - xctr)**2 + + (myg.y2d - yctr)**2) e_rate = rp.get_param("heating.e_rate") r_src = rp.get_param("heating.r_src") - e_src[:, :] = e_rate * np.exp(-(dist / r_src)**2) + S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_src)**2) + return S def finalize(): diff --git a/pyro/compressible/problems/inputs.heating b/pyro/compressible/problems/inputs.heating new file mode 100644 index 000000000..c8e8139ae --- /dev/null +++ b/pyro/compressible/problems/inputs.heating @@ -0,0 +1,41 @@ +[driver] +max_steps = 5000 +tmax = 1.0 + + +[compressible] +limiter = 2 +cvisc = 0.1 + + +[io] +basename = heating_ +dt_out = 0.1 + + +[eos] +gamma = 1.4 + + +[mesh] +nx = 64 +ny = 64 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[heating] +rho_ambient = 1.0 +p_ambient = 10.0 +r_src = 0.05 +e_rate = 0.1 + + +[vis] +dovis = 1 diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 8d210900f..6009329d9 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -91,7 +91,7 @@ def prim_to_cons(q, gamma, ivars, myg): return U -def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): +def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None, problem_source=None): """compute the external sources, including gravity""" _ = t # maybe unused @@ -126,6 +126,11 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): S[:, :, ivars.iener] = ymom_new * grav + # now add the heating + if problem_source: + S_heating = problem_source(myg, U, ivars, rp) + S[...] += S_heating + return S @@ -350,7 +355,8 @@ def evolve(self): # * correct: U^{n+1} = U^{n+1,*} + dt/2 (S^{n+1} - S^n) S_old = get_external_sources(self.cc_data.t, self.dt, U_old, - self.ivars, self.rp, myg) + self.ivars, self.rp, myg, + problem_source=self.problem_source) for n in range(self.ivars.nvar): var = self.cc_data.get_var_by_index(n) @@ -359,7 +365,8 @@ def evolve(self): # now get the new time source S_new = get_external_sources(self.cc_data.t, self.dt, self.cc_data.data, - self.ivars, self.rp, myg, U_old=U_old) + self.ivars, self.rp, myg, U_old=U_old, + problem_source=self.problem_source) # and correct for n in range(self.ivars.nvar): diff --git a/pyro/pyro_sim.py b/pyro/pyro_sim.py index c29c9163a..bdaf6efa9 100755 --- a/pyro/pyro_sim.py +++ b/pyro/pyro_sim.py @@ -69,6 +69,7 @@ def __init__(self, solver_name, *, from_commandline=False): self.problem_name = None self.problem_func = None + self.problem_source = None self.problem_params = None self.problem_finalize = None @@ -124,6 +125,7 @@ def initialize_problem(self, problem_name, *, inputs_file=None, inputs_dict=None self.problem_name = problem_name self.problem_func, self.problem_params = self.custom_problems[problem_name] self.problem_finalize = None + self.problem_source = None else: problem = importlib.import_module("pyro.{}.problems.{}".format(self.solver_name, problem_name)) @@ -131,6 +133,10 @@ def initialize_problem(self, problem_name, *, inputs_file=None, inputs_dict=None self.problem_func = problem.init_data self.problem_params = problem.PROBLEM_PARAMS self.problem_finalize = problem.finalize + try: + self.problem_source = problem.source_terms + except AttributeError: + self.problem_source = None if inputs_file is None: inputs_file = problem.DEFAULT_INPUTS @@ -175,7 +181,9 @@ def initialize_problem(self, problem_name, *, inputs_file=None, inputs_dict=None # are running self.sim = self.solver.Simulation( self.solver_name, self.problem_name, self.problem_func, self.rp, - problem_finalize_func=self.problem_finalize, timers=self.tc) + problem_finalize_func=self.problem_finalize, + problem_source_func=self.problem_source, + timers=self.tc) self.sim.initialize() self.sim.preevolve() diff --git a/pyro/simulation_null.py b/pyro/simulation_null.py index 81c86cc0d..c10fc75ce 100644 --- a/pyro/simulation_null.py +++ b/pyro/simulation_null.py @@ -101,7 +101,7 @@ def bc_setup(rp): class NullSimulation: def __init__(self, solver_name, problem_name, problem_func, rp, *, - problem_finalize_func=None, + problem_finalize_func=None, problem_source_func=None, timers=None, data_class=patch.CellCenterData2d): """ Initialize the Simulation object @@ -119,6 +119,9 @@ def __init__(self, solver_name, problem_name, problem_func, rp, *, problem_finalize_func : function An (optional) function to call when the simulation is over. + problem_source_func : function + An (optional) function to call to get source terms + for the state. timers : TimerCollection object, optional The timers used for profiling this simulation. data_class : CellCenterData2d or FV2d @@ -151,6 +154,7 @@ def __init__(self, solver_name, problem_name, problem_func, rp, *, self.problem_name = problem_name self.problem_func = problem_func self.problem_finalize = problem_finalize_func + self.problem_source = problem_source_func if timers is None: self.tc = profile.TimerCollection() From cdea66de2239898be7beb39a426d7225bbf1b60d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 11:57:27 -0400 Subject: [PATCH 11/27] update classes --- pyro/compressible_fv4/simulation.py | 4 +++- pyro/lm_atm/simulation.py | 7 +++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index b47d19821..6c085475c 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -6,9 +6,11 @@ class Simulation(compressible_rk.Simulation): def __init__(self, solver_name, problem_name, problem_func, rp, *, - problem_finalize_func=None, timers=None, data_class=fv.FV2d): + problem_finalize_func=None, problem_source_func=None, + timers=None, data_class=fv.FV2d): super().__init__(solver_name, problem_name, problem_func, rp, problem_finalize_func=problem_finalize_func, + problem_sourcE_func=problem_source_func, timers=timers, data_class=data_class) def substep(self, myd): diff --git a/pyro/lm_atm/simulation.py b/pyro/lm_atm/simulation.py index 91fa35438..04ac8e7d1 100644 --- a/pyro/lm_atm/simulation.py +++ b/pyro/lm_atm/simulation.py @@ -36,10 +36,13 @@ def jp(self, shift, buf=0): class Simulation(NullSimulation): def __init__(self, solver_name, problem_name, problem_func, rp, *, - problem_finalize_func=None, timers=None): + problem_finalize_func=None, problem_source_func=None, + timers=None): super().__init__(solver_name, problem_name, problem_func, rp, - problem_finalize_func=problem_finalize_func, timers=timers) + problem_finalize_func=problem_finalize_func, + problem_source_func=problem_source_func, + timers=timers) self.base = {} self.aux_data = None From 22cefacf343d97d6c8290ebbfa32e1d28e07e953 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 11:57:56 -0400 Subject: [PATCH 12/27] update modules --- pyro/compressible/problems/heating.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyro/compressible/problems/heating.py b/pyro/compressible/problems/heating.py index c70519669..a454ee77e 100644 --- a/pyro/compressible/problems/heating.py +++ b/pyro/compressible/problems/heating.py @@ -1,8 +1,6 @@ """A test of the energy sources. This uses a uniform domain and slowly adds heat to the center over time.""" -import math - import numpy as np from pyro.util import msg From 4352ef0834f12497bd2b0c427feb386864dc57e0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 12:01:16 -0400 Subject: [PATCH 13/27] update source --- pyro/compressible_fv4/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index 6c085475c..8efe4af59 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -10,7 +10,7 @@ def __init__(self, solver_name, problem_name, problem_func, rp, *, timers=None, data_class=fv.FV2d): super().__init__(solver_name, problem_name, problem_func, rp, problem_finalize_func=problem_finalize_func, - problem_sourcE_func=problem_source_func, + problem_source_func=problem_source_func, timers=timers, data_class=data_class) def substep(self, myd): From c8ef31be5e7f16f5b45ed55eb16ded38caf8ad7b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 13:51:41 -0400 Subject: [PATCH 14/27] add the plume problem --- pyro/compressible/problems/inputs.plume | 39 +++++++++++ pyro/compressible/problems/plume.py | 89 +++++++++++++++++++++++++ 2 files changed, 128 insertions(+) create mode 100644 pyro/compressible/problems/inputs.plume create mode 100644 pyro/compressible/problems/plume.py diff --git a/pyro/compressible/problems/inputs.plume b/pyro/compressible/problems/inputs.plume new file mode 100644 index 000000000..623452f98 --- /dev/null +++ b/pyro/compressible/problems/inputs.plume @@ -0,0 +1,39 @@ +# simple inputs files for the four-corner problem. + +[driver] +max_steps = 1000 +tmax = 100.0 + + +[io] +basename = bubble_ +n_out = 100 + + +[mesh] +nx = 128 +ny = 256 +xmax = 4.0 +ymax = 8.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = hse +yrboundary = hse + + +[plume] +scale_height = 3.0 +dens_base = 1000.0 + +x_pert = 2.0 +y_pert = 2.0 +r_pert = 0.25 +e_rate = 0.5 + + +[compressible] +grav = -2.0 + +limiter = 2 diff --git a/pyro/compressible/problems/plume.py b/pyro/compressible/problems/plume.py new file mode 100644 index 000000000..fa98f0d38 --- /dev/null +++ b/pyro/compressible/problems/plume.py @@ -0,0 +1,89 @@ +"""A heat source at a point creates a plume that buoynantly rises""" + +import numpy as np + +from pyro.util import msg + +DEFAULT_INPUTS = "inputs.plume" + +PROBLEM_PARAMS = {"plume.dens_base": 10.0, # density at the base of the atmosphere + "plume.scale_height": 4.0, # scale height of the isothermal atmosphere + "plume.x_pert": 2.0, + "plume.y_pert": 2.0, + "plume.r_pert": 0.25, + "plume.e_rate": 0.1, + "plume.dens_cutoff": 0.01} + + +def init_data(my_data, rp): + """ initialize the bubble problem """ + + if rp.get_param("driver.verbose"): + msg.bold("initializing the bubble problem...") + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + gamma = rp.get_param("eos.gamma") + + grav = rp.get_param("compressible.grav") + + scale_height = rp.get_param("plume.scale_height") + dens_base = rp.get_param("plume.dens_base") + dens_cutoff = rp.get_param("plume.dens_cutoff") + + # initialize the components, remember, that ener here is + # rho*eint + 0.5*rho*v**2, where eint is the specific + # internal energy (erg/g) + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + dens[:, :] = dens_cutoff + + # set the density to be stratified in the y-direction + myg = my_data.grid + + p = myg.scratch_array() + + pres_base = scale_height*dens_base*abs(grav) + + for j in range(myg.jlo, myg.jhi+1): + profile = 1.0 - (gamma-1.0)/gamma * myg.y[j]/scale_height + if profile > 0.0: + dens[:, j] = max(dens_base*(profile)**(1.0/(gamma-1.0)), + dens_cutoff) + else: + dens[:, j] = dens_cutoff + + if j == myg.jlo: + p[:, j] = pres_base + else: + p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*grav + + # set the energy (P = cs2*dens) + ener[:, :] = p[:, :]/(gamma - 1.0) + \ + 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :] + + +def source_terms(myg, U, ivars, rp): + """source terms to be added to the evolution""" + + S = myg.scratch_array(nvar=ivars.nvar) + + x_pert = rp.get_param("plume.x_pert") + y_pert = rp.get_param("plume.y_pert") + + dist = np.sqrt((myg.x2d - x_pert)**2 + + (myg.y2d - y_pert)**2) + + e_rate = rp.get_param("plume.e_rate") + r_pert = rp.get_param("plume.r_pert") + + S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_pert)**2) + return S + + +def finalize(): + """ print out any information to the user at the end of the run """ From 9bce42ac4932566fe751b170c9d6c2a8ebff8966 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 14:12:46 -0400 Subject: [PATCH 15/27] add more doc --- pyro/compressible/problems/plume.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pyro/compressible/problems/plume.py b/pyro/compressible/problems/plume.py index fa98f0d38..f37e8111d 100644 --- a/pyro/compressible/problems/plume.py +++ b/pyro/compressible/problems/plume.py @@ -1,4 +1,5 @@ -"""A heat source at a point creates a plume that buoynantly rises""" +"""A heat source at a point creates a plume that buoynantly rises in +an adiabatically stratified atmosphere.""" import numpy as np From 012a9e5c62468457e73532f4efdb7789523b8b45 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 14:27:01 -0400 Subject: [PATCH 16/27] fix spherical grav --- pyro/compressible/simulation.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 8d210900f..1d7493e25 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -118,13 +118,23 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): # we assume that U is an approximation to U^{n+1}, which includes # a full dt * S_old - S[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav - S_old_ymom = U_old[:, :, ivars.idens] * grav + if myg.coord_type == 1: + S[:, :, ivars.ixmom] = U[:, :, ivars.idens] * grav + S_old_xmom = U_old[:, :, ivars.idens] * grav + + # we want the corrected xmom that has a time-centered source + xmom_new = U[:, :, ivars.ixmom] + 0.5 * dt * (S[:, :, ivars.ixmom] - S_old_xmom) + + S[:, :, ivars.iener] = xmom_new * grav + + else: + S[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav + S_old_ymom = U_old[:, :, ivars.idens] * grav - # we want the corrected ymom that has a time-centered source - ymom_new = U[:, :, ivars.iymom] + 0.5 * dt * (S[:, :, ivars.iymom] - S_old_ymom) + # we want the corrected ymom that has a time-centered source + ymom_new = U[:, :, ivars.iymom] + 0.5 * dt * (S[:, :, ivars.iymom] - S_old_ymom) - S[:, :, ivars.iener] = ymom_new * grav + S[:, :, ivars.iener] = ymom_new * grav return S From 3ab90fe8fb4620c816e684306823a017e175abdd Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 14:43:21 -0400 Subject: [PATCH 17/27] update coord --- pyro/compressible/simulation.py | 26 ++++++++++++++------------ pyro/compressible/unsplit_fluxes.py | 8 +++----- 2 files changed, 17 insertions(+), 17 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 1d7493e25..ca7a8920d 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -100,6 +100,7 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): grav = rp.get_param("compressible.grav") + if U_old is None: # we are just computing the sources from the current state U @@ -108,6 +109,9 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): S[:, :, ivars.ixmom] = U[:, :, ivars.idens] * grav S[:, :, ivars.iener] = U[:, :, ivars.ixmom] * grav + S[:, :, ivars.ixmom] += U[:, :, ivars.iymom]**2 / (U[:, :, ivars.idens] * myg.x2d) + S[:, :, ivars.iymom] += -U[:, :, ivars.ixmom] * U[:, :, ivars.iymom] / U[:, :, ivars.idens] + else: # gravity points in the vertical (y) direction for Cartesian S[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav @@ -127,6 +131,9 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): S[:, :, ivars.iener] = xmom_new * grav + S[:, :, ivars.ixmom] += U[:, :, ivars.iymom]**2 / (U[:, :, ivars.idens] * myg.x2d) + S[:, :, ivars.iymom] += -U[:, :, ivars.ixmom] * U[:, :, ivars.iymom] / U[:, :, ivars.idens] + else: S[:, :, ivars.iymom] = U[:, :, ivars.idens] * grav S_old_ymom = U_old[:, :, ivars.idens] * grav @@ -336,20 +343,15 @@ def evolve(self): # Now apply external sources - # For SphericalPolar (coord_type == 1) there are geometric - # terms and the pressure gradient since we don't include - # pressure in xmom and ymom fluxes + # For SphericalPolar (coord_type == 1) there are pressure + # gradients since we don't include pressure in xmom and ymom + # fluxes if myg.coord_type == 1: - xmom.v()[:, :] += 0.5*self.dt * \ - ((ymom.v()**2 / dens.v() + - U_old.v(n=self.ivars.iymom)**2 / U_old(n=self.ivars.idens)) / myg.x2d.v()) - \ - self.dt * (qx.ip(1, n=self.ivars.ip) - qx.v(n=self.ivars.ip)) / myg.Lx.v() - - ymom.v()[:, :] += 0.5*self.dt * \ - (-xmom.v()*ymom.v() / dens.v() - - U_old.v(n=self.ivars.ixmom) * U_old.v(n=self.ivars.iymom) / U_old.v(n=self.ivars.idens)) / myg.x2d.v() - \ - self.dt * (qy.jp(1, n=self.ivars.ip) - qy.v(n=self.ivars.ip)) / myg.Ly.v() + xmom.v()[:, :] += self.dt * (qx.ip(1, n=self.ivars.ip) - + qx.v(n=self.ivars.ip)) / myg.Lx.v() + ymom.v()[:, :] += self.dt * (qy.jp(1, n=self.ivars.ip) - + qy.v(n=self.ivars.ip)) / myg.Ly.v() # now the external sources (including gravity). We are going # to do a predictor-corrector here: diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index 4123e97b6..48330fb27 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -300,12 +300,10 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr, ymom_src[:, :] = U_src[:, :, ivars.iymom] E_src[:, :] = U_src[:, :, ivars.iener] - # add geometric and pressure terms + # apply non-conservative pressure gradient if myg.coord_type == 1: - xmom_src.v()[:, :] += ymom.v()**2 / (dens.v()*myg.x2d.v()) - \ - (pres.ip(1) - pres.v()) / myg.Lx.v() - ymom_src.v()[:, :] += -(pres.jp(1) - pres.v()) / myg.Ly.v() - \ - xmom.v()*ymom.v() / (dens.v()*myg.x2d.v()) + xmom_src.v()[:, :] += -(pres.ip(1) - pres.v()) / myg.Lx.v() + ymom_src.v()[:, :] += -(pres.jp(1) - pres.v()) / myg.Ly.v() my_aux.fill_BC("dens_src") my_aux.fill_BC("xmom_src") From 71a022416f5a071340e0aef60d1e32518ce5df7c Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 14:46:57 -0400 Subject: [PATCH 18/27] fix flake8 --- pyro/compressible/simulation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index ca7a8920d..009b187d7 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -100,7 +100,6 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None): grav = rp.get_param("compressible.grav") - if U_old is None: # we are just computing the sources from the current state U From 26ffe6f90554ea94af37b432d35f8411432e9830 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 18 Oct 2024 14:48:04 -0400 Subject: [PATCH 19/27] fix pylint --- pyro/compressible/unsplit_fluxes.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index 48330fb27..e1a91033a 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -283,9 +283,6 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr, myg = my_data.grid - dens = my_data.get_var("density") - xmom = my_data.get_var("x-momentum") - ymom = my_data.get_var("y-momentum") pres = my_data.get_var("pressure") dens_src = my_aux.get_var("dens_src") From 0bde429630aeec32106f92bc9d02331d4e1d53de Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Mon, 21 Oct 2024 08:26:00 -0400 Subject: [PATCH 20/27] fix p gradient sign --- pyro/compressible/simulation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 009b187d7..69dd36448 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -347,9 +347,9 @@ def evolve(self): # fluxes if myg.coord_type == 1: - xmom.v()[:, :] += self.dt * (qx.ip(1, n=self.ivars.ip) - + xmom.v()[:, :] -= self.dt * (qx.ip(1, n=self.ivars.ip) - qx.v(n=self.ivars.ip)) / myg.Lx.v() - ymom.v()[:, :] += self.dt * (qy.jp(1, n=self.ivars.ip) - + ymom.v()[:, :] -= self.dt * (qy.jp(1, n=self.ivars.ip) - qy.v(n=self.ivars.ip)) / myg.Ly.v() # now the external sources (including gravity). We are going From 4ce25453b4a9f5d1324d4738d1a0b809acceb4bc Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 09:54:48 -0400 Subject: [PATCH 21/27] use the same compressible external source function for RK and SDC solvers --- pyro/compressible_fv4/simulation.py | 34 ++++++++++++++--------------- pyro/compressible_rk/simulation.py | 19 +++++----------- 2 files changed, 21 insertions(+), 32 deletions(-) diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index b47d19821..81785c2d9 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -1,5 +1,6 @@ import pyro.compressible_fv4.fluxes as flx from pyro import compressible_rk +from pyro.compressible import get_external_sources from pyro.mesh import fv, integration @@ -17,24 +18,23 @@ def substep(self, myd): """ myg = myd.grid - grav = self.rp.get_param("compressible.grav") - # compute the source terms -- we need to do these to 4th - # order. Start by evaluating the sources using the - # cell-center quantities (including one ghost cell. - dens_cc = myd.to_centers("density") - ymom_cc = myd.to_centers("y-momentum") + # compute the source terms -- we need to do this first + # using the cell-center data and then convert it back to + # averages - ymom_src = myg.scratch_array() - ymom_src.v(buf=1)[:, :] = dens_cc.v(buf=1)[:, :]*grav + U_cc = myg.scratch_array(nvar=self.ivars.nvar) + for n in range(self.ivars.nvar): + U_cc[:, :, n] = myd.to_centers(myg.get_var_by_index(n)) - E_src = myg.scratch_array() - E_src.v(buf=1)[:, :] = ymom_cc.v(buf=1)[:, :]*grav + # cell-centered sources + S = get_external_sources(self.cc_data.t, self.dt, U_cc, + self.ivars, self.rp, myg) - # now bring back to averages -- we only need this in the - # interior (no ghost cells) - ymom_src.v()[:, :] = ymom_src.v()[:, :] - myg.dx**2*ymom_src.lap()/24.0 - E_src.v()[:, :] = E_src.v()[:, :] - myg.dx**2*E_src.lap()/24.0 + # bring the sources back to averages -- we only care about + # the interior (no ghost cells) + for n in range(self.ivars.nvar): + S.v(n=n)[:, :] -= myg.dx**2 * S.lap(n=n) / 24.0 k = myg.scratch_array(nvar=self.ivars.nvar) @@ -43,10 +43,8 @@ def substep(self, myd): for n in range(self.ivars.nvar): k.v(n=n)[:, :] = \ (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \ - (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy - - k.v(n=self.ivars.iymom)[:, :] += ymom_src.v()[:, :] - k.v(n=self.ivars.iener)[:, :] += E_src.v()[:, :] + (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + \ + S.v(n=n) return k diff --git a/pyro/compressible_rk/simulation.py b/pyro/compressible_rk/simulation.py index 56be8d8f8..8fe5465bd 100644 --- a/pyro/compressible_rk/simulation.py +++ b/pyro/compressible_rk/simulation.py @@ -16,17 +16,10 @@ def substep(self, myd): """ myg = myd.grid - grav = self.rp.get_param("compressible.grav") - # compute the source terms - dens = myd.get_var("density") - ymom = myd.get_var("y-momentum") - - ymom_src = myg.scratch_array() - ymom_src.v()[:, :] = dens.v()[:, :]*grav - - E_src = myg.scratch_array() - E_src.v()[:, :] = ymom.v()[:, :]*grav + # source terms + S = compressible.get_external_sources(self.cc_data.t, self.dt, self.cc_data.data, + self.ivars, self.rp, myg) k = myg.scratch_array(nvar=self.ivars.nvar) @@ -36,10 +29,8 @@ def substep(self, myd): for n in range(self.ivars.nvar): k.v(n=n)[:, :] = \ (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \ - (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy - - k.v(n=self.ivars.iymom)[:, :] += ymom_src.v()[:, :] - k.v(n=self.ivars.iener)[:, :] += E_src.v()[:, :] + (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + \ + S.v(n=n) return k From f40525fca1742780552da1c5cc355bcc9e5b1cf0 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 09:58:05 -0400 Subject: [PATCH 22/27] fix flake8 --- pyro/compressible_fv4/simulation.py | 3 +-- pyro/compressible_rk/simulation.py | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index 81785c2d9..ba013cfbc 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -43,8 +43,7 @@ def substep(self, myd): for n in range(self.ivars.nvar): k.v(n=n)[:, :] = \ (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \ - (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + \ - S.v(n=n) + (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + S.v(n=n) return k diff --git a/pyro/compressible_rk/simulation.py b/pyro/compressible_rk/simulation.py index 8fe5465bd..e186abd0d 100644 --- a/pyro/compressible_rk/simulation.py +++ b/pyro/compressible_rk/simulation.py @@ -29,8 +29,7 @@ def substep(self, myd): for n in range(self.ivars.nvar): k.v(n=n)[:, :] = \ (flux_x.v(n=n) - flux_x.ip(1, n=n))/myg.dx + \ - (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + \ - S.v(n=n) + (flux_y.v(n=n) - flux_y.jp(1, n=n))/myg.dy + S.v(n=n) return k From e3e2ad7caf06fa2e5231cbea8f8d740233f3180d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 10:03:20 -0400 Subject: [PATCH 23/27] fix crash --- pyro/compressible_fv4/simulation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index ba013cfbc..53aa70adf 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -25,7 +25,7 @@ def substep(self, myd): U_cc = myg.scratch_array(nvar=self.ivars.nvar) for n in range(self.ivars.nvar): - U_cc[:, :, n] = myd.to_centers(myg.get_var_by_index(n)) + U_cc[:, :, n] = myd.to_centers(self.cc_data.names[n]) # cell-centered sources S = get_external_sources(self.cc_data.t, self.dt, U_cc, From f3cbca1201a63ce26cdb66f7348fefb2716fd3d3 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 10:22:57 -0400 Subject: [PATCH 24/27] fix stage data --- pyro/compressible_fv4/simulation.py | 4 ++-- pyro/compressible_rk/simulation.py | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index 53aa70adf..5f096be68 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -25,10 +25,10 @@ def substep(self, myd): U_cc = myg.scratch_array(nvar=self.ivars.nvar) for n in range(self.ivars.nvar): - U_cc[:, :, n] = myd.to_centers(self.cc_data.names[n]) + U_cc[:, :, n] = myd.to_centers(myd.names[n]) # cell-centered sources - S = get_external_sources(self.cc_data.t, self.dt, U_cc, + S = get_external_sources(myd.t, self.dt, U_cc, self.ivars, self.rp, myg) # bring the sources back to averages -- we only care about diff --git a/pyro/compressible_rk/simulation.py b/pyro/compressible_rk/simulation.py index e186abd0d..a4d5aad35 100644 --- a/pyro/compressible_rk/simulation.py +++ b/pyro/compressible_rk/simulation.py @@ -17,8 +17,9 @@ def substep(self, myd): myg = myd.grid - # source terms - S = compressible.get_external_sources(self.cc_data.t, self.dt, self.cc_data.data, + # source terms -- note: this dt is the entire dt, not the + # stage's dt + S = compressible.get_external_sources(myd.t, self.dt, myd.data, self.ivars, self.rp, myg) k = myg.scratch_array(nvar=self.ivars.nvar) From 87d148fd50b3e23cd73af64b8176526f2ef07f03 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 11:05:41 -0400 Subject: [PATCH 25/27] hook problem sources into RK --- pyro/compressible/problems/inputs.plume | 6 +++--- pyro/compressible_fv4/simulation.py | 3 ++- pyro/compressible_rk/simulation.py | 3 ++- 3 files changed, 7 insertions(+), 5 deletions(-) diff --git a/pyro/compressible/problems/inputs.plume b/pyro/compressible/problems/inputs.plume index 623452f98..d35725942 100644 --- a/pyro/compressible/problems/inputs.plume +++ b/pyro/compressible/problems/inputs.plume @@ -1,12 +1,12 @@ # simple inputs files for the four-corner problem. [driver] -max_steps = 1000 -tmax = 100.0 +max_steps = 10000 +tmax = 10.0 [io] -basename = bubble_ +basename = plume_ n_out = 100 diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index ea9202bf7..df662ae4a 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -31,7 +31,8 @@ def substep(self, myd): # cell-centered sources S = get_external_sources(myd.t, self.dt, U_cc, - self.ivars, self.rp, myg) + self.ivars, self.rp, myg, + problem_source=self.problem_source) # bring the sources back to averages -- we only care about # the interior (no ghost cells) diff --git a/pyro/compressible_rk/simulation.py b/pyro/compressible_rk/simulation.py index a4d5aad35..f0306403a 100644 --- a/pyro/compressible_rk/simulation.py +++ b/pyro/compressible_rk/simulation.py @@ -20,7 +20,8 @@ def substep(self, myd): # source terms -- note: this dt is the entire dt, not the # stage's dt S = compressible.get_external_sources(myd.t, self.dt, myd.data, - self.ivars, self.rp, myg) + self.ivars, self.rp, myg, + problem_source=self.problem_source) k = myg.scratch_array(nvar=self.ivars.nvar) From d4a1486ae195a0a3d750cb1feb7da6842723d1a1 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 11:59:57 -0400 Subject: [PATCH 26/27] add source to the predictor --- pyro/compressible/simulation.py | 3 ++- pyro/compressible/unsplit_fluxes.py | 8 ++++++-- 2 files changed, 8 insertions(+), 3 deletions(-) diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index fb5977742..89be94994 100644 --- a/pyro/compressible/simulation.py +++ b/pyro/compressible/simulation.py @@ -279,7 +279,8 @@ def evolve(self): # Only gravitional source for Cartesian2d U_xl, U_xr, U_yl, U_yr = flx.apply_source_terms(U_xl, U_xr, U_yl, U_yr, self.cc_data, self.aux_data, self.rp, - self.ivars, self.tc, self.dt) + self.ivars, self.tc, self.dt, + problem_source=self.problem_source) # Apply transverse corrections. U_xl, U_xr, U_yl, U_yr = flx.apply_transverse_flux(U_xl, U_xr, U_yl, U_yr, diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index e1a91033a..64c656a9b 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -243,7 +243,8 @@ def interface_states(my_data, rp, ivars, tc, dt): def apply_source_terms(U_xl, U_xr, U_yl, U_yr, - my_data, my_aux, rp, ivars, tc, dt): + my_data, my_aux, rp, ivars, tc, dt, *, + problem_source=None): """ This function applies source terms including external (gravity), geometric terms, and pressure terms to the left and right @@ -270,6 +271,8 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr, The timers we are using to profile dt : float The timestep we are advancing through. + problem_source : function (optional) + A problem-specific source function to call Returns ------- @@ -290,7 +293,8 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr, ymom_src = my_aux.get_var("ymom_src") E_src = my_aux.get_var("E_src") - U_src = comp.get_external_sources(my_data.t, dt, my_data.data, ivars, rp, myg) + U_src = comp.get_external_sources(my_data.t, dt, my_data.data, ivars, rp, myg, + problem_source=problem_source) dens_src[:, :] = U_src[:, :, ivars.idens] xmom_src[:, :] = U_src[:, :, ivars.ixmom] From 5ef519413ef761bc13a19c32cab419b3031b2e1b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 1 Nov 2024 13:44:30 -0400 Subject: [PATCH 27/27] add a convection problem --- pyro/compressible/problems/convection.py | 104 +++++++++++++++++++ pyro/compressible/problems/inputs.convection | 39 +++++++ 2 files changed, 143 insertions(+) create mode 100644 pyro/compressible/problems/convection.py create mode 100644 pyro/compressible/problems/inputs.convection diff --git a/pyro/compressible/problems/convection.py b/pyro/compressible/problems/convection.py new file mode 100644 index 000000000..f5816b320 --- /dev/null +++ b/pyro/compressible/problems/convection.py @@ -0,0 +1,104 @@ +"""A heat source in a layer at some height above the bottom will drive +convection in an adiabatically stratified atmosphere.""" + +import numpy as np + +from pyro.util import msg + +DEFAULT_INPUTS = "inputs.convection" + +PROBLEM_PARAMS = {"convection.dens_base": 10.0, # density at the base of the atmosphere + "convection.scale_height": 4.0, # scale height of the isothermal atmosphere + "convection.y_height": 2.0, + "convection.thickness": 0.25, + "convection.e_rate": 0.1, + "convection.dens_cutoff": 0.01} + + +def init_data(my_data, rp): + """ initialize the bubble problem """ + + if rp.get_param("driver.verbose"): + msg.bold("initializing the bubble problem...") + + # get the density, momenta, and energy as separate variables + dens = my_data.get_var("density") + xmom = my_data.get_var("x-momentum") + ymom = my_data.get_var("y-momentum") + ener = my_data.get_var("energy") + + gamma = rp.get_param("eos.gamma") + + grav = rp.get_param("compressible.grav") + + scale_height = rp.get_param("convection.scale_height") + dens_base = rp.get_param("convection.dens_base") + dens_cutoff = rp.get_param("convection.dens_cutoff") + + # initialize the components, remember, that ener here is + # rho*eint + 0.5*rho*v**2, where eint is the specific + # internal energy (erg/g) + xmom[:, :] = 0.0 + ymom[:, :] = 0.0 + dens[:, :] = dens_cutoff + + # set the density to be stratified in the y-direction + myg = my_data.grid + + p = myg.scratch_array() + + pres_base = scale_height*dens_base*abs(grav) + + for j in range(myg.jlo, myg.jhi+1): + profile = 1.0 - (gamma-1.0)/gamma * myg.y[j]/scale_height + if profile > 0.0: + dens[:, j] = max(dens_base*(profile)**(1.0/(gamma-1.0)), + dens_cutoff) + else: + dens[:, j] = dens_cutoff + + if j == myg.jlo: + p[:, j] = pres_base + elif dens[0, j] <= dens_cutoff + 1.e-30: + p[:, j] = p[:, j-1] + else: + p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*grav + + # set the energy (P = cs2*dens) + ener[:, :] = p[:, :]/(gamma - 1.0) + \ + 0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :] + + # pairs of random numbers between [-1, 1] + vel_pert = 2.0 * np.random.random_sample((myg.qx, myg.qy, 2)) - 1 + + cs = np.sqrt(gamma * p / dens) + + # make vel_pert have M < 0.05 + vel_pert[:, :, 0] *= 0.05 * cs + vel_pert[:, :, 1] *= 0.05 * cs + + idx = dens > 2 * dens_cutoff + xmom[idx] = dens[idx] * vel_pert[idx, 0] + ymom[idx] = dens[idx] * vel_pert[idx, 1] + + ener[:, :] += 0.5 * (xmom[:, :]**2 + ymom[:, :]**2) / dens[:, :] + + +def source_terms(myg, U, ivars, rp): + """source terms to be added to the evolution""" + + S = myg.scratch_array(nvar=ivars.nvar) + + y_height = rp.get_param("convection.y_height") + + dist = np.abs(myg.y2d - y_height) + + e_rate = rp.get_param("convection.e_rate") + thick = rp.get_param("convection.thickness") + + S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / thick)**2) + return S + + +def finalize(): + """ print out any information to the user at the end of the run """ diff --git a/pyro/compressible/problems/inputs.convection b/pyro/compressible/problems/inputs.convection new file mode 100644 index 000000000..36389589e --- /dev/null +++ b/pyro/compressible/problems/inputs.convection @@ -0,0 +1,39 @@ +# simple inputs files for the four-corner problem. + +[driver] +max_steps = 10000 +tmax = 10.0 + + +[io] +basename = convection_ +n_out = 100 + + +[mesh] +nx = 128 +ny = 256 +xmax = 4.0 +ymax = 8.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = reflect +yrboundary = outflow + + +[convection] +scale_height = 2.0 +dens_base = 1000.0 + +x_pert = 2.0 +y_pert = 2.0 +r_pert = 0.25 +e_rate = 0.5 + + +[compressible] +grav = -2.0 + +limiter = 2