python-hydro/pyro2

A framework for hydrodynamics explorations and prototyping

A simple python-based tutorial on computational methods for hydrodynamics

pyro is a computational hydrodynamics code that presents two-dimensional solvers for advection, compressible hydrodynamics, diffusion, incompressible hydrodynamics, and multigrid, all in a finite-volume framework. The code is mainly written in python and is designed with simplicity in mind. The algorithms are written to encourage experimentation and allow for self-learning of these code methods.

The latest version of pyro is always available at:

https://github.com/python-hydro/pyro2

The project webpage, where you'll find documentation, plots, notes, etc. is here:

https://python-hydro.github.io/pyro2

(Note: there is an outdated page at readthedocs.org that is no longer updated)

Table of Contents

• Getting started
• Core Data Structures
• Solvers
• Working with data
• Understanding the algorithms
• Regression and unit testing
• Acknowledgements
• Getting help

Getting started

• By default, we assume python 3.9 or later.

• We require numpy, numba, matplotlib, and h5py for running pyro and setuptools_scm for the install.

• There are several ways to install pyro. The simplest way is via PyPI/pip:

  pip install pyro-hydro
  

  Alternately, you can install directly from source, via

  pip install .
  

  you can optionally add the -e argument to install if you are planning on developing pyro solvers directly.

• Not all matplotlib backends allow for the interactive plotting as pyro is run. One that does is the TkAgg backend. This can be made the default by creating a file ~/.matplotlib/matplotlibrc with the content:

  backend: TkAgg

  You can check what backend is your current default in python via:

  import matplotlib.pyplot
  print matplotlib.pyplot.get_backend()

• If you want to run the unit tests, you need to have pytest installed.

• Finally, you can run a quick test of the advection solver:

  pyro_sim.py advection smooth inputs.smooth
  

  you should see a graphing window pop up with a smooth pulse advecting diagonally through the periodic domain.

Core Data Structures

The main data structures that describe the grid and the data the lives on the grid are described in a jupyter notebook:

https://github.com/python-hydro/pyro2/blob/main/pyro/mesh/mesh-examples.ipynb

Many of the methods here rely on multigrid. The basic multigrid solver is demonstrated in the juputer notebook:

https://github.com/python-hydro/pyro2/blob/main/pyro/multigrid/multigrid-constant-coefficients.ipynb

Solvers

pyro provides the following solvers (all in 2-d):

• advection: a second-order unsplit linear advection solver. This uses characteristic tracing and corner coupling for the prediction of the interface states. This is the basic method to understand hydrodynamics.

• advection_fv4: a fourth-order accurate finite-volume advection solver that uses RK4 time integration.

• advection_nonuniform: a solver for advection with a non-uniform velocity field.

• advection_rk: a second-order unsplit solver for linear advection that uses Runge-Kutta integration instead of characteristic tracing.

• advection_weno: a method-of-lines WENO solver for linear advection.

• 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 states and a 2-shock or HLLC approximate Riemann solver.

• compressible_fv4: a fourth-order accurate finite-volume compressible hydro solver that uses RK4 time integration. This is built from the method of McCourquodale and Colella (2011).

• compressible_rk: a second-order unsplit solver for Euler equations that uses Runge-Kutta integration instead of characteristic tracing.

• compressible_sdc: a fourth-order compressible solver, using spectral-deferred correction (SDC) for the time integration.

• diffusion: a Crank-Nicolson time-discretized solver for the constant-coefficient diffusion equation.

• incompressible: a second-order cell-centered approximate projection method for the incompressible equations of hydrodynamics.

• incompressible_viscous: an extension of the incompressible solver including a diffusion term for viscosity.

• lm_atm: a solver for the equations of low Mach number hydrodynamics for atmospheric flows.

• lm_combustion: (in development) a solver for the equations of low Mach number hydrodynamics for smallscale combustion.

• multigrid: a cell-centered multigrid solver for a constant-coefficient Helmholtz equation, as well as a variable-coefficient Poisson equation (which inherits from the constant-coefficient solver).

• particles: a solver for Lagrangian tracer particles.

• swe: a solver for the shallow water equations.

Working with data

In addition to the main pyro program, there are many analysis tools that we describe here. Note: some problems write a report at the end of the simulation specifying the analysis routines that can be used with their data.

• pyro/util/compare.py: this takes two simulation output files as input and compares zone-by-zone for exact agreement. This is used as part of the regression testing.

  usage: ./compare.py file1 file2

• pyro/plot.py: this takes the an output file as input and plots the data using the solver's dovis method.

  usage: ./plot.py file

