Skip to content

Commit 279776b

Browse files
authored
viscous burgers solver (#171)
this adds a viscous Burgers' solver by adding parabolic solves for the viscous terms
1 parent 113a50d commit 279776b

25 files changed

+683
-1
lines changed

README.md

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -134,6 +134,10 @@ pyro provides the following solvers (all in 2-d):
134134
135135
- `burgers`: a second-order unsplit solver for invsicid Burgers' equation.
136136
137+
- `viscous_burgers`: a second-order unsplit solver for viscous Burgers' equation
138+
with constant-coefficient diffusion. It uses Crank-Nicolson time-discretized
139+
solver for solving diffusion.
140+
137141
- `compressible`: a second-order unsplit solver for the Euler
138142
equations of compressible hydrodynamics. This uses characteristic
139143
tracing and corner coupling for the prediction of the interface

docs/source/burgers_basics.rst

Lines changed: 28 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -6,7 +6,7 @@ Burgers' Equation is a nonlinear hyperbolic equation. It has the same form as th
66
``Inviscid Burgers``
77
--------------------------------
88

9-
A 2D Burgers' Equation has the following form:
9+
A 2D inviscid Burgers' Equation has the following form:
1010

1111
.. math::
1212
@@ -34,3 +34,30 @@ The parameters for this solver are:
3434
:align: center
3535

3636
The figure above is generated using ``burgers/problems/test.py``, which is used to test the validity of the solver. Bottom-left of the domain has a higher velocity than the top-right domain. With :math:`u_{i,j}=v_{i,j}`, the wave travels diagonally to the top-right with a constant velocity that is equal to the shock speed. ``burgers/problem/verify.py`` can be used to calculate the wave speed using outputs from ``test.py`` and compare to the theoretical shock speed.
37+
38+
``Viscous Burgers``
39+
--------------------------------
40+
41+
A 2D viscous Burgers' Equation has the following form:
42+
43+
.. math::
44+
45+
u_t + u u_x + v u_y = \epsilon \left( u_{xx} + u_{yy}\right) \\
46+
v_t + u v_x + v v_y = \epsilon \left( v_{xx} + v_{yy}\right)
47+
48+
The viscous Burgers' equation has an additional velocity diffusion term on the RHS compared to the inviscid Burgers' equation. Here :math:`\epsilon` represents the constant viscosity.
49+
50+
:py:mod:`pyro.viscous_burgers` is inherited from :py:mod:`pyro.burgers`, where we added an additional diffusion term when constructing the interface states. We then solve for diffusion along with the extra advective source to the Helmholtz equation by using the Crank-Nicolson discretization and multigrid solvers.
51+
52+
53+
The parameters for this solver are:
54+
55+
.. include:: viscous_burgers_defaults.inc
56+
57+
58+
.. image:: viscous_burgers.png
59+
:align: center
60+
61+
The figure above is generated using ``viscous_burgers/problems/test.py``, which has the identical setup as in ``burgers/problems/test.py``. With diffusion added to the system, we see the shock (discontinuity) is smeared out as system evolves.
62+
63+

docs/source/pyro.rst

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,7 @@ Subpackages
1818
pyro.advection_rk
1919
pyro.advection_weno
2020
pyro.burgers
21+
pyro.viscous_burgers
2122
pyro.compressible
2223
pyro.compressible_fv4
2324
pyro.compressible_react
Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
pyro.viscous\_burgers.problems package
2+
=======================================
3+
4+
.. automodule:: pyro.viscous_burgers.problems
5+
:members:
6+
:undoc-members:
7+
:show-inheritance:
8+
9+
Submodules
10+
----------
11+
12+
pyro.viscous\_burgers.problems.converge module
13+
-----------------------------------------------
14+
15+
.. automodule:: pyro.viscous_burgers.problems.converge
16+
:members:
17+
:undoc-members:
18+
:show-inheritance:
19+
20+
pyro.viscous\_burgers.problems.tophat module
21+
--------------------------------------------
22+
23+
.. automodule:: pyro.viscous_burgers.problems.tophat
24+
:members:
25+
:undoc-members:
26+
:show-inheritance:

docs/source/pyro.viscous_burgers.rst

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
pyro.viscous\_burgers package
2+
==============================
3+
4+
.. automodule:: pyro.viscous_burgers
5+
:members:
6+
:undoc-members:
7+
:show-inheritance:
8+
9+
Subpackages
10+
-----------
11+
12+
.. toctree::
13+
:maxdepth: 4
14+
15+
pyro.viscous_burgers.problems
16+
17+
Submodules
18+
----------
19+
20+
pyro.viscous\_burgers.interface module
21+
---------------------------------------
22+
23+
.. automodule:: pyro.viscous_burgers.interface
24+
:members:
25+
:undoc-members:
26+
:show-inheritance:
27+
28+
pyro.viscous_burgers.simulation module
29+
---------------------------------------
30+
31+
.. automodule:: pyro.viscous_burgers.simulation
32+
:members:
33+
:undoc-members:
34+
:show-inheritance:
Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
* section: [advection]
2+
3+
+----------------------------------+----------------+----------------------------------------------------+
4+
| option | value | description |
5+
+==================================+================+====================================================+
6+
| limiter | ``2`` | limiter (0 = none, 1 = 2nd order, 2 = 4th order) |
7+
+----------------------------------+----------------+----------------------------------------------------+
8+
9+
* section: [diffusion]
10+
11+
+----------------------------------+----------------+----------------------------------------------------+
12+
| option | value | description |
13+
+==================================+================+====================================================+
14+
| eps | ``0.005`` | Viscosity for diffusion |
15+
+----------------------------------+----------------+----------------------------------------------------+
16+
17+
* section: [driver]
18+
19+
+----------------------------------+----------------+----------------------------------------------------+
20+
| option | value | description |
21+
+==================================+================+====================================================+
22+
| cfl | ``0.8`` | advective CFL number |
23+
+----------------------------------+----------------+----------------------------------------------------+
24+
25+
* section: [particles]
26+
27+
+----------------------------------+----------------+----------------------------------------------------+
28+
| option | value | description |
29+
+==================================+================+====================================================+
30+
| do_particles | ``0`` | |
31+
+----------------------------------+----------------+----------------------------------------------------+
32+
| particle_generator | ``grid`` | |
33+
+----------------------------------+----------------+----------------------------------------------------+
34+

pyro/pyro_sim.py

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
"advection_fv4",
1818
"advection_weno",
1919
"burgers",
20+
"viscous_burgers",
2021
"compressible",
2122
"compressible_rk",
2223
"compressible_fv4",

pyro/viscous_burgers/__init__.py

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
"""
2+
The pyro viscous burgers solver. This implements a second-order,
3+
unsplit method for viscous burgers equations.
4+
"""
5+
6+
__all__ = ['simulation']
7+
from .simulation import Simulation

pyro/viscous_burgers/_defaults

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
[driver]
2+
cfl = 0.8 ; advective CFL number
3+
4+
5+
[advection]
6+
limiter = 2 ; limiter (0 = none, 1 = 2nd order, 2 = 4th order)
7+
8+
9+
[particles]
10+
do_particles = 0
11+
particle_generator = grid
12+
13+
[diffusion]
14+
eps = 0.005 ; Viscosity for diffusion

pyro/viscous_burgers/interface.py

Lines changed: 171 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,171 @@
1+
from pyro.multigrid import MG
2+
3+
4+
def get_lap(grid, a):
5+
r"""
6+
Parameters
7+
----------
8+
grid: grid object
9+
a: ndarray
10+
the variable that we want to find the laplacian of
11+
12+
Returns
13+
-------
14+
out : ndarray (laplacian of state a)
15+
"""
16+
17+
lap = grid.scratch_array()
18+
19+
lap.v(buf=2)[:, :] = (a.ip(1, buf=2) - 2.0*a.v(buf=2) + a.ip(-1, buf=2)) / grid.dx**2 \
20+
+ (a.jp(1, buf=2) - 2.0*a.v(buf=2) + a.jp(-1, buf=2)) / grid.dy**2
21+
22+
return lap
23+
24+
25+
def diffuse(my_data, rp, dt, scalar_name, A):
26+
r"""
27+
A routine to solve the Helmhotlz equation with constant coefficient
28+
and update the state.
29+
30+
(a + b \lap) phi = f
31+
32+
using Crank-Nicolson discretization with multigrid V-cycle.
33+
34+
Parameters
35+
----------
36+
my_data : CellCenterData2d object
37+
The data object containing the grid and advective scalar that
38+
we are advecting.
39+
rp : RuntimeParameters object
40+
The runtime parameters for the simulation
41+
dt : float
42+
The timestep we are advancing through.
43+
scalar_name : str
44+
The name of the variable contained in my_data that we are
45+
advecting
46+
A: ndarray
47+
The advective source term for diffusion
48+
49+
Returns
50+
-------
51+
out : ndarray (solution of the Helmholtz equation)
52+
53+
"""
54+
55+
myg = my_data.grid
56+
57+
a = my_data.get_var(scalar_name)
58+
eps = rp.get_param("diffusion.eps")
59+
60+
# Create the multigrid with a constant diffusion coefficient
61+
# the equation has form:
62+
# (alpha - beta L) phi = f
63+
#
64+
# with alpha = 1
65+
# beta = (dt/2) k
66+
# f = phi + (dt/2) k L phi - A
67+
#
68+
# Same as Crank-Nicolson discretization except with an extra
69+
# advection source term)
70+
71+
mg = MG.CellCenterMG2d(myg.nx, myg.ny,
72+
xmin=myg.xmin, xmax=myg.xmax,
73+
ymin=myg.ymin, ymax=myg.ymax,
74+
xl_BC_type=my_data.BCs[scalar_name].xlb,
75+
xr_BC_type=my_data.BCs[scalar_name].xrb,
76+
yl_BC_type=my_data.BCs[scalar_name].ylb,
77+
yr_BC_type=my_data.BCs[scalar_name].yrb,
78+
alpha=1.0, beta=0.5*dt*eps,
79+
verbose=0)
80+
81+
# Compute the RHS: f
82+
83+
f = mg.soln_grid.scratch_array()
84+
85+
lap = get_lap(myg, a)
86+
87+
f.v()[:, :] = a.v() + 0.5*dt*eps * lap.v() - dt*A.v()
88+
89+
mg.init_RHS(f)
90+
91+
# initial guess is zeros
92+
93+
mg.init_zeros()
94+
mg.solve(rtol=1.e-12)
95+
96+
# perform the diffusion update
97+
98+
a.v()[:, :] = mg.get_solution().v()
99+
100+
101+
def apply_diffusion_corrections(grid, dt, eps,
102+
u, v,
103+
u_xl, u_xr,
104+
u_yl, u_yr,
105+
v_xl, v_xr,
106+
v_yl, v_yr):
107+
r"""
108+
Apply diffusion correction term to the interface state
109+
110+
.. math::
111+
112+
u_t + u u_x + v u_y = eps (u_xx + u_yy)
113+
v_t + u v_x + v v_y = eps (v_xx + v_yy)
114+
115+
We use a second-order (piecewise linear) unsplit Godunov method
116+
(following Colella 1990).
117+
118+
Our convection is that the fluxes are going to be defined on the
119+
left edge of the computational zones::
120+
121+
| | | |
122+
| | | |
123+
-+------+------+------+------+------+------+--
124+
| i-1 | i | i+1 |
125+
126+
a_l,i a_r,i a_l,i+1
127+
128+
129+
a_r,i and a_l,i+1 are computed using the information in
130+
zone i,j.
131+
132+
Parameters
133+
----------
134+
grid : Grid2d
135+
The grid object
136+
dt : float
137+
The timestep
138+
eps : float
139+
The viscosity
140+
u_xl, u_xr : ndarray ndarray
141+
left and right states of x-velocity in x-interface.
142+
u_yl, u_yr : ndarray ndarray
143+
left and right states of x-velocity in y-interface.
144+
v_xl, v_xr : ndarray ndarray
145+
left and right states of y-velocity in x-interface.
146+
v_yl, u_yr : ndarray ndarray
147+
left and right states of y-velocity in y-interface.
148+
149+
Returns
150+
-------
151+
out : ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray
152+
unsplit predictions of the left and right states of u and v on both the x- and
153+
y-interfaces along with diffusion correction terms.
154+
"""
155+
156+
#apply diffusion correction to the interface
157+
158+
lap_u = get_lap(grid, u)
159+
lap_v = get_lap(grid, v)
160+
161+
u_xl.ip(1, buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
162+
u_yl.jp(1, buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
163+
u_xr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
164+
u_yr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
165+
166+
v_xl.ip(1, buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
167+
v_yl.jp(1, buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
168+
v_xr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
169+
v_yr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
170+
171+
return u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr

0 commit comments

Comments
 (0)