Skip to content

Commit 3cda63b

Browse files
authored
compressible: add the ability for a problem-dependent external source (#289)
this also adds a "heating" test problem.
1 parent 13eb252 commit 3cda63b

File tree

11 files changed

+279
-13
lines changed

11 files changed

+279
-13
lines changed

pyro/compressible/problems/heating.py

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
"""A test of the energy sources. This uses a uniform domain and
2+
slowly adds heat to the center over time."""
3+
4+
import numpy as np
5+
6+
from pyro.util import msg
7+
8+
DEFAULT_INPUTS = "inputs.heating"
9+
10+
PROBLEM_PARAMS = {"heating.rho_ambient": 1.0, # ambient density
11+
"heating.p_ambient": 10.0, # ambient pressure
12+
"heating.r_src": 0.1, # physical size of the heating src
13+
"heating.e_rate": 0.1} # energy generation rate (energy / mass / time)
14+
15+
16+
def init_data(my_data, rp):
17+
""" initialize the heating problem """
18+
19+
if rp.get_param("driver.verbose"):
20+
msg.bold("initializing the heating problem...")
21+
22+
# get the density, momenta, and energy as separate variables
23+
dens = my_data.get_var("density")
24+
xmom = my_data.get_var("x-momentum")
25+
ymom = my_data.get_var("y-momentum")
26+
ener = my_data.get_var("energy")
27+
28+
gamma = rp.get_param("eos.gamma")
29+
30+
# initialize the components, remember, that ener here is rho*eint
31+
# + 0.5*rho*v**2, where eint is the specific internal energy
32+
# (erg/g)
33+
dens[:, :] = rp.get_param("heating.rho_ambient")
34+
xmom[:, :] = 0.0
35+
ymom[:, :] = 0.0
36+
ener[:, :] = rp.get_param("heating.p_ambient") / (gamma - 1.0)
37+
38+
39+
def source_terms(myg, U, ivars, rp):
40+
"""source terms to be added to the evolution"""
41+
42+
S = myg.scratch_array(nvar=ivars.nvar)
43+
44+
xctr = 0.5 * (myg.xmin + myg.xmax)
45+
yctr = 0.5 * (myg.ymin + myg.ymax)
46+
47+
dist = np.sqrt((myg.x2d - xctr)**2 +
48+
(myg.y2d - yctr)**2)
49+
50+
e_rate = rp.get_param("heating.e_rate")
51+
r_src = rp.get_param("heating.r_src")
52+
53+
S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_src)**2)
54+
return S
55+
56+
57+
def finalize():
58+
""" print out any information to the user at the end of the run """
59+
60+
print("""
61+
The script analysis/sedov_compare.py can be used to analyze these
62+
results. That will perform an average at constant radius and
63+
compare the radial profiles to the exact solution. Sample exact
64+
data is provided as analysis/cylindrical-sedov.out
65+
""")
Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
[driver]
2+
max_steps = 5000
3+
tmax = 1.0
4+
5+
6+
[compressible]
7+
limiter = 2
8+
cvisc = 0.1
9+
10+
11+
[io]
12+
basename = heating_
13+
dt_out = 0.1
14+
15+
16+
[eos]
17+
gamma = 1.4
18+
19+
20+
[mesh]
21+
nx = 64
22+
ny = 64
23+
xmax = 1.0
24+
ymax = 1.0
25+
26+
xlboundary = outflow
27+
xrboundary = outflow
28+
29+
ylboundary = outflow
30+
yrboundary = outflow
31+
32+
33+
[heating]
34+
rho_ambient = 1.0
35+
p_ambient = 10.0
36+
r_src = 0.05
37+
e_rate = 0.1
38+
39+
40+
[vis]
41+
dovis = 1
Lines changed: 39 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
# simple inputs files for the four-corner problem.
2+
3+
[driver]
4+
max_steps = 10000
5+
tmax = 10.0
6+
7+
8+
[io]
9+
basename = plume_
10+
n_out = 100
11+
12+
13+
[mesh]
14+
nx = 128
15+
ny = 256
16+
xmax = 4.0
17+
ymax = 8.0
18+
19+
xlboundary = outflow
20+
xrboundary = outflow
21+
22+
ylboundary = hse
23+
yrboundary = hse
24+
25+
26+
[plume]
27+
scale_height = 3.0
28+
dens_base = 1000.0
29+
30+
x_pert = 2.0
31+
y_pert = 2.0
32+
r_pert = 0.25
33+
e_rate = 0.5
34+
35+
36+
[compressible]
37+
grav = -2.0
38+
39+
limiter = 2

