diff --git a/pyro/compressible/problems/heating.py b/pyro/compressible/problems/heating.py new file mode 100644 index 000000000..a454ee77e --- /dev/null +++ b/pyro/compressible/problems/heating.py @@ -0,0 +1,65 @@ +"""A test of the energy sources. This uses a uniform domain and +slowly adds heat to the center over time.""" + +import numpy as np + +from pyro.util import msg + +DEFAULT_INPUTS = "inputs.heating" + +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(myg, U, ivars, rp): + """source terms to be added to the evolution""" + + S = myg.scratch_array(nvar=ivars.nvar) + + xctr = 0.5 * (myg.xmin + myg.xmax) + yctr = 0.5 * (myg.ymin + myg.ymax) + + 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") + + S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_src)**2) + return S + + +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 + """) 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/problems/inputs.plume b/pyro/compressible/problems/inputs.plume new file mode 100644 index 000000000..d35725942 --- /dev/null +++ b/pyro/compressible/problems/inputs.plume @@ -0,0 +1,39 @@ +# simple inputs files for the four-corner problem. + +[driver] +max_steps = 10000 +tmax = 10.0 + + +[io] +basename = plume_ +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..f37e8111d --- /dev/null +++ b/pyro/compressible/problems/plume.py @@ -0,0 +1,90 @@ +"""A heat source at a point creates a plume that buoynantly rises in +an adiabatically stratified atmosphere.""" + +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 """ diff --git a/pyro/compressible/simulation.py b/pyro/compressible/simulation.py index 69dd36448..89be94994 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 @@ -142,6 +142,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 @@ -274,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, @@ -361,7 +367,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) @@ -370,7 +377,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/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] diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index 5f096be68..df662ae4a 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -7,9 +7,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): @@ -29,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) 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 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()