Setting up your own problem¶
The best way to set up a new problem is to find an existing problem setup that is similar. The main steps are:
Write the initialization script
Write routines for source terms, custom boundary conditions, or other customizations
Write a Riemann solver (if solving a new system of equations)
Write a Makefile if using any custom Fortran code
Write a setplot.py file for visualization
If needed for your problem, custom Riemann solvers, boundary condition routines, source term routines, and other functions can all be written in Python but you may prefer to write some of them in Fortran for performance reasons. The latter approach requires direct use of f2py. See Porting a problem from Clawpack 4.6.x to PyClaw for more details.
Writing the initialization script¶
This script should:
Usually the script then instantiates a Controller, sets the
initial solution and solver, and calls run().
Setting initial conditions¶
Once you have initialized a Solution object, it contains a member state.q whose first dimension is num_eqn and whose remaining dimensions are those of the grid. Now you must set the initial condition. For instance
>>> import numpy as np
>>> Y,X = np.meshgrid(state.grid.y.centers,state.grid.x.centers)
>>> r = np.sqrt(X**2 + Y**2)
>>> width = 0.2
>>> state.q[0,:,:] = (np.abs(r-0.5)<=width)*(1.+np.cos(np.pi*(r-0.5)/width))
>>> state.q[1,:,:] = 0.
>>> state.q[2,:,:] = 0.
Note that in a parallel run we only wish to set the local values of the state
so the appropriate geometry object to use here is the
Grid class.
Setting auxiliary variables¶
If the problem involves coefficients that vary in space or a mapped grid, the required fields are stored in state.aux. In order to use such fields, you must pass the num_aux argument to the State initialization
>>> state = pyclaw.State(domain,solver.num_eqn,num_aux)
The number of fields in state.aux (i.e., the length of its first dimension) is set equal to num_aux. The values of state.aux are set in the same way as those of state.q.
Setting boundary conditions¶
The boundary conditions are specified through solver.bc_lower and
solver.bc_upper, each of which is a list of length solver.num_dim. The
ordering of the boundary conditions in each list is the same as the ordering of
the Dimensions in the Grid; typically \((x,y)\). Thus
solver.bc_lower[0] specifies the boundary condition at the left boundary
and solver.bc_upper[0] specifies the condition at the right boundary.
Similarly, solver.bc_lower[1] and solver.bc_upper[1] specify the
boundary conditions at the top and bottom of the domain.
PyClaw includes the following built-in boundary condition implementations:
pyclaw.BC.periodic- periodic
pyclaw.BC.extrap- zero-order extrapolation
pyclaw.BC.wall- solid wall conditions, assuming that the 2nd/3rd component of q is the normal velocity in x/y.
Other boundary conditions can be implemented by using pyclaw.BC.custom, and
providing a custom BC function. The attribute solver.user_bc_lower/upper must
be set to the corresponding function handle. For instance
>>> def custom_bc(state,dim,t,qbc,num_ghost):
... for i in xrange(num_ghost):
... qbc[0,i,:] = q0
>>> solver.bc_lower[0] = pyclaw.BC.custom
>>> solver.user_bc_lower = custom_bc
If the state.aux array is used, boundary conditions must be set for it
in a similar way, using solver.aux_bc_lower and solver.aux_bc_upper.
Note that although state is passed to the BC routines, they should
NEVER modify state. Rather, they should modify qbc/auxbc.
Setting solver options¶
Using your own Riemann solver¶
The Riemann package has solvers for many hyperbolic systems. If your problem involves a new system, you will need to write your own Riemann solver. A nice example of how to compile and import your own Riemann solver can be seen `here https://github.com/damiansra/empyclaw/tree/master/maxwell_1d_homogeneous`_. You will need to:
Put the Riemann solver in the same directory as your Python script
Write a short makefile calling f2py
import the Riemann solver module in your Python script
Here are some tips if you are converting an old Clawpack 4.5 or earlier Riemann solver:
Rename the file from .f to .f90 and switch to free-format Fortran
Move the spatial index (i) to the last place in all array indexing
Please do contribute your solver to the package by sending a pull request on Github or e-mailing one of the developers. To add your Riemann solver to the Clawpack Riemann package, you will need to:
Place the .f90 file(s) in clawpack/riemann/src.
Add the solver to the list in clawpack/riemann/setup.py
Add the solver to the list in clawpack/riemann/src/python/riemann/setup.py
Add the solver to the list in clawpack/riemann/src/python/riemann/Makefile
Add the solver to the list in clawpack/riemann/src/python/riemann/__init__.py
For very simple problems in one dimension, it may be worthwhile to write the Riemann solver in Python, especially if you are more comfortable with Python than with Fortran. For two-dimensional problems, or one-dimensional problems requiring fine grids (or if you are impatient) the solver should be written in Fortran. The best approach is generally to find a similar solver in the Riemann package and modify it to solve your system.
Adding source terms¶
Non-hyperbolic terms (representing, e.g., reaction or diffusion) can be included in a PyClaw simulation by providing an appropriate function handle to
solver.step_source if using Classic Clawpack. In this case, the function specified should modify q by taking a step on the equation \(q_t = \psi(q)\).
solver.dq_src if using SharpClaw. In this case, the function should return \(\Delta t \cdot \psi(q)\).
For an example, see pyclaw/examples/euler_2d/shockbubble.py.
Setting up the Makefile¶
Generally you can just copy the Makefile from an example in pyclaw/examples and replace the value of RP_SOURCES. Make sure the example you choose has the same dimensionality. Also be sure to use the f-wave targets if your Riemann solver is an f-wave solver.
