From 279776bda92577615437de5a475cb188f4fb6c0f Mon Sep 17 00:00:00 2001 From: Zhi Chen <62574124+zhichen3@users.noreply.github.com> Date: Wed, 9 Aug 2023 19:57:23 -0400 Subject: [PATCH] viscous burgers solver (#171) this adds a viscous Burgers' solver by adding parabolic solves for the viscous terms --- README.md | 4 + docs/source/burgers_basics.rst | 29 ++- docs/source/pyro.rst | 1 + docs/source/pyro.viscous_burgers.problems.rst | 26 +++ docs/source/pyro.viscous_burgers.rst | 34 ++++ docs/source/viscous_burgers_defaults.inc | 34 ++++ pyro/pyro_sim.py | 1 + pyro/viscous_burgers/__init__.py | 7 + pyro/viscous_burgers/_defaults | 14 ++ pyro/viscous_burgers/interface.py | 171 ++++++++++++++++++ pyro/viscous_burgers/problems/__init__.py | 1 + .../problems/_converge.defaults | 1 + pyro/viscous_burgers/problems/_test.defaults | 1 + .../viscous_burgers/problems/_tophat.defaults | 1 + pyro/viscous_burgers/problems/converge.py | 1 + .../problems/inputs.converge.128 | 40 ++++ .../problems/inputs.converge.256 | 40 ++++ .../problems/inputs.converge.32 | 39 ++++ .../problems/inputs.converge.64 | 40 ++++ pyro/viscous_burgers/problems/inputs.smooth | 36 ++++ pyro/viscous_burgers/problems/inputs.test | 38 ++++ pyro/viscous_burgers/problems/inputs.tophat | 34 ++++ pyro/viscous_burgers/problems/test.py | 1 + pyro/viscous_burgers/problems/tophat.py | 1 + pyro/viscous_burgers/simulation.py | 89 +++++++++ 25 files changed, 683 insertions(+), 1 deletion(-) create mode 100644 docs/source/pyro.viscous_burgers.problems.rst create mode 100644 docs/source/pyro.viscous_burgers.rst create mode 100644 docs/source/viscous_burgers_defaults.inc create mode 100644 pyro/viscous_burgers/__init__.py create mode 100644 pyro/viscous_burgers/_defaults create mode 100644 pyro/viscous_burgers/interface.py create mode 120000 pyro/viscous_burgers/problems/__init__.py create mode 120000 pyro/viscous_burgers/problems/_converge.defaults create mode 120000 pyro/viscous_burgers/problems/_test.defaults create mode 120000 pyro/viscous_burgers/problems/_tophat.defaults create mode 120000 pyro/viscous_burgers/problems/converge.py create mode 100644 pyro/viscous_burgers/problems/inputs.converge.128 create mode 100644 pyro/viscous_burgers/problems/inputs.converge.256 create mode 100644 pyro/viscous_burgers/problems/inputs.converge.32 create mode 100644 pyro/viscous_burgers/problems/inputs.converge.64 create mode 100644 pyro/viscous_burgers/problems/inputs.smooth create mode 100644 pyro/viscous_burgers/problems/inputs.test create mode 100644 pyro/viscous_burgers/problems/inputs.tophat create mode 120000 pyro/viscous_burgers/problems/test.py create mode 120000 pyro/viscous_burgers/problems/tophat.py create mode 100644 pyro/viscous_burgers/simulation.py diff --git a/README.md b/README.md index 17813d7e7..4fe0b512b 100644 --- a/README.md +++ b/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 diff --git a/docs/source/burgers_basics.rst b/docs/source/burgers_basics.rst index 8cc5af11e..51162edbe 100644 --- a/docs/source/burgers_basics.rst +++ b/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. + + diff --git a/docs/source/pyro.rst b/docs/source/pyro.rst index da63fd3d8..1d8831dcb 100644 --- a/docs/source/pyro.rst +++ b/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 diff --git a/docs/source/pyro.viscous_burgers.problems.rst b/docs/source/pyro.viscous_burgers.problems.rst new file mode 100644 index 000000000..12c3ecb41 --- /dev/null +++ b/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: diff --git a/docs/source/pyro.viscous_burgers.rst b/docs/source/pyro.viscous_burgers.rst new file mode 100644 index 000000000..4df56a120 --- /dev/null +++ b/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: diff --git a/docs/source/viscous_burgers_defaults.inc b/docs/source/viscous_burgers_defaults.inc new file mode 100644 index 000000000..6308266a3 --- /dev/null +++ b/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`` | | + +----------------------------------+----------------+----------------------------------------------------+ + diff --git a/pyro/pyro_sim.py b/pyro/pyro_sim.py index 4d5ff8d2b..28973c39c 100755 --- a/pyro/pyro_sim.py +++ b/pyro/pyro_sim.py @@ -17,6 +17,7 @@ "advection_fv4", "advection_weno", "burgers", + "viscous_burgers", "compressible", "compressible_rk", "compressible_fv4", diff --git a/pyro/viscous_burgers/__init__.py b/pyro/viscous_burgers/__init__.py new file mode 100644 index 000000000..cb2105614 --- /dev/null +++ b/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 diff --git a/pyro/viscous_burgers/_defaults b/pyro/viscous_burgers/_defaults new file mode 100644 index 000000000..3b5842e5c --- /dev/null +++ b/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 diff --git a/pyro/viscous_burgers/interface.py b/pyro/viscous_burgers/interface.py new file mode 100644 index 000000000..131193e6e --- /dev/null +++ b/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 diff --git a/pyro/viscous_burgers/problems/__init__.py b/pyro/viscous_burgers/problems/__init__.py new file mode 120000 index 000000000..9245b6099 --- /dev/null +++ b/pyro/viscous_burgers/problems/__init__.py @@ -0,0 +1 @@ +../../burgers/problems/__init__.py \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/_converge.defaults b/pyro/viscous_burgers/problems/_converge.defaults new file mode 120000 index 000000000..c44b31278 --- /dev/null +++ b/pyro/viscous_burgers/problems/_converge.defaults @@ -0,0 +1 @@ +../../burgers/problems/_converge.defaults \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/_test.defaults b/pyro/viscous_burgers/problems/_test.defaults new file mode 120000 index 000000000..963fcf562 --- /dev/null +++ b/pyro/viscous_burgers/problems/_test.defaults @@ -0,0 +1 @@ +../../burgers/problems/_test.defaults \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/_tophat.defaults b/pyro/viscous_burgers/problems/_tophat.defaults new file mode 120000 index 000000000..0b932730d --- /dev/null +++ b/pyro/viscous_burgers/problems/_tophat.defaults @@ -0,0 +1 @@ +../../burgers/problems/_tophat.defaults \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/converge.py b/pyro/viscous_burgers/problems/converge.py new file mode 120000 index 000000000..53b975cbf --- /dev/null +++ b/pyro/viscous_burgers/problems/converge.py @@ -0,0 +1 @@ +../../burgers/problems/converge.py \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/inputs.converge.128 b/pyro/viscous_burgers/problems/inputs.converge.128 new file mode 100644 index 000000000..7716e4851 --- /dev/null +++ b/pyro/viscous_burgers/problems/inputs.converge.128 @@ -0,0 +1,40 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 500 +tmax = 1.0 + +max_dt_change = 1.e33 +init_tstep_factor = 1.0 +fix_dt = 0.0025 + +[driver] +cfl = 0.8 + +[io] +basename = converge.128_ +dt_out = 0.2 + +[mesh] +nx = 128 +ny = 128 +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = periodic +yrboundary = periodic + + +[advection] +limiter = 0 + + +[particles] +do_particles = 1 +particle_generator = grid + +[diffusion] +eps = 0.05 \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/inputs.converge.256 b/pyro/viscous_burgers/problems/inputs.converge.256 new file mode 100644 index 000000000..2f067fb6c --- /dev/null +++ b/pyro/viscous_burgers/problems/inputs.converge.256 @@ -0,0 +1,40 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 1000 +tmax = 1.0 + +max_dt_change = 1.e33 +init_tstep_factor = 1.0 +fix_dt = 0.00125 + +[driver] +cfl = 0.8 + +[io] +basename = converge.256_ +dt_out = 0.2 + +[mesh] +nx = 256 +ny = 256 +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = periodic +yrboundary = periodic + + +[advection] +limiter = 0 + + +[particles] +do_particles = 1 +particle_generator = grid + +[diffusion] +eps = 0.05 \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/inputs.converge.32 b/pyro/viscous_burgers/problems/inputs.converge.32 new file mode 100644 index 000000000..0b1bb24f1 --- /dev/null +++ b/pyro/viscous_burgers/problems/inputs.converge.32 @@ -0,0 +1,39 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 500 +tmax = 1.0 + +max_dt_change = 1.e33 +init_tstep_factor = 1.0 +fix_dt = 0.01 + +[driver] +cfl = 0.8 + +[io] +basename = converge.32_ +dt_out = 0.2 + +[mesh] +nx = 32 +ny = 32 +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = periodic +yrboundary = periodic + + +[advection] +limiter = 0 + +[particles] +do_particles = 1 +particle_generator = grid + +[diffusion] +eps = 0.05 diff --git a/pyro/viscous_burgers/problems/inputs.converge.64 b/pyro/viscous_burgers/problems/inputs.converge.64 new file mode 100644 index 000000000..bf6513932 --- /dev/null +++ b/pyro/viscous_burgers/problems/inputs.converge.64 @@ -0,0 +1,40 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 500 +tmax = 1.0 + +max_dt_change = 1.e33 +init_tstep_factor = 1.0 +fix_dt = 0.005 + +[driver] +cfl = 0.8 + +[io] +basename = converge.64_ +dt_out = 0.2 + +[mesh] +nx = 64 +ny = 64 +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = periodic +yrboundary = periodic + + +[advection] +limiter = 0 + + +[particles] +do_particles = 1 +particle_generator = grid + +[diffusion] +eps = 0.05 \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/inputs.smooth b/pyro/viscous_burgers/problems/inputs.smooth new file mode 100644 index 000000000..9b1f5ff76 --- /dev/null +++ b/pyro/viscous_burgers/problems/inputs.smooth @@ -0,0 +1,36 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 500 +tmax = 1.0 + +max_dt_change = 1.e33 +init_tstep_factor = 1.0 + + +[driver] +cfl = 0.8 + +[io] +basename = smooth_ +dt_out = 0.2 + +[mesh] +nx = 32 +ny = 32 +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = periodic +yrboundary = periodic + + +[advection] +limiter = 2 + +[particles] +do_particles = 1 +particle_generator = grid diff --git a/pyro/viscous_burgers/problems/inputs.test b/pyro/viscous_burgers/problems/inputs.test new file mode 100644 index 000000000..6b0afb863 --- /dev/null +++ b/pyro/viscous_burgers/problems/inputs.test @@ -0,0 +1,38 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 500 +tmax = 0.1 + +max_dt_change = 1.e33 +init_tstep_factor = 1.0 + +cfl = 0.8 + + +[io] +basename = test_ +n_out = 10 + +[mesh] +nx = 128 +ny = 128 +xmax = 1.0 +ymax = 1.0 + +xlboundary = outflow +xrboundary = outflow + +ylboundary = outflow +yrboundary = outflow + + +[advection] +limiter = 2 + +[particles] +do_particles = 1 +particle_generator = grid + +[diffusion] +eps = 0.05 \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/inputs.tophat b/pyro/viscous_burgers/problems/inputs.tophat new file mode 100644 index 000000000..1766bff48 --- /dev/null +++ b/pyro/viscous_burgers/problems/inputs.tophat @@ -0,0 +1,34 @@ +# simple inputs files for the unsplit CTU hydro scheme + +[driver] +max_steps = 500 +tmax = 1.0 + +max_dt_change = 1.e33 +init_tstep_factor = 1.0 + +cfl = 0.8 + + +[io] +basename = tophat_ +n_out = 10 + +[mesh] +nx = 32 +ny = 32 +xmax = 1.0 +ymax = 1.0 + +xlboundary = periodic +xrboundary = periodic + +ylboundary = periodic +yrboundary = periodic + + +[advection] +limiter = 2 + +[diffusion] +eps = 0.05 diff --git a/pyro/viscous_burgers/problems/test.py b/pyro/viscous_burgers/problems/test.py new file mode 120000 index 000000000..74e979bdd --- /dev/null +++ b/pyro/viscous_burgers/problems/test.py @@ -0,0 +1 @@ +../../burgers/problems/test.py \ No newline at end of file diff --git a/pyro/viscous_burgers/problems/tophat.py b/pyro/viscous_burgers/problems/tophat.py new file mode 120000 index 000000000..99fb40d04 --- /dev/null +++ b/pyro/viscous_burgers/problems/tophat.py @@ -0,0 +1 @@ +../../burgers/problems/tophat.py \ No newline at end of file diff --git a/pyro/viscous_burgers/simulation.py b/pyro/viscous_burgers/simulation.py new file mode 100644 index 000000000..bfbc7bfc4 --- /dev/null +++ b/pyro/viscous_burgers/simulation.py @@ -0,0 +1,89 @@ +from pyro.burgers import Simulation as burgers_sim +from pyro.burgers import burgers_interface +from pyro.mesh import reconstruction +from pyro.viscous_burgers import interface + + +class Simulation(burgers_sim): + + def evolve(self): + """ + Evolve the viscous burgers equation through one timestep. + """ + + myg = self.cc_data.grid + + u = self.cc_data.get_var("x-velocity") + v = self.cc_data.get_var("y-velocity") + + # -------------------------------------------------------------------------- + # monotonized central differences + # -------------------------------------------------------------------------- + + limiter = self.rp.get_param("advection.limiter") + eps = self.rp.get_param("diffusion.eps") + + # Give da/dx and da/dy using input: (state, grid, direction, limiter) + + ldelta_ux = reconstruction.limit(u, myg, 1, limiter) + ldelta_uy = reconstruction.limit(u, myg, 2, limiter) + ldelta_vx = reconstruction.limit(v, myg, 1, limiter) + ldelta_vy = reconstruction.limit(v, myg, 2, limiter) + + # Compute the advective fluxes + # Get u, v fluxes + + # Get the interface states without transverse or diffusion corrections + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = burgers_interface.get_interface_states(myg, self.dt, + u, v, + ldelta_ux, ldelta_vx, + ldelta_uy, ldelta_vy) + + # Apply diffusion correction terms to the interface states + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = interface.apply_diffusion_corrections(myg, self.dt, eps, + u, v, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr) + + # Apply transverse correction terms to the interface states + u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = burgers_interface.apply_transverse_corrections(myg, self.dt, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr) + + # Construct the interface fluxes + u_flux_x, u_flux_y, v_flux_x, v_flux_y = burgers_interface.construct_unsplit_fluxes(myg, + u_xl, u_xr, + u_yl, u_yr, + v_xl, v_xr, + v_yl, v_yr) + + # Compute the advective source terms for diffusion using flux computed above + + A_u = myg.scratch_array() + A_v = myg.scratch_array() + + A_u.v()[:, :] = (u_flux_x.ip(1) - u_flux_x.v()) / myg.dx + \ + (u_flux_y.jp(1) - u_flux_y.v()) / myg.dy + + A_v.v()[:, :] = (v_flux_x.ip(1) - v_flux_x.v()) / myg.dx + \ + (v_flux_y.jp(1) - v_flux_y.v()) / myg.dy + + # Update state by doing diffusion update with extra advective source term. + + interface.diffuse(self.cc_data, self.rp, self.dt, "x-velocity", A_u) + interface.diffuse(self.cc_data, self.rp, self.dt, "y-velocity", A_v) + + if self.particles is not None: + + u2d = self.cc_data.get_var("x-velocity") + v2d = self.cc_data.get_var("y-velocity") + + self.particles.update_particles(self.dt, u2d, v2d) + + # increment the time + self.cc_data.t += self.dt + self.n += 1