viscous burgers solver

README.md

@@ -134,6 +134,10 @@ pyro provides the following solvers (all in 2-d):
 
   - `burgers`: a second-order unsplit solver for invsicid Burgers' equation.
 
+  - `viscous_burgers`: a second-order unsplit solver for viscous Burgers' equation
+    with constant-coefficient diffusion. It uses Crank-Nicolson time-discretized
+    solver for solving diffusion.
+
   - `compressible`: a second-order unsplit solver for the Euler
     equations of compressible hydrodynamics.  This uses characteristic
     tracing and corner coupling for the prediction of the interface

docs/source/burgers_basics.rst

@@ -6,7 +6,7 @@ Burgers' Equation is a nonlinear hyperbolic equation. It has the same form as th
 ``Inviscid Burgers``
 --------------------------------
 
-A 2D Burgers' Equation has the following form:
+A 2D inviscid Burgers' Equation has the following form:
 
 .. math::
 
@@ -34,3 +34,30 @@ The parameters for this solver are:
    :align: center
 
 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.
+
+``Viscous Burgers``
+--------------------------------
+
+A 2D viscous Burgers' Equation has the following form:
+
+.. math::
+
+   u_t + u u_x + v u_y = \epsilon \left( u_{xx} + u_{yy}\right) \\
+   v_t + u v_x + v v_y = \epsilon \left( v_{xx} + v_{yy}\right)
+
+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.
+
+: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.
+
+
+The parameters for this solver are:
+
+.. include:: viscous_burgers_defaults.inc
+
+
+.. image:: viscous_burgers.png
+   :align: center
+
+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.
+
+

docs/source/pyro.rst

@@ -18,6 +18,7 @@ Subpackages
    pyro.advection_rk
    pyro.advection_weno
    pyro.burgers
+   pyro.viscous_burgers
    pyro.compressible
    pyro.compressible_fv4
    pyro.compressible_react

docs/source/pyro.viscous_burgers.problems.rst

@@ -0,0 +1,26 @@
+pyro.viscous\_burgers.problems package
+=======================================
+
+.. automodule:: pyro.viscous_burgers.problems
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Submodules
+----------
+
+pyro.viscous\_burgers.problems.converge module
+-----------------------------------------------
+
+.. automodule:: pyro.viscous_burgers.problems.converge
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+pyro.viscous\_burgers.problems.tophat module
+--------------------------------------------
+
+.. automodule:: pyro.viscous_burgers.problems.tophat
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/source/pyro.viscous_burgers.rst

@@ -0,0 +1,34 @@
+pyro.viscous\_burgers package
+==============================
+
+.. automodule:: pyro.viscous_burgers
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Subpackages
+-----------
+
+.. toctree::
+   :maxdepth: 4
+
+   pyro.viscous_burgers.problems
+
+Submodules
+----------
+
+pyro.viscous\_burgers.interface module
+---------------------------------------
+
+.. automodule:: pyro.viscous_burgers.interface
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+pyro.viscous_burgers.simulation module
+---------------------------------------
+
+.. automodule:: pyro.viscous_burgers.simulation
+   :members:
+   :undoc-members:
+   :show-inheritance:

docs/source/viscous_burgers_defaults.inc

@@ -0,0 +1,34 @@
+* section: [advection]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | limiter                          | ``2``          | limiter (0 = none, 1 = 2nd order, 2 = 4th order)   |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [diffusion]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | eps                              | ``0.005``      | Viscosity for diffusion                            |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [driver]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | cfl                              | ``0.8``        | advective CFL number                               |
+  +----------------------------------+----------------+----------------------------------------------------+
+
+* section: [particles]
+
+  +----------------------------------+----------------+----------------------------------------------------+
+  | option                           | value          | description                                        |
+  +==================================+================+====================================================+
+  | do_particles                     | ``0``          |                                                    |
+  +----------------------------------+----------------+----------------------------------------------------+
+  | particle_generator               | ``grid``       |                                                    |
+  +----------------------------------+----------------+----------------------------------------------------+
+

pyro/pyro_sim.py

@@ -17,6 +17,7 @@
                  "advection_fv4",
                  "advection_weno",
                  "burgers",
+                 "viscous_burgers",
                  "compressible",
                  "compressible_rk",
                  "compressible_fv4",

pyro/viscous_burgers/__init__.py

@@ -0,0 +1,7 @@
+"""
+The pyro viscous burgers solver.  This implements a second-order,
+unsplit method for viscous burgers equations.
+"""
+
+__all__ = ['simulation']
+from .simulation import Simulation

pyro/viscous_burgers/_defaults

@@ -0,0 +1,14 @@
+[driver]
+cfl       = 0.8           ; advective CFL number
+
+
+[advection]
+limiter = 2  ; limiter (0 = none, 1 = 2nd order, 2 = 4th order)
+
+
+[particles]
+do_particles = 0
+particle_generator = grid
+
+[diffusion]
+eps = 0.005               ; Viscosity for diffusion

pyro/viscous_burgers/interface.py