• pyro/analysis/

  • dam_compare.py: this takes an output file from the shallow water dam break problem and plots a slice through the domain together with the analytic solution (calculated in the script).

    usage: ./dam_compare.py file

  • gauss_diffusion_compare.py: this is for the diffusion solver's Gaussian diffusion problem. It takes a sequence of output files as arguments, computes the angle-average, and the plots the resulting points over the analytic solution for comparison with the exact result.

    usage: ./gauss_diffusion_compare.py file*

  • incomp_converge_error.py: this is for the incompressible solver's converge problem. This takes a single output file as input and compares the velocity field to the analytic solution, reporting the L2 norm of the error.

    usage: ./incomp_converge_error.py file

  • plotvar.py: this takes a single output file and a variable name and plots the data for that variable.

    usage: ./plotvar.py file variable

  • sedov_compare.py: this takes an output file from the compressible Sedov problem, computes the angle-average profile of the solution and plots it together with the analytic data (read in from cylindrical-sedov.out).

    usage: ./sedov_compare.py file

  • smooth_error.py: this takes an output file from the advection solver's smooth problem and compares to the analytic solution, outputting the L2 norm of the error.

    usage: ./smooth_error.py file

  • sod_compare.py: this takes an output file from the compressible Sod problem and plots a slice through the domain over the analytic solution (read in from sod-exact.out).

    usage: ./sod_compare.py file

Understanding the algorithms

There is a set of notes that describe the background and details of the algorithms that pyro implements:

http://open-astrophysics-bookshelf.github.io/numerical_exercises/

The source for these notes is also available on github:

https://github.com/Open-Astrophysics-Bookshelf/numerical_exercises

Regression and unit testing

The pyro/test.py script will run several of the problems (as well as some stand-alone multigrid tests) and compare the solution to stored benchmarks (in each solver's tests/ subdirectory). The return value of the script is the number of tests that failed.

Unit tests are controlled by pytest and can be run simply via

pytest

Acknowledgements

If you use pyro in a class or workshop, please e-mail us to let us know (we'd like to start listing these on the website).

If pyro was used for a publication, please cite the article found in the CITATION file.

Getting help

We use github discussions as a way to ask about the code:

https://github.com/python-hydro/pyro2/discussions

pyro2's Issues

subclassed ArrayIndexer returns 0-dimensional array instead of scalar for some operations

When we work on an ArrayIndexer object a, and do a.min(), we expect a scalar. NumPy gives a scalar, but our wrapped class returns a 0-dimensional ArrayIndexer object.

This is discussed here at bit:

numpy/numpy#5819

create FaceCenterData class

This should be nodal in one direction or both, e.g. to support MHD B fields

Multigrid projection docs

The sphinx docs (multigrid_basics.rst) have a projection example that is missing a figure. It is meant to be the output of project-periodic.py, but that script is also missing.

There is this script in your hydro_examples repo, but I'm not sure it's doing the same thing?

the limiting for 4th order FV advection is not written

The advection_fv4 method does not implement limiting yet. Thus is only works for smooth problems.

add slotted advection test problem

Add this test problem from the Zalesak paper

projection module?

Should there be a separate projection module that all the low speed solvers can use?

improve energy conservation with gravity

Energy conservation with body forces can be substantially improved if the mass fluxes are used, instead of the cell-centered estimates of the work. For general body forces, I have tested:

       energy_update  = 0.5*dtdy*(Flux_y.v(n=self.vars.idens) + Flux_y.jp(1,n=self.vars.idens)) * \
                         0.5*(old_yacc_src.v() + yacc_src.v())*myg.dy

        energy_update += 0.5*dtdx*(Flux_x.v(n=self.vars.idens) + Flux_x.ip(1,n=self.vars.idens)) * \
                         0.5*(old_xacc_src.v() + xacc_src.v())*myg.dx

        ener.v()[:,:] += energy_update

This is identical to the formulation in Springel 2010 (the Arepo paper).

I also switched the half-step prediction using the primitive variables:

    V_l[1:,:,vars.iu] += 0.5*dt*0.5*(xacc_src.d[0:-1,:]+xacc_src.d[1:,:])
    V_l[1:,:,vars.iv] += 0.5*dt*0.5*(yacc_src.d[0:-1,:]+yacc_src.d[1:,:])

    V_r[1:-1,:,vars.iu] += 0.5*dt*0.5*(xacc_src.d[2:,:]+xacc_src.d[1:-1,:])
    V_r[1:-1,:,vars.iv] += 0.5*dt*0.5*(yacc_src.d[2:,:]+yacc_src.d[1:-1,:])

and likewise for the y-direction. This seems to work well, though there may be a better way to do the indexing for the half-step variables.

Setting PYTHONPATH is unnecessary?

It looks like at least the simple advection problem suggested in the readme runs fine on my Python3.7 setup without first setting the PYTHONPATH environment variable to point at the pyro2 directory. Is it still necessary to set that for some other purposes?

Add units

We should attach units to quantities using the unyt library

store data as ArrayIndexer()

In patch.py, we create the main data as a NumPy array and then case it into an ArrayIndexer() when we do a get_var(). Why not make it ArrayIndexer() from the start?

cannot plot RT data after the simulation

plot.py complains about boundary conditions when run on an RT plotfile

"Latest release" on GitHub release page points at 2014 paper version

If you go to https://github.com/python-hydro/pyro2/releases, you'll see that a tag for the version of the code corresponding to the 2014 paper submission is somehow pinned to be prominently displayed in the releases page. Given the new paper and the several releases that have come out since the 2014 paper, you may want to consider adjusting the repository settings so that newer releases are visible without clicking a button on the page.

create a mapped grids class and solvers

A simple MOL advection solver on mapped grids would be a first start

add python 3.7 to testing

our travis only considers 3.6.

fourth order source terms

The gravity source term for compressible_fv4 is not treated as 4th order

well balanced scheme doesn't work anymore?

It seems that the well-balanced scheme doesn't work anymore. If I run:

./pyro.py compressible_rk hse inputs.hse vis.dovis=1 compressible.well_balanced=1

it goes 16 steps, and has all NaNs in the velocity in the output

automate docs building

at the moment, the docs live on the gh-pages branch, but must be manually pushed. We should automate the docs build

create MOL MHD solver

We should create a simple MHD solver using constrained transport and a method of lines integration. This will need:

• Face-centered data (#63)
• a Riemann solver (#67)
• a MHD solver that inherits from compressible

create a nodal MG solver?

This could have applications to the projection stuff

logo

We need a vector graphics pyro logo

add particle support to pyro

We want a Particles class that can store and manage particles and update their positions based on the velocity on the grid.

Add _skeleton directory to help make new solvers/problems

To help creators of new solvers/problems, we could have a directory containing files with outlines of the functions required, similar to yt's frontends/_skeleton directory

HSE BC issue with compressible_rk?

The RT problem seems to get pushed up from the boundary when run with compressible_rk. Other hydro problems are fine, so this is a gravity or HSE BC problem.

move to standalone organization

We should move pyro to it's own organization

ordering of component differs in `Grid2d` and `CellCenterData2d`

In mesh.patch, we allocate the state data as (nvar, nx, ny), but in the compressible unsplit fluxes routines, we have it swapped, as (nx, ny, nvar) (first constructed in interface_f.states()).

This is also the ordering that Grid2d.scratch_array() does (nvar at the end).

We need to make all of this consistence.

switch ArrayIndexer over to subclassing np.adarray

This will clean-up the code substantially, since I won't need to reimplement all the ufuncs.

A sample implementation of this is already done in the rk branch

fix pyro binder

the pyro binder launch is not working

create a JOSS paper

There have been a lot of updates since the original pyro paper, and more importantly, more contributors, so it is a good time for a new pub, this time via JOSS.

create a Pyro class

At the moment, we run a simulation through the script pyro.py. It would be nice to be able to interact better in the Jupyter notebook. To make this easier, it would be useful to move most of the driver loop in pyro.py into a new class, Pyro, and have evolve methods that let us evolve the whole simulation, or perhaps just take a single step. This would be useful for interactive explorations in notebooks

Consider making pyro2 pip-installable

This is unrelated to my JOSS review so don't take this as a blocker for paper acceptance.

You may want to consider making pyro a pip-installable package, rather than relying on the directory structure of the git repository to make it runnable. This will make it easier for users to install pyro (people could just do pip install pyro2 or conda install -c conda-forge pyro2) and also make it easier for developers who have worked with other python packages to get set up initially. If you did this you could make the pyro.py script an executable that gets installed to a user's UNIX PATH out of the box, so there would be no need to always work in the pyro2 directory or refer to the pyro.py executable using its full path.

switch to Fortran order for arrays?

Since we are interfacing with Fortran, then it is more natural to use ordering="F" with our NumPy arrays. We like to have "x" as the first index (row) and "y" as the second index (column).

Also, this will eliminate a lot of transposes for plotting and the pretty_printing

and when we just print an ndarray to the screen, x will vary horizontally instead of vertically.

tutorials

we should have some tutorials

add Burgers solver

Both viscous and inviscid

documentation

Need to better document the flow of the code, including:

-- how to add custom BCs
-- how to add solvers
-- how to add new problems

AttributeError: 'list' object has no attribute 'copy'

I'm getting this error when trying to simulate the advection for advection_rk,fv4,weno solvers.

File "/home/administrator/Downloads/pyro2-master/mesh/patch.py", line 848, in cell_center_data_clone
new.derives = old.derives.copy()
AttributeError: 'list' object has no attribute 'copy'

I changed the "old.derives.copy()" to just "old.derives" the solver ran normally.

I'm using pyro2 in python 2.7 environment.

advection_nonuniform is not documented

We should add a section about this solver to the advection docs

MG solvers should check that incoming coefficient array has right # of ghost cells

Right now we don't check that the coefficient array has the same # of ghost cells in the variable_coefficient_MG solver. We should create a CCData2d object for the coefficients and pass that through, like we do with the general solver.

But then we should also make sure that we check the number of ghost cells before copying.

initial iterations in `incompressible` and `lm_atm`

These solvers call method_compute_timestep to get the intiial timestep for the initial iterations, but this doesn't have the safety checks that the main compute_timestep provides -- should they call the main version?

Also, since the pressure is lagged as p^{n-1/2}, is there a correction needed to account for the fact that dt is not the same between the previous step as it is for the next? (sources don't seem to deal with this)

y-velocity and energy are same just before plotting

Pardon me if this is silly,
If I place a breakpoint just before line 284 in simulation.py of compressible solver, to check the variables u, v, p, rho, e..
the arrays of v and energy are same values.

And also ivars.iu = ivars.iener = 1

make python3 the default

switch mk.sh and pyro.py to default to python 3

lm_combustion should subclass simulation_null

there is a lot of redundancy in the lm_combustion solver -- it should subclass simulation_null like the others do

pickle doesn't work for external functions

We save by pickling, but with inhomogeneous BCs in the MG solvers, we pass a function reference into the BC object and that gets pickled. When we unpickle, it is invalid and throws an error.

get rid of wildcard imports

there are some places where we do

import compressible.unsplitfluxes as *

e.g. We should not do this...

add reactive flow

There is a start of such a solver as compressible_react using Strang splitting

reflecting boundary conditions do not ensure zero mass flux

Using reflecting boundary conditions in a (nearly) isothermal exponential atmosphere, I find that pyro does not ensure a zero mass flux at the lower boundary, whether or not I include gravity, and whether or not I reflect the gravity vector at the boundary (it should, of course, be reflected to match the boundary conditions, and I do so by adding a new gravitational acceleration state variable and setting its boundary condition to reflect y-odd).

I have tested this setup in Castro and Athena, which both produce identically zero mass flux at the reflecting boundary.

The initial conditions are prescribed by the following:

xmom.d[:,:] = 0.0
ymom.d[:,:] = 0.0
dens.d[:,:] = 0.0
grav.d[:,:] = grav_const # this will be reflected via the boundary conditions routines              

# set the density to be stratified in the y-direction                                               
myg = my_data.grid

p = myg.scratch_array()

dens.d[:,:] = dens1*numpy.exp(-myg.y2d/scale_height)
p.d[:,:] = dens.d*cs*cs/gamma

# set the energy (P = cs2*dens)                                                                     
ener.d[:,:] = p.d[:,:]/(gamma - 1.0) + \
    0.5*(xmom.d[:,:]**2 + ymom.d[:,:]**2)/dens.d[:,:]

compare.py should check errors to some tolerance

compare.py currently checks the outputs of benchmark files are exactly equal. Instead, it should check that they agree up to some tolerance.

This will help with the new Numba implementation as the results now seem to be more sensitive to the exact system setup/compiler.

consider migrating from f2py to Cython

back when pyro started, there was no cython. Now we should consider migrating to cython cause it is easier and can be used more easily on some platforms (e.g. Macs)

create MHD HLLD Riemann solver

Consider adding support for pyro to yt

Again, not a blocker for accepting the pyro paper.

It would be cool to add support for pyro's output format to yt. If you ever wanted to add AMR support to pyro this would make the visualization and analysis problem much easier. You might also be able to do cool in-memory analysis and visualization stuff by passing data back and forth between pyro and yt via the python process.

stored benchmarks work in python2 but not python3

pyro simply uses pickle to store the output, since this is easy to implement, but it is not portable. The test benchmarks were output with python3, and they are in its pickle format. As a result, python2 cannot unpickle them. This discusses some of the issues:

http://stackoverflow.com/questions/28218466/unpickling-a-python-2-object-with-python-3

Potential solutions are

• convert to h5py to write out a portable output file. This is probably best in the long run, but will add another dependency and require a read routine and reinitialization routine to give the features that unpickle does for us

• have 2 benchmarks directories, 2/ and 3/, storing the python2 and python3 output respectively.