pyro/compressible/problems/plume.py

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
"""A heat source at a point creates a plume that buoynantly rises in
2+
an adiabatically stratified atmosphere."""
3+
4+
import numpy as np
5+
6+
from pyro.util import msg
7+
8+
DEFAULT_INPUTS = "inputs.plume"
9+
10+
PROBLEM_PARAMS = {"plume.dens_base": 10.0, # density at the base of the atmosphere
11+
"plume.scale_height": 4.0, # scale height of the isothermal atmosphere
12+
"plume.x_pert": 2.0,
13+
"plume.y_pert": 2.0,
14+
"plume.r_pert": 0.25,
15+
"plume.e_rate": 0.1,
16+
"plume.dens_cutoff": 0.01}
17+
18+
19+
def init_data(my_data, rp):
20+
""" initialize the bubble problem """
21+
22+
if rp.get_param("driver.verbose"):
23+
msg.bold("initializing the bubble problem...")
24+
25+
# get the density, momenta, and energy as separate variables
26+
dens = my_data.get_var("density")
27+
xmom = my_data.get_var("x-momentum")
28+
ymom = my_data.get_var("y-momentum")
29+
ener = my_data.get_var("energy")
30+
31+
gamma = rp.get_param("eos.gamma")
32+
33+
grav = rp.get_param("compressible.grav")
34+
35+
scale_height = rp.get_param("plume.scale_height")
36+
dens_base = rp.get_param("plume.dens_base")
37+
dens_cutoff = rp.get_param("plume.dens_cutoff")
38+
39+
# initialize the components, remember, that ener here is
40+
# rho*eint + 0.5*rho*v**2, where eint is the specific
41+
# internal energy (erg/g)
42+
xmom[:, :] = 0.0
43+
ymom[:, :] = 0.0
44+
dens[:, :] = dens_cutoff
45+
46+
# set the density to be stratified in the y-direction
47+
myg = my_data.grid
48+
49+
p = myg.scratch_array()
50+
51+
pres_base = scale_height*dens_base*abs(grav)
52+
53+
for j in range(myg.jlo, myg.jhi+1):
54+
profile = 1.0 - (gamma-1.0)/gamma * myg.y[j]/scale_height
55+
if profile > 0.0:
56+
dens[:, j] = max(dens_base*(profile)**(1.0/(gamma-1.0)),
57+
dens_cutoff)
58+
else:
59+
dens[:, j] = dens_cutoff
60+
61+
if j == myg.jlo:
62+
p[:, j] = pres_base
63+
else:
64+
p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*grav
65+
66+
# set the energy (P = cs2*dens)
67+
ener[:, :] = p[:, :]/(gamma - 1.0) + \
68+
0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]
69+
70+
71+
def source_terms(myg, U, ivars, rp):
72+
"""source terms to be added to the evolution"""
73+
74+
S = myg.scratch_array(nvar=ivars.nvar)
75+
76+
x_pert = rp.get_param("plume.x_pert")
77+
y_pert = rp.get_param("plume.y_pert")
78+
79+
dist = np.sqrt((myg.x2d - x_pert)**2 +
80+
(myg.y2d - y_pert)**2)
81+
82+
e_rate = rp.get_param("plume.e_rate")
83+
r_pert = rp.get_param("plume.r_pert")
84+
85+
S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_pert)**2)
86+
return S
87+
88+
89+
def finalize():
90+
""" print out any information to the user at the end of the run """

pyro/compressible/simulation.py

Lines changed: 12 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -91,7 +91,7 @@ def prim_to_cons(q, gamma, ivars, myg):
9191
return U
9292

9393

94-
def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None):
94+
def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None, problem_source=None):
9595
"""compute the external sources, including gravity"""
9696

9797
_ = t # maybe unused
@@ -142,6 +142,11 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None):
142142

143143
S[:, :, ivars.iener] = ymom_new * grav
144144

145+
# now add the heating
146+
if problem_source:
147+
S_heating = problem_source(myg, U, ivars, rp)
148+
S[...] += S_heating
149+
145150
return S
146151

147152

@@ -274,7 +279,8 @@ def evolve(self):
274279
# Only gravitional source for Cartesian2d
275280
U_xl, U_xr, U_yl, U_yr = flx.apply_source_terms(U_xl, U_xr, U_yl, U_yr,
276281
self.cc_data, self.aux_data, self.rp,
277-
self.ivars, self.tc, self.dt)
282+
self.ivars, self.tc, self.dt,
283+
problem_source=self.problem_source)
278284

279285
# Apply transverse corrections.
280286
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):
361367
# * correct: U^{n+1} = U^{n+1,*} + dt/2 (S^{n+1} - S^n)
362368

363369
S_old = get_external_sources(self.cc_data.t, self.dt, U_old,
364-
self.ivars, self.rp, myg)
370+
self.ivars, self.rp, myg,
371+
problem_source=self.problem_source)
365372

366373
for n in range(self.ivars.nvar):
367374
var = self.cc_data.get_var_by_index(n)
@@ -370,7 +377,8 @@ def evolve(self):
370377
# now get the new time source
371378

372379
S_new = get_external_sources(self.cc_data.t, self.dt, self.cc_data.data,
373-
self.ivars, self.rp, myg, U_old=U_old)
380+
self.ivars, self.rp, myg, U_old=U_old,
381+
problem_source=self.problem_source)
374382

375383
# and correct
376384
for n in range(self.ivars.nvar):

pyro/compressible/unsplit_fluxes.py

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -243,7 +243,8 @@ def interface_states(my_data, rp, ivars, tc, dt):
243243

244244

245245
def apply_source_terms(U_xl, U_xr, U_yl, U_yr,
246-
my_data, my_aux, rp, ivars, tc, dt):
246+
my_data, my_aux, rp, ivars, tc, dt, *,
247+
problem_source=None):
247248
"""
248249
This function applies source terms including external (gravity),
249250
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,
270271
The timers we are using to profile
271272
dt : float
272273
The timestep we are advancing through.
274+
problem_source : function (optional)
275+
A problem-specific source function to call
273276
274277
Returns
275278
-------
@@ -290,7 +293,8 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr,
290293
ymom_src = my_aux.get_var("ymom_src")
291294
E_src = my_aux.get_var("E_src")
292295

293-
U_src = comp.get_external_sources(my_data.t, dt, my_data.data, ivars, rp, myg)
296+
U_src = comp.get_external_sources(my_data.t, dt, my_data.data, ivars, rp, myg,
297+
problem_source=problem_source)
294298

295299
dens_src[:, :] = U_src[:, :, ivars.idens]
296300
xmom_src[:, :] = U_src[:, :, ivars.ixmom]

pyro/compressible_fv4/simulation.py

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,11 @@
77
class Simulation(compressible_rk.Simulation):
88

99
def __init__(self, solver_name, problem_name, problem_func, rp, *,
10-
problem_finalize_func=None, timers=None, data_class=fv.FV2d):
10+
problem_finalize_func=None, problem_source_func=None,
11+
timers=None, data_class=fv.FV2d):
1112
super().__init__(solver_name, problem_name, problem_func, rp,
1213
problem_finalize_func=problem_finalize_func,
14+
problem_source_func=problem_source_func,
1315
timers=timers, data_class=data_class)
1416

1517
def substep(self, myd):
@@ -29,7 +31,8 @@ def substep(self, myd):
2931

3032
# cell-centered sources
3133
S = get_external_sources(myd.t, self.dt, U_cc,
32-
self.ivars, self.rp, myg)
34+
self.ivars, self.rp, myg,
35+
problem_source=self.problem_source)
3336

3437
# bring the sources back to averages -- we only care about
3538
# the interior (no ghost cells)

pyro/compressible_rk/simulation.py

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -20,7 +20,8 @@ def substep(self, myd):
2020
# source terms -- note: this dt is the entire dt, not the
2121
# stage's dt
2222
S = compressible.get_external_sources(myd.t, self.dt, myd.data,
23-
self.ivars, self.rp, myg)
23+
self.ivars, self.rp, myg,
24+
problem_source=self.problem_source)
2425

2526
k = myg.scratch_array(nvar=self.ivars.nvar)
2627

pyro/lm_atm/simulation.py

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -36,10 +36,13 @@ def jp(self, shift, buf=0):
3636
class Simulation(NullSimulation):
3737

3838
def __init__(self, solver_name, problem_name, problem_func, rp, *,
39-
problem_finalize_func=None, timers=None):
39+
problem_finalize_func=None, problem_source_func=None,
40+
timers=None):
4041

4142
super().__init__(solver_name, problem_name, problem_func, rp,
42-
problem_finalize_func=problem_finalize_func, timers=timers)
43+
problem_finalize_func=problem_finalize_func,
44+
problem_source_func=problem_source_func,
45+
timers=timers)
4346

4447
self.base = {}
4548
self.aux_data = None

0 commit comments

Comments
 (0)