@@ -0,0 +1,171 @@
+from pyro.multigrid import MG
+
+
+def get_lap(grid, a):
+    r"""
+    Parameters
+    ----------
+    grid: grid object
+    a: ndarray
+       the variable that we want to find the laplacian of
+
+    Returns
+    -------
+    out : ndarray (laplacian of state a)
+    """
+
+    lap = grid.scratch_array()
+
+    lap.v(buf=2)[:, :] = (a.ip(1, buf=2) - 2.0*a.v(buf=2) + a.ip(-1, buf=2)) / grid.dx**2 \
+                       + (a.jp(1, buf=2) - 2.0*a.v(buf=2) + a.jp(-1, buf=2)) / grid.dy**2
+
+    return lap
+
+
+def diffuse(my_data, rp, dt, scalar_name, A):
+    r"""
+    A routine to solve the Helmhotlz equation with constant coefficient
+    and update the state.
+
+    (a + b \lap) phi = f
+
+    using Crank-Nicolson discretization with multigrid V-cycle.
+
+    Parameters
+    ----------
+    my_data : CellCenterData2d object
+        The data object containing the grid and advective scalar that
+        we are advecting.
+    rp : RuntimeParameters object
+        The runtime parameters for the simulation
+    dt : float
+        The timestep we are advancing through.
+    scalar_name : str
+        The name of the variable contained in my_data that we are
+        advecting
+    A: ndarray
+       The advective source term for diffusion
+
+    Returns
+    -------
+    out : ndarray (solution of the Helmholtz equation)
+
+    """
+
+    myg = my_data.grid
+
+    a = my_data.get_var(scalar_name)
+    eps = rp.get_param("diffusion.eps")
+
+    # Create the multigrid with a constant diffusion coefficient
+    # the equation has form:
+    # (alpha - beta L) phi = f
+    #
+    # with alpha = 1
+    #      beta  = (dt/2) k
+    #      f     = phi + (dt/2) k L phi - A
+    #
+    # Same as Crank-Nicolson discretization except with an extra
+    # advection source term)
+
+    mg = MG.CellCenterMG2d(myg.nx, myg.ny,
+                           xmin=myg.xmin, xmax=myg.xmax,
+                           ymin=myg.ymin, ymax=myg.ymax,
+                           xl_BC_type=my_data.BCs[scalar_name].xlb,
+                           xr_BC_type=my_data.BCs[scalar_name].xrb,
+                           yl_BC_type=my_data.BCs[scalar_name].ylb,
+                           yr_BC_type=my_data.BCs[scalar_name].yrb,
+                           alpha=1.0, beta=0.5*dt*eps,
+                           verbose=0)
+
+    # Compute the RHS: f
+
+    f = mg.soln_grid.scratch_array()
+
+    lap = get_lap(myg, a)
+
+    f.v()[:, :] = a.v() + 0.5*dt*eps * lap.v() - dt*A.v()
+
+    mg.init_RHS(f)
+
+    # initial guess is zeros
+
+    mg.init_zeros()
+    mg.solve(rtol=1.e-12)
+
+    # perform the diffusion update
+
+    a.v()[:, :] = mg.get_solution().v()
+
+
+def apply_diffusion_corrections(grid, dt, eps,
+                                u, v,
+                                u_xl, u_xr,
+                                u_yl, u_yr,
+                                v_xl, v_xr,
+                                v_yl, v_yr):
+    r"""
+    Apply diffusion correction term to the interface state
+
+    .. math::
+
+       u_t  + u u_x  + v u_y  = eps (u_xx + u_yy)
+       v_t  + u v_x  + v v_y  = eps (v_xx + v_yy)
+
+    We use a second-order (piecewise linear) unsplit Godunov method
+    (following Colella 1990).
+
+    Our convection is that the fluxes are going to be defined on the
+    left edge of the computational zones::
+
+        |             |             |             |
+        |             |             |             |
+       -+------+------+------+------+------+------+--
+        |     i-1     |      i      |     i+1     |
+
+                 a_l,i  a_r,i   a_l,i+1
+
+
+    a_r,i and a_l,i+1 are computed using the information in
+    zone i,j.
+
+    Parameters
+    ----------
+    grid : Grid2d
+        The grid object
+    dt : float
+        The timestep
+    eps : float
+         The viscosity
+    u_xl, u_xr : ndarray ndarray
+        left and right states of x-velocity in x-interface.
+    u_yl, u_yr : ndarray ndarray
+        left and right states of x-velocity in y-interface.
+    v_xl, v_xr : ndarray ndarray
+        left and right states of y-velocity in x-interface.
+    v_yl, u_yr : ndarray ndarray
+        left and right states of y-velocity in y-interface.
+
+    Returns
+    -------
+    out : ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray
+        unsplit predictions of the left and right states of u and v on both the x- and
+        y-interfaces along with diffusion correction terms.
+    """
+
+    #apply diffusion correction to the interface
+
+    lap_u = get_lap(grid, u)
+    lap_v = get_lap(grid, v)
+
+    u_xl.ip(1, buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
+    u_yl.jp(1, buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
+    u_xr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
+    u_yr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_u.v(buf=2)
+
+    v_xl.ip(1, buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
+    v_yl.jp(1, buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
+    v_xr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
+    v_yr.v(buf=2)[:, :] += 0.5 * eps * dt * lap_v.v(buf=2)
+
+    return u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr

pyro/viscous_burgers/problems/__init__.py

@@ -0,0 +1 @@
+../../burgers/problems/__init__.py

pyro/viscous_burgers/problems/_converge.defaults

@@ -0,0 +1 @@
+../../burgers/problems/_converge.defaults

pyro/viscous_burgers/problems/_test.defaults

@@ -0,0 +1 @@
+../../burgers/problems/_test.defaults

pyro/viscous_burgers/problems/_tophat.defaults

@@ -0,0 +1 @@
+../../burgers/problems/_tophat.defaults

pyro/viscous_burgers/problems/converge.py

@@ -0,0 +1 @@
+../../burgers/problems/converge.py