diff --git a/.gitignore b/.gitignore index 3c78a9482..861dd3981 100755 --- a/.gitignore +++ b/.gitignore @@ -9,6 +9,8 @@ dist/ inputs.auto .noseids +*problems.inc + .ipynb_checkpoints/ *.pdf diff --git a/docs/Makefile b/docs/Makefile index a10024564..db840c4cc 100644 --- a/docs/Makefile +++ b/docs/Makefile @@ -17,4 +17,5 @@ help: # Catch-all target: route all unknown targets to Sphinx using the new # "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). %: Makefile + ./document_problems.py @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/document_problems.py b/docs/document_problems.py new file mode 100755 index 000000000..0585fab58 --- /dev/null +++ b/docs/document_problems.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python + +import importlib +from pathlib import Path + +SOLVERS = ["advection", + "advection_fv4", + "advection_nonuniform", + "advection_rk", + "advection_weno", + "burgers", + "burgers_viscous", + "compressible", + "compressible_fv4", + "compressible_react", + "compressible_rk", + "compressible_sdc", + "diffusion", + "incompressible", + "incompressible_viscous", + "lm_atm", + "swe"] + +MAX_LEN = 36 + + +def doit(pyro_home): + + for s in SOLVERS: + + # check it the problems/ directory is not a softlink to + # a different solver + + p = Path(f"{pyro_home}/pyro/{s}/problems").resolve() + + with open(f"{pyro_home}/docs/source/{s}_problems.inc", "w") as finc: + + finc.write("supported problems\n") + finc.write("------------------\n") + + if (parent_solver := p.parts[-2]) == s: + + # find all the problems + for prob in p.glob("*.py"): + if prob.name == "__init__.py": + continue + + mprob = importlib.import_module(f"pyro.{s}.problems.{prob.stem}") + + if "init_data" not in dir(mprob): + # not a problem setup + continue + + finc.write(f"``{prob.stem}``\n") + finc.write("^" * (len(prob.stem)+4) + "\n\n") + + if mprob.__doc__: + finc.write(mprob.__doc__) + + finc.write("\n") + + try: + params = mprob.PROBLEM_PARAMS + except AttributeError: + params = {} + + if params: + finc.write("parameters: \n\n") + + finc.write("="*MAX_LEN + " " + "="*MAX_LEN + "\n") + finc.write(f"{'name':{MAX_LEN}} {'default':{MAX_LEN}}\n") + finc.write("="*MAX_LEN + " " + "="*MAX_LEN + "\n") + + for k, v in params.items(): + pname = "``" + k + "``" + finc.write(f"{pname:{MAX_LEN}} {v:{MAX_LEN}}\n") + + finc.write("="*MAX_LEN + " " + "="*MAX_LEN + "\n") + + + finc.write("\n\n") + + else: + finc.write(f"``{s}`` uses the problems defined by ``{parent_solver}``.") + + +if __name__ == "__main__": + doit("..") diff --git a/docs/source/advection_basics.rst b/docs/source/advection_basics.rst index 2f467d038..c51cf8067 100644 --- a/docs/source/advection_basics.rst +++ b/docs/source/advection_basics.rst @@ -1,5 +1,6 @@ -Advection solvers -================= +********* +Advection +********* The linear advection equation: @@ -14,7 +15,7 @@ pyro has several solvers for linear advection, which solve the equation with different spatial and temporal integration schemes. ``advection`` solver --------------------- +==================== :py:mod:`pyro.advection` implements the directionally unsplit corner transport upwind algorithm :cite:`colella:1990` with piecewise linear reconstruction. @@ -29,9 +30,10 @@ The parameters for this solver are: .. include:: advection_defaults.inc +.. include:: advection_problems.inc ``advection_fv4`` solver ------------------------- +======================== :py:mod:`pyro.advection_fv4` uses a fourth-order accurate finite-volume method with RK4 time integration, following the ideas in @@ -50,8 +52,10 @@ The parameters for this solver are: .. include:: advection_fv4_defaults.inc +.. include:: advection_fv4_problems.inc + ``advection_nonuniform`` solver -------------------------------- +=============================== :py:mod:`pyro.advection_nonuniform` models advection with a non-uniform velocity field. This is used to implement the slotted disk problem @@ -62,8 +66,11 @@ The parameters for this solver are: .. include:: advection_nonuniform_defaults.inc +.. include:: advection_nonuniform_problems.inc + + ``advection_rk`` solver ------------------------ +======================= :py:mod:`pyro.advection_rk` uses a method of lines time-integration approach with piecewise linear spatial reconstruction for linear @@ -76,8 +83,10 @@ The parameter for this solver are: .. include:: advection_rk_defaults.inc +.. include:: advection_rk_problems.inc + ``advection_weno`` solver -------------------------- +========================= :py:mod:`pyro.advection_weno` uses a WENO reconstruction and method of lines time-integration @@ -87,9 +96,11 @@ The main parameters that affect this solver are: .. include:: advection_weno_defaults.inc +.. include:: advection_weno_problems.inc + General ideas -------------- +============= The main use for the advection solver is to understand how Godunov techniques work for hyperbolic problems. These same ideas will be used @@ -108,10 +119,10 @@ reconstruction, evolution, and averaging steps: Examples --------- +======== smooth -^^^^^^ +------ The smooth problem initializes a Gaussian profile and advects it with :math:`u = v = 1` through periodic boundaries for a period. The result is that @@ -147,22 +158,15 @@ with the ``advection_fv4`` solver. Departures from perfect scaling are likely due to the use of limiters. -tophat -^^^^^^ - -The tophat problem initializes a circle in the center of the domain -with value 1, and 0 outside. This has very steep jumps, and the -limiters will kick in strongly here. - Exercises ---------- +========= The best way to learn these methods is to play with them yourself. The exercises below are suggestions for explorations and features to add to the advection solver. Explorations -^^^^^^^^^^^^ +------------ * Test the convergence of the solver for a variety of initial conditions (tophat hat will differ from the smooth case because of @@ -175,7 +179,7 @@ Explorations problem?) Extensions -^^^^^^^^^^ +---------- * Implement a dimensionally split version of the advection algorithm. How does the solution compare between the unsplit and diff --git a/docs/source/burgers_basics.rst b/docs/source/burgers_basics.rst index 9f299d1ea..523ad6a9b 100644 --- a/docs/source/burgers_basics.rst +++ b/docs/source/burgers_basics.rst @@ -1,10 +1,13 @@ +***************** Burgers' Equation -================== +***************** -Burgers' Equation is a nonlinear hyperbolic equation. It has the same form as the advection equation, except that the quantity being advected is the velocity itself. +Burgers' Equation is a nonlinear hyperbolic equation. It has the same +form as the advection equation, except that the quantity being +advected is the velocity itself. ``Inviscid Burgers`` --------------------------------- +==================== A 2D inviscid Burgers' Equation has the following form: @@ -29,14 +32,25 @@ The parameters for this solver are: .. include:: burgers_defaults.inc +.. include:: burgers_problems.inc + +Example +------- .. image:: burgers.png :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. +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: @@ -60,6 +74,10 @@ The parameters for this solver are: .. include:: burgers_viscous_defaults.inc +.. include:: burgers_problems.inc + +Example +------- .. image:: viscous_burgers.png :align: center diff --git a/docs/source/compressible_basics.rst b/docs/source/compressible_basics.rst index deb7c1f44..9ceb8259b 100644 --- a/docs/source/compressible_basics.rst +++ b/docs/source/compressible_basics.rst @@ -1,5 +1,6 @@ -Compressible hydrodynamics solvers -================================== +************************** +Compressible hydrodynamics +************************** The Euler equations of compressible hydrodynamics take the form: @@ -27,7 +28,7 @@ direction is allowed. directory in the solver's directory. ``compressible`` solver ------------------------ +======================= :py:mod:`pyro.compressible` is based on a directionally unsplit (the corner transport upwind algorithm) piecewise linear method for the Euler @@ -38,8 +39,11 @@ The parameters for this solver are: .. include:: compressible_defaults.inc +.. include:: compressible_problems.inc + + ``compressible_rk`` solver --------------------------- +========================== :py:mod:`pyro.compressible_rk` uses a method of lines time-integration approach with piecewise linear spatial reconstruction for the Euler @@ -49,8 +53,10 @@ The parameters for this solver are: .. include:: compressible_rk_defaults.inc +.. include:: compressible_rk_problems.inc + ``compressible_fv4`` solver ---------------------------- +=========================== :py:mod:`pyro.compressible_fv4` uses a 4th order accurate method with RK4 time integration, following :cite:`mccorquodalecolella`. @@ -59,9 +65,10 @@ The parameter for this solver are: .. include:: compressible_fv4_defaults.inc +.. include:: compressible_fv4_problems.inc ``compressible_sdc`` solver ---------------------------- +=========================== :py:mod:`pyro.compressible_sdc` uses a 4th order accurate method with spectral-deferred correction (SDC) for the time integration. This @@ -72,9 +79,10 @@ The parameters for this solver are: .. include:: compressible_sdc_defaults.inc +.. include:: compressible_sdc_problems.inc Example problems ----------------- +================ .. note:: @@ -88,7 +96,7 @@ Example problems Sod -^^^ +--- The Sod problem is a standard hydrodynamics problem. It is a one-dimensional shock tube (two states separated by an interface), @@ -119,7 +127,7 @@ is discussed in the notes above, and can be improved in the PPM method with contact steepening. Sedov -^^^^^ +----- The Sedov blast wave problem is another standard test with an analytic solution (Sedov 1959). A lot of energy is point into a point in a @@ -154,7 +162,7 @@ This shows good agreement with the analytic solution. quad -^^^^ +---- The quad problem sets up different states in four regions of the domain and watches the complex interfaces that develop as shocks @@ -172,7 +180,7 @@ online by Pawel Artymowicz). It is run as: rt -^^ +-- The Rayleigh-Taylor problem puts a dense fluid over a lighter one and perturbs the interface with a sinusoidal velocity. Hydrostatic @@ -192,7 +200,7 @@ escape the domain. It is run as: bubble -^^^^^^ +------ The bubble problem initializes a hot spot in a stratified domain and watches it buoyantly rise and roll up. This is run as: @@ -211,10 +219,10 @@ above that rains down on our atmosphere. Also note the acoustic signal propagating outward from the bubble (visible in the U and e panels). Exercises ---------- +========= Explorations -^^^^^^^^^^^^ +------------ * Measure the growth rate of the Rayleigh-Taylor instability for different wavenumbers. @@ -228,7 +236,7 @@ Explorations Extensions -^^^^^^^^^^ +---------- * Limit on the characteristic variables instead of the primitive variables. What changes do you see? (the notes show how to implement @@ -246,16 +254,3 @@ Extensions * Swap the piecewise linear reconstruction for piecewise parabolic (PPM). The notes and the Miller and Colella paper provide a good basis for this. Research the Roe Riemann solver and implement it in pyro. - - -Going further -------------- - -The compressible algorithm presented here is essentially the -single-grid hydrodynamics algorithm used in the `Castro code `_—an -adaptive mesh radiation hydrodynamics code developed at -CCSE/LBNL. `Castro is freely available for download `_. - -A simple, pure Fortran, 1-d compressible hydrodynamics code that does -piecewise constant, linear, or parabolic (PPM) reconstruction is also -available. See the `hydro1d `_ page. diff --git a/docs/source/diffusion_basics.rst b/docs/source/diffusion_basics.rst index 77b71d3b6..050800a8f 100644 --- a/docs/source/diffusion_basics.rst +++ b/docs/source/diffusion_basics.rst @@ -1,5 +1,9 @@ +********* Diffusion -========= +********* + +``diffusion`` solver +==================== pyro solves the constant-conductivity diffusion equation: @@ -16,11 +20,14 @@ multigrid class. The main parameters that affect this solver are: .. include:: diffusion_defaults.inc +.. include:: diffusion_problems.inc + + Examples --------- +======== gaussian -^^^^^^^^ +-------- The gaussian problem initializes a strongly peaked Gaussian centered in the domain. The analytic solution for this shows that the profile @@ -46,14 +53,14 @@ restricted in range to bring out the detail at later times. Exercises ---------- +========= The best way to learn these methods is to play with them yourself. The exercises below are suggestions for explorations and features to add to the advection solver. Explorations -^^^^^^^^^^^^ +------------ * Test the convergence of the solver by varying the resolution and comparing to the analytic solution. @@ -64,7 +71,7 @@ Explorations * Setup some other profiles and experiment with different boundary conditions. Extensions -^^^^^^^^^^ +---------- * Switch from Crank-Nicolson (2nd order in time) to backward Euler (1st order in time) and compare the solution and convergence. This diff --git a/docs/source/incompressible_basics.rst b/docs/source/incompressible_basics.rst index c1ed33c00..053044263 100644 --- a/docs/source/incompressible_basics.rst +++ b/docs/source/incompressible_basics.rst @@ -1,5 +1,13 @@ -Incompressible hydrodynamics solver -=================================== +**************************** +Incompressible hydrodynamics +**************************** + +pyro has two different incompressible solvers: ``incompressible`` is +inviscid and ``incompressible_viscous`` has viscosity. + + +``incompressible`` solver +========================= pyro's incompressible solver solves: @@ -20,6 +28,8 @@ The main parameters that affect this solver are: .. include:: incompressible_defaults.inc +.. include:: incompressible_problems.inc + Examples -------- @@ -82,11 +92,119 @@ The dashed line is second order convergence. We see almost second order behavior with the limiters enabled and slightly better than second order with no limiting. +``incompressible_viscous`` solver +================================= + +pyro's incompressible viscous solver solves: + +.. math:: + + \frac{\partial U}{\partial t} + U \cdot \nabla U + \nabla p &= \nu \nabla^2 U \\ + \nabla \cdot U &= 0 + +This is based on the ``incompressible`` solver, but modifies the +velocity update step to take viscosity into account, by solving two +parabolic equations (one for each velocity component) using multigrid. + +The main parameters that affect this solver are: + +.. include:: incompressible_viscous_defaults.inc + +.. include:: incompressible_viscous_problems.inc + +Examples +-------- + +shear +^^^^^ + +The same shear problem as in incompressible solver, here with viscosity +added. + +.. prompt:: bash + + pyro_sim.py incompressible_viscous shear inputs.shear + +.. image:: shear_viscous.png + :align: center + +Compare this with the inviscid result. Notice how the velocities have +diffused in all directions. + +cavity +^^^^^^ + +The lid-driven cavity is a well-known benchmark problem for hydro codes +(see e.g. :cite:t:`ghia1982`, :cite:t:`Kuhlmann2019`). In a unit square box +with initially static fluid, motion is initiated by a "lid" at the top +boundary, moving to the right with unit velocity. The basic command is: + +.. prompt:: bash + + pyro_sim.py incompressible_viscous cavity inputs.cavity + +It is interesting to observe what happens when varying the viscosity, or, +equivalently the Reynolds number (in this case :math:`\rm{Re}=1/\nu` since +the characteristic length and velocity scales are 1 by default). + +|pic1| |pic2| |pic3| + +.. |pic1| image:: cavity_Re100.png + :width: 32% + +.. |pic2| image:: cavity_Re400.png + :width: 32% + +.. |pic3| image:: cavity_Re1000.png + :width: 32% + +These plots were made by allowing the code to run for longer and approach a +steady-state with the option ``driver.max_steps=1000``, then running +(e.g. for the Re=100 case): + +.. prompt:: bash + + python incompressible_viscous/problems/plot_cavity.py cavity_n64_Re100_0406.h5 -Re 100 -o cavity_Re100.png + +convergence +^^^^^^^^^^^ + +This is the same test as in the incompressible solver. With viscosity, +an exponential term is added to the solution. Limiting can again be +disabled by adding ``incompressible.limiter=0`` to the run command. +The basic set of tests shown below are run as: + +.. prompt:: bash + + pyro_sim.py incompressible_viscous converge inputs.converge.32 vis.dovis=0 + pyro_sim.py incompressible_viscous converge inputs.converge.64 vis.dovis=0 + pyro_sim.py incompressible_viscous converge inputs.converge.128 vis.dovis=0 + +The error is measured by comparing with the analytic solution using +the routine ``incomp_viscous_converge_error.py`` in ``analysis/``. To +generate the plot below, run + +.. prompt:: bash + + python incompressible_viscous/tests/convergence_errors.py convergence_errors.txt + +or ``convergence_errors_no_limiter.txt`` after running with that option. Then: + +.. prompt:: bash + + python incompressible_viscous/tests/convergence_plot.py + +.. image:: incomp_viscous_converge.png + :align: center + +The solver is converging but below second-order, unlike the inviscid case. Limiting +does not seem to make a difference here. + Exercises ---------- +========= Explorations -^^^^^^^^^^^^ +------------ * Disable the MAC projection and run the converge problem—is the method still 2nd order? @@ -94,29 +212,20 @@ Explorations * Experiment with what is projected. Try projecting :math:`U_t` to see if that makes a difference. +* In the lid-driven cavity problem, when does the solution reach a steady-state? + +* :cite:`ghia1982` give benchmark velocities at different Reynolds number for the + lid-driven cavity problem (see their Table I). Do we agree with their results? + Extensions -^^^^^^^^^^ +---------- * Switch the final projection from a cell-centered approximate projection to a nodal projection. This will require writing a new multigrid solver that operates on nodal data. -* Add viscosity to the system. This will require doing 2 parabolic - solves (one for each velocity component). These solves will look - like the diffusion operation, and will update the provisional - velocity field. - * Switch to a variable density system. This will require adding a mass continuity equation that is advected and switching the projections to a variable-coefficient form (since ρ now enters). -Going further -------------- - -The incompressible algorithm presented here is a simplified version of -the projection methods used in the `Maestro low Mach number -hydrodynamics code `_. Maestro -can do variable-density incompressible, anelastic, and low Mach number -stratified flows in stellar (and terrestrial) environments in close -hydrostatic equilibrium. diff --git a/docs/source/incompressible_viscous_basics.rst b/docs/source/incompressible_viscous_basics.rst deleted file mode 100644 index d5476dcc1..000000000 --- a/docs/source/incompressible_viscous_basics.rst +++ /dev/null @@ -1,115 +0,0 @@ -Incompressible viscous hydrodynamics solver -=========================================== - -pyro's incompressible viscous solver solves: - -.. math:: - - \frac{\partial U}{\partial t} + U \cdot \nabla U + \nabla p &= \nu \nabla^2 U \\ - \nabla \cdot U &= 0 - -This solver is based on pyro's incompressible solver, but modifies the -velocity update step to take viscosity into account. - -The main parameters that affect this solver are: - -.. include:: incompressible_viscous_defaults.inc - -Examples --------- - -shear -^^^^^ - -The same shear problem as in incompressible solver, here with viscosity -added. - -.. prompt:: bash - - pyro_sim.py incompressible_viscous shear inputs.shear - -.. image:: shear_viscous.png - :align: center - -Compare this with the inviscid result. Notice how the velocities have -diffused in all directions. - -cavity -^^^^^^ - -The lid-driven cavity is a well-known benchmark problem for hydro codes -(see e.g. :cite:t:`ghia1982`, :cite:t:`Kuhlmann2019`). In a unit square box -with initially static fluid, motion is initiated by a "lid" at the top -boundary, moving to the right with unit velocity. The basic command is: - -.. prompt:: bash - - pyro_sim.py incompressible_viscous cavity inputs.cavity - -It is interesting to observe what happens when varying the viscosity, or, -equivalently the Reynolds number (in this case :math:`\rm{Re}=1/\nu` since -the characteristic length and velocity scales are 1 by default). - -|pic1| |pic2| |pic3| - -.. |pic1| image:: cavity_Re100.png - :width: 32% - -.. |pic2| image:: cavity_Re400.png - :width: 32% - -.. |pic3| image:: cavity_Re1000.png - :width: 32% - -These plots were made by allowing the code to run for longer and approach a -steady-state with the option ``driver.max_steps=1000``, then running -(e.g. for the Re=100 case): - -.. prompt:: bash - - python incompressible_viscous/problems/plot_cavity.py cavity_n64_Re100_0406.h5 -Re 100 -o cavity_Re100.png - -convergence -^^^^^^^^^^^ - -This is the same test as in the incompressible solver. With viscosity, -an exponential term is added to the solution. Limiting can again be -disabled by adding ``incompressible.limiter=0`` to the run command. -The basic set of tests shown below are run as: - -.. prompt:: bash - - pyro_sim.py incompressible_viscous converge inputs.converge.32 vis.dovis=0 - pyro_sim.py incompressible_viscous converge inputs.converge.64 vis.dovis=0 - pyro_sim.py incompressible_viscous converge inputs.converge.128 vis.dovis=0 - -The error is measured by comparing with the analytic solution using -the routine ``incomp_viscous_converge_error.py`` in ``analysis/``. To -generate the plot below, run - -.. prompt:: bash - - python incompressible_viscous/tests/convergence_errors.py convergence_errors.txt - -or ``convergence_errors_no_limiter.txt`` after running with that option. Then: - -.. prompt:: bash - - python incompressible_viscous/tests/convergence_plot.py - -.. image:: incomp_viscous_converge.png - :align: center - -The solver is converging but below second-order, unlike the inviscid case. Limiting -does not seem to make a difference here. - -Exercises ---------- - -Explorations -^^^^^^^^^^^^ - -* In the lid-driven cavity problem, when does the solution reach a steady-state? - -* :cite:`ghia1982` give benchmark velocities at different Reynolds number for the - lid-driven cavity problem (see their Table I). Do we agree with their results? diff --git a/docs/source/index.rst b/docs/source/index.rst index 283979fe2..688308baa 100644 --- a/docs/source/index.rst +++ b/docs/source/index.rst @@ -61,7 +61,6 @@ new ideas. compressible_basics diffusion_basics incompressible_basics - incompressible_viscous_basics lowmach_basics swe_basics particles_basics diff --git a/docs/source/lowmach_basics.rst b/docs/source/lowmach_basics.rst index 2f647feab..70fa5e5b3 100644 --- a/docs/source/lowmach_basics.rst +++ b/docs/source/lowmach_basics.rst @@ -1,5 +1,6 @@ -Low Mach number hydrodynamics solver -==================================== +***************************** +Low Mach number hydrodynamics +***************************** pyro's low Mach hydrodynamics solver is designed for atmospheric flows. It captures the effects of stratification on a fluid element by @@ -14,6 +15,8 @@ governing equations are: with :math:`\nabla p_0 = \rho_0 g` and :math:`\beta_0 = p_0^{1/\gamma}`. +``lm_atm`` solver +================= As with the incompressible solver, we implement a cell-centered approximate projection method. @@ -21,17 +24,4 @@ The main parameters that affect this solver are: .. include:: lm_atm_defaults.inc -Examples --------- - -bubble -^^^^^^ - -The bubble problem places a buoyant bubble in a stratified atmosphere -and watches the development of the roll-up due to shear as it -rises. This is run as: - -.. prompt:: bash - - pyro_sim.py lm_atm bubble inputs.bubble - +.. include:: lm_atm_problems.inc diff --git a/docs/source/swe_basics.rst b/docs/source/swe_basics.rst index 85537485b..b17c5a97c 100644 --- a/docs/source/swe_basics.rst +++ b/docs/source/swe_basics.rst @@ -1,5 +1,6 @@ -Shallow water solver -==================== +*************************** +Shallow water hydrodynamics +*************************** The (augmented) shallow water equations take the form: @@ -13,26 +14,32 @@ with :math:`h` is the fluid height, :math:`U` the fluid velocity, :math:`g` the gravitational acceleration and :math:`\psi = \psi(x, t)` represents some passive scalar. +``swe`` solver +============== -The implementation here has flattening at shocks and a choice of Riemann solvers. +The pyro ``swe`` implementation has flattening at shocks and a choice of Riemann solvers. The main parameters that affect this solver are: .. include:: swe_defaults.inc -Example problems ----------------- +.. include:: swe_problems.inc -dam -^^^ +Examples +======== -The dam break problem is a standard hydrodynamics problem, analogous to the Sod -shock tube problem in compressible hydrodynamics. It considers a one-multidimensional -problem of two regions of fluid at different heights, initially separated by a dam. -The problem then models the evolution of the system when this dam is removed. -As for the Sod problem, there exists an exact solution for the dam break problem, -so we can check our solution against the exact solutions. See Toro's shallow water -equations book for details on this problem and the exact Riemann solver. +dam +--- + +The dam break problem is a standard hydrodynamics problem, analogous +to the Sod shock tube problem in compressible hydrodynamics. It +considers a one-multidimensional problem of two regions of fluid at +different heights, initially separated by a dam. The problem then +models the evolution of the system when this dam is removed. As for +the Sod problem, there exists an exact solution for the dam break +problem, so we can check our solution against the exact solutions. See +Toro's shallow water equations book for details on this problem and +the exact Riemann solver. Because it is one-dimensional, we run it in narrow domains in the x- or y-directions. It can be run as: @@ -57,40 +64,11 @@ slightly better than the HLLC solver, with less smearing at the shock and head/tail of the rarefaction. -quad -^^^^ - -The quad problem sets up different states in four regions of the -domain and watches the complex interfaces that develop as shocks -interact. This problem has appeared in several places (and a `detailed -investigation -`_ is -online by Pawel Artymowicz). It is run as: - -.. prompt:: bash - - pyro_sim.py swe quad inputs.quad - - -kh -^^ - -The Kelvin-Helmholtz problem models three layers of fluid: two at the top and -bottom of the domain travelling in one direction, one in the central part of the -domain travelling in the opposite direction. At the interface of the layers, -shearing produces the characteristic Kelvin-Helmholtz instabilities, just as -is seen in the standard compressible problem. It is run as: - -.. prompt:: bash - - pyro_sim.py swe kh inputs.kh - - Exercises ---------- +========= Explorations -^^^^^^^^^^^^ +------------ * There are multiple Riemann solvers in the swe algorithm. Run the same problem with the different Riemann solvers @@ -101,7 +79,7 @@ Explorations Extensions -^^^^^^^^^^ +---------- * Limit on the characteristic variables instead of the primitive variables. What changes do you see? (the notes show how to implement diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index d48166307..5c1fd5462 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -14,27 +14,29 @@ * delta, z0, z1: flattening parameters (we use Colella 1990 defaults) -The grid indices look like:: - - j+3/2--+---------+---------+---------+ - | | | | - j+1 _| | | | - | | | | - | | | | - j+1/2--+---------XXXXXXXXXXX---------+ - | X X | - j _| X X | - | X X | - | X X | - j-1/2--+---------XXXXXXXXXXX---------+ - | | | | - j-1 _| | | | - | | | | - | | | | - j-3/2--+---------+---------+---------+ - | | | | | | | - i-1 i i+1 - i-3/2 i-1/2 i+1/2 i+3/2 +The grid indices look like + +:: + + j+3/2--+---------+---------+---------+ + | | | | + j+1 _| | | | + | | | | + | | | | + j+1/2--+---------XXXXXXXXXXX---------+ + | X X | + j _| X X | + | X X | + | X X | + j-1/2--+---------XXXXXXXXXXX---------+ + | | | | + j-1 _| | | | + | | | | + | | | | + j-3/2--+---------+---------+---------+ + | | | | | | | + i-1 i i+1 + i-3/2 i-1/2 i+1/2 i+3/2 We wish to solve @@ -354,38 +356,39 @@ def apply_transverse_flux(U_xl, U_xr, U_yl, U_yr, The states that we represent by indices i,j are shown below (1,2,3,4): - - j+3/2--+----------+----------+----------+ - | | | | - | | | | - j+1 -+ | | | - | | | | - | | | | 1: U_xl[i,j,:] = U - j+1/2--+----------XXXXXXXXXXXX----------+ i-1/2,j,L - | X X | - | X X | - j -+ 1 X 2 X | 2: U_xr[i,j,:] = U - | X X | i-1/2,j,R - | X 4 X | - j-1/2--+----------XXXXXXXXXXXX----------+ - | | 3 | | 3: U_yl[i,j,:] = U - | | | | i,j-1/2,L - j-1 -+ | | | - | | | | - | | | | 4: U_yr[i,j,:] = U - j-3/2--+----------+----------+----------+ i,j-1/2,R - | | | | | | | - i-1 i i+1 - i-3/2 i-1/2 i+1/2 i+3/2 - - - remember that the fluxes are stored on the left edge, so - - F_x[i,j,:] = F_x - i-1/2, j - - F_y[i,j,:] = F_y - i, j-1/2 + :: + + j+3/2--+----------+----------+----------+ + | | | | + | | | | + j+1 -+ | | | + | | | | + | | | | 1: U_xl[i,j,:] = U + j+1/2--+----------XXXXXXXXXXXX----------+ i-1/2,j,L + | X X | + | X X | + j -+ 1 X 2 X | 2: U_xr[i,j,:] = U + | X X | i-1/2,j,R + | X 4 X | + j-1/2--+----------XXXXXXXXXXXX----------+ + | | 3 | | 3: U_yl[i,j,:] = U + | | | | i,j-1/2,L + j-1 -+ | | | + | | | | + | | | | 4: U_yr[i,j,:] = U + j-3/2--+----------+----------+----------+ i,j-1/2,R + | | | | | | | + i-1 i i+1 + i-3/2 i-1/2 i+1/2 i+3/2 + + + remember that the fluxes are stored on the left edge, so:: + + F_x[i,j,:] = F_x + i-1/2, j + + F_y[i,j,:] = F_y + i, j-1/2 Parameters ---------- diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py index dcb9aafbb..721bfc9d1 100644 --- a/pyro/mesh/patch.py +++ b/pyro/mesh/patch.py @@ -943,10 +943,6 @@ def cell_center_data_clone(old): old : CellCenterData2d object The CellCenterData2d object we wish to copy - Note - ---- - It may be that this whole thing can be replaced with a copy.deepcopy() - """ if not isinstance(old, CellCenterData2d):