The Clawpack 5.X Software

1 Clawpack is a software package designed to solve nonlinear hyperbolic partial differential equa2 tions using high-resolution finite volume methods based on Riemann solvers and limiters. The pack3 age includes a number of variants aimed at different applications and user communities. Clawpack 4 has been actively developed as an open source project for over 20 years. The latest major release, 5 Clawpack 5, introduces a number of new features and changes to the code base and a new devel6 opment model based on GitHub and Git submodules. This article provides a summary of the most 7 significant changes, the rationale behind some of these changes, and a description of our current 8 development model. 9


Introduction
The Clawpack software suite [12] is designed for the solution of nonlinear conservation laws, balance laws, and other first-order hyperbolic partial differential equations not necessarily in conservation form.
The underlying solvers are based on the wave propagation algorithms described by LeVeque in [34], and are designed for logically Cartesian uniform or mapped grids or an adaptive hierarchy of such grids.
The original Clawpack was first released as a software package in 1994 and since then has made major strides in both capability and interface.More recently a major refactoring of the code and a move to GitHub for development has resulted in the release of Clawpack 5.0 in January, 2014.A significant number of additional improvements have been made since then.
Because scientific software has become central to many advances made in science, engineering, resource management, natural hazards modeling and other fields, it is increasingly important to describe and document changes made to widely used packages.Such documentation efforts serve to orient new and existing users to the strategies taken by developers of the software, place the software package in the context of other packages, document major code changes, and provide a concrete, citable reference for users of the software.
With this in mind, the goals of this paper are to: • Summarize the development history of Clawpack, • Summarize some of the major changes made between the early Clawpack 4.x versions and the most recent version, Clawpack 5.3, • Summarize the development model we have adopted, for managing open source scientific software projects with many contributors, and • Identify how users can contribute to the Clawpack suite of tools.
This paper provides a brief history of Clawpack in Section 1.1, a background of the mathematical concerns in Section 1.2, the modern development approach now being used in Section 2, the major feature additions in the Clawpack 5.x major release up until Version 5.3 in Section 3. Some concluding thoughts and future plans for Clawpack are mentioned in Section 4.

History of Clawpack
The first version of Clawpack was released by LeVeque in 1994 [32] and consisted of Fortran code for solving problems on a single, uniform Cartesian grid in one or two space dimensions, together with some Matlab [40] scripts for plotting solutions.The wave-propagation method implemented in this code provided a general way to apply recently developed high-resolution shock capturing methods to general hyperbolic systems and required only that the user provide a "Riemann solver" to specify a new hyperbolic problem.Collaboration with Berger [7] soon led to the incorporation of adaptive mesh refinement (AMR) in two space dimensions, and work with Langseth [31,30] led to three-dimensional versions of the wave-propagation algorithm and the software, with three-dimensional AMR then added by Berger.Version 4.3 of Clawpack contained a number of other improvements to the code and formed the basis for the examples presented in a textbook [34] published in 2003.That text not only provided a complete description of the wave propagation algorithm, developed by LeVeque, but also is notable in that the codes used to produce virtually all of figures in the text were made available online [34].These examples are available at http://depts.washington.edu/clawpack/clawpack-4.3/book.html.
In 2009, Clawpack Version 4.4 was released with a major change from Matlab to Python as the recommended visualization tool, and the development of a Python user interface for specifying the input data.
Since then, a number of other features were added to handle new applications, to provide a better user interface and visualization tools, to incorporate higher-order accurate algorithms, to parallelize through MPI and OpenMP, and other enhancements.The Clawpack 4.x line of code ended with Version 4.6.3(released in January 2013) with many of the changes from 4.3 to 4.6 1 .Version 5 of Clawpack introduces a number of modern approaches to code development, interfacing with other codes, and adding new capabilities.These changes are the subject of the rest of this paper.

Hyperbolic problems
In one space dimension, the hyperbolic systems solved with Clawpack typically take the form of conservation laws or non-conservative linear systems where subscripts denote partial derivatives and q(x, t) is a vector with m ≥ 1 components.The coefficient matrix A in (2) or the Jacobian matrix f (q) in ( 1) is assumed to be diagonalizable with real eigenvalues for all relevant values of q, x, and t.This condition guarantees that the system is hyperbolic, with solutions that are wave-like.The eigenvectors of the system determine the relation between the different components of the system, or waves, and the eigenvalues determine the speeds at which these waves travel.The right hand side of these equations could be replaced by a "source term" ψ(q, x, t) to give a non-homogeneous equation that is sometimes called a "balance law" rather than a conservation law.Spatially-varying flux functions f (q, x) in (1) can also be handled using the f-wave approach [3].
Examples of equations solved by Clawpack include: • Advection equation(s) for one or more tracers.The velocity field is typically prescribed from the solution to another fluid flow problem, such as wind.Typical applications include transport of heat, energy, pollution, smoke or another species that does not influence the velocity field.
• The Euler equations of compressible, inviscid fluid dynamics, consist of conservation laws for mass, momentum, and energy.The wave speeds depend on the local fluid velocity and the acoustic wave velocity (sound speed).Source terms can be added to include the effect of gravity, viscosity or heat transfer.These systems have important applications in aerodynamics, climate and weather modeling, and astrophysics.
• The shallow water equations, describing the velocity and surface height of a liquid whose depth is small relative to typical wavelengths.In this case source terms may include the effect of varying bathymetry and of bottom friction.These equations are used, for instance, to model inundation caused tsunamis and dam breaks, as well as atmospheric flows.
• Elastic wave equations, used to model waves in solid materials.Here even a linear problem can be complex due to varying material properties on multiple scales that then effect the wave speeds.
Discontinuities (shock waves) can arise in the solution of nonlinear hyperbolic equations, causing difficulties for traditional numerical methods based on discretizing derivatives directly.Modern shock capturing methods are often based on solutions to the Riemann problem that consists of equations (1) or (2) together with piecewise constant initial data with a single jump discontinuity.The solution to the Riemann problem is a similarity solution (a function of x/t only), typically consisting of m waves (for a system of m equations) propagating at constant speed.This is true even for nonlinear problems, where the waves may be shocks or rarefaction waves (through which the solution varies continuously in a self-similar manner).
The main theoretical and numerical difficulties of hyperbolic problems involve the prescription of physically correct weak solutions and understanding the behavior of the solution at discontinuities.The Riemann solver is an algorithm that encodes the specifics of the hyperbolic system to be solved, and it is the only routine (other than problem-specific setup such as initial conditions) that needs to be changed in order to apply the code to different hyperbolic systems.In some cases, the Riemann solver may also be designed to enforce physical properties like positivity (e.g., for the water depth in GeoClaw) or to account for forces (like that of gravity) that may be balanced by flux terms.
Clawpack is based on Godunov-type finite volume methods in which the solution is represented by cell averages.Riemann problems between the cell averages in neighboring states are used as the fundamental building block of the algorithm.The wave-propagation algorithm originally implemented in Clawpack (and still used in much of the code) is based on using the waves resulting from each Riemann solution together with limiter functions to achieve second-order accuracy where the solution is smooth together with sharp resolution of discontinuities without spurious numerical oscillations (see [34] for a detailed description of the algorithms).The recently developed SharpClaw algorithms (see Section 3.6), now incorporated into PyClaw, use higher-order WENO methods but rely on the same Riemann solvers.
Problem-specific boundary conditions must also be imposed, which are implemented by a subroutine that sets the solution value in ghost cells exterior to the domain each time step.The Clawpack software contains library routines that implement several sets of boundary conditions that are commonly used, e.g.periodic boundary conditions, reflecting solid wall boundary conditions for problems such as acoustics, Euler, or shallow water equations, and non-reflecting (absorbing) extrapolation boundary conditions.As with all Clawpack library routines, the boundary condition routine can be copied and modified by the user to implement other boundary conditions needed for a particular application.
In two space dimensions, hyperbolic equations might take the form q t (x, y, t) + f (q(x, y, t)) x + g(q(x, y, t)) y = 0 (3) or q t (x, y, t) + A(x, y)q(x, y, t) x + B(x, y)q(x, y, t) y = 0 (4) In order to be hyperbolic, the coefficient matrices A and B in (3) or the Jacobian matrices f (q) and g (q) in (4) must have the property that any linear combination gives a diagonalizable matrix with real eigenvalues.The extension to three space dimensions is similar.
In two or three space dimensions, the wave-propagation methods are extended using either dimensional splitting, so that only one-dimensional Riemann solvers are needed, or by a multi-dimensional algorithm based on transverse Riemann solvers introduced in [33].Both approaches are supported in Clawpack.A variety of Riemann solvers have been developed for Clawpack, many of which are collected in the riemann repository, see Section 3.2.
Adaptive mesh refinement is essential for many problems and has been available in two space dimensions since 1995, when Marsha Berger joined the project team and her AMR code for the Euler equations of compressible flow was generalized to fit into the software which became AMRClaw [8].
AMRClaw was carried over to three space dimensions using the unsplit algorithms introduced in [31].
Starting in Version 5.3.0,dimensional splitting is also supported in AMRClaw, which can be particularly useful in three space dimensions where the unsplit algorithms are much more expensive.Other recent improvements to AMRClaw are discussed in Section 3.4.

Development Approach
Clawpack's development model is driven by the needs of its developer community.The Clawpack project consists of several interdependent projects: core solver functionality, a visualization suite, a general adaptive mesh refinement code, a specialized geophysical flow code, and a massively parallel Python framework.Changes to the core solvers and visualization suite have a downstream effect on the other codes, and the developers largely work in an independent, asynchronous manner across continents and time zones.
The core Clawpack software repositories are: • clawpack -responsible for installation and coordination of other repositories, • riemann -Riemann solvers used by all the other projects, • visclaw -a visualization suite used by all the other projects, • clawutil -utility functions used by most other projects, • classic -the original single grid methods in 1, 2, and 3 space dimensions, • amrclaw -the general adaptive mesh refinement framework in 2 and 3 dimensions, • geoclaw -solvers for depth-averaged geophysical flows which employs the framework in amrclaw,

and
• pyclaw -a Python implementation and interface to the Clawpack algorithms including highorder methods and massively parallel capabilities.
A release of Clawpack downloaded by users contains all of the above.The repositories riemann, visclaw, and clawutil are sometimes referred to as upstream projects, since their changes affect all the remaining projects in the above list, commonly referred to as downstream projects.There are some variations on this, for instance AMRClaw is upstream of GeoClaw, which uses many of the algorithms and software base from AMRClaw.To coordinate this the clawpack repository points to the compatible version of each repository, described later in this section.
Beyond the major code containing repositories additional repositories contain documentation and extended examples for using the packages: • doc -the primary documentation source files, developed using Sphinx2 , • clawpack.github.com-a host repository for the documentation html files that appear at http: //www.clawpack.org,and • apps -applications contributed by developers and users that go beyond the introductory examples included in the core repositories.
The Clawpack 4.x code is also available in the repository clawpack-4.xbut is no longer under development.

Version Control
The Clawpack team uses the Git distributed version control system to coordinate development of each major project.The repositories are publicly coordinated under the Clawpack organization on GitHub3 with the top-level clawpack super-repository responsible for hosting build and installation tools, as well as providing a synchronization point for the other repositories.The remaining "core Clawpack repositories" listed above are subrepositories of the main clawpack organization.
GitHub itself is a free provider of public Git repositories.In addition to repository hosting, the Clawpack team uses GitHub for issue tracking, code review, automated continuous integration via Travis CI 4 , and test coverage tracking via Coveralls5 for the Python-based modules.The issue tracker on GitHub supports cross-repository references, simplifying communication between Clawpack developer sub-teams.The Travis CI service, which provides free continuous integration for publicly developed repositories on GitHub, runs Clawpack's test suites through nose6 on proposed changes to the code base, and through a connection to the Coveralls service, reports on any test failures as well as changes to test coverage.

Submodules
The clawpack "super-repository" serves as an installation and synchronization point for the project repositories: each of the other core Clawpack repositories listed above is a submodule of the clawpack repository.A commit to the clawpack repository serves to keep track of the versions of each submodule that are meant to function together.Git submodules provide an invaluable mechanism for allowing Clawpack team members to work asynchronously on independent projects while reusing and maintaining common software infrastructure.
Typically the Clawpack developers advance the master development branch of the top-level clawpack repository any time a major feature is added or a bug is fixed in one of the upstream projects that might affect code in other repositories.By checking out a particular commit in the clawpack repository and performing a git submodule update, all repositories can be updated to versions that are intended to be consistent and functional.
In particular, when Travis runs the regression tests in any project repository (performed automatically for any pull request), it starts by installing Clawpack on a virtual machine and the current head of the clawpack/master branch indicates the commit from each of the other projects that must be checked out before performing the tests.If the clawpack repository has not been properly updated following changes in other upstream projects, these tests may fail.
Any new release of Clawpack is a snapshot of one particular commit on clawpack and the related commits on all submodules.These particular commits are also tagged for future reference with consistent names, such as v5.3.1.(Git tags simply provide a descriptive name for a particular commit rather than having to refer to a Git hash code.)

Contributing
Scientists who program are often discouraged from sharing code due to existing reward mechanisms and the fear of being "scooped".However, recent studies indicate that scientific communities that openly share and develop code have an advantage because each researcher can leverage the work of many others [43], and that paper citation rates can be increased by sharing code [44] and/or data [41].
Moreover, journals and funding agencies are increasingly requiring investigators to share code used to obtain published results.One of the goals of the Clawpack project is to facilitate code sharing by users, by providing an easy mechanism to refer to a specific version of the Clawpack software and ensuring that past versions of the software remain available on a stable and citable platform.
On the development side, we expect that the open source development model with important discussions conducted in public will lead to further growth of the developer community and additional contributions from users.Over the past twenty years, many users have written code extending Clawpack with new Riemann solvers, algorithms, or domain-specific problem tools.Unfortunately, much of this code did not make it back into the core software for others to use.Many of the development changes in Clawpack 5.x were done to encourage contributions from a broader community.We have begun to see an increase in contributions from outside the developers' groups, and hope to encourage more of this in the future.
The primary development model is typical for GitHub projects: a contributor forks the repository on GitHub, then develops improvements in a branch that is pushed to her own fork.She issues a "pull request" (PR) when the branch is ready to be merged into the main repository.Increasingly, contributors are also using PRs as a way to conveniently post preliminary or prototype code for discussion prior to further development, often marked WIP for "work in progress" to signal that it is not ready to merge.
After a PR is issued, other developers, including one or more of the maintainers for the corresponding project, review the code.The Travis CI server also automatically runs the tests on the proposed new code.The test results are visible on the GitHub page for the PR.Usually there is some iteration as developers suggest improvements or discuss implementation choices in the code.Once the tests are passing and it is agreed that the code is acceptable, a maintainer merges it.

Releases
Although Clawpack is continuously developed, it is convenient for users to be able to install stable versions of the software.The Clawpack developers provide these releases through two distribution channels: GitHub and the Python Package Index (PyPI).Full source releases are available on GitHub.
Alternatively, the PyClaw subproject and its dependencies can be installed automatically using a PyPI client such as pip.
Clawpack does not follow a calendar release cycle.Instead, releases emerge when the developer community feels enough changes have accumulated since the last release to justify the cost of switching to a new release.For the most part, Clawpack releases are versioned using an M.m.p triplet, representing the major (M), minor (m), and patch (p) versions respectively.In the broader software engineering community, this is often referred to as semantic versioning.Small changes that fix bugs and cosmetic issues result in increments to the patch-level.Backwards-compatible changes result in an increase to the minor version.The introduction of backwards-incompatible changes require that the major version be incremented.In addition, the implementation of significant new algorithms or capability will also justify the increment of major release number, and is often an impetus for providing another release to the public.In practice, the Clawpack software has frequently included changes in minor version releases that were not entirely backwards compatible, but these have been relatively minor and documented in the release notes.Major version numbers have changed infrequently and related to major refactoring of the code as in going from 4.x to 5.0.

Advances
This section describes the major changes in each of the code repositories in moving from Clawpack Clawpack Git repositories for a more complete list.

Global Changes
Substantial redesign of the Clawpack code base was performed in the move from Clawpack 4.x to 5.x.Major changes that affected all aspects of the code include: • The interface to the Clawpack Riemann solvers was changed so that one set of solvers can be used for all versions of the code (including PyClaw via f2py 8 ).Rather than appearing in scattered example directories, these Riemann solvers have all been collected into the new riemann repository.Modifications to the calling sequences were made to accommodate this increased generality.
• Calling sequences for a number of other Fortran subroutines were also modified based on experiences with the Clawpack 4.x code.These can also be used as a stand-alone product for those who only want the Riemann solvers.
• Python front-ends were redesigned to more easily specify run-time options for the solver and visualization.The Fortran variants (ClassicClaw, AMRClaw, and GeoClaw) all use a Python script to facilitate setting input variables.These scripts create text files with a rigidly specified format that are then read in when the Fortran code is run.The interface now allows updates to the input parameters while maintaining backwards compatibility.
• The indices of the primary conserved quantities were reordered.In Clawpack 4.x, the mth component of a system of equations in grid cell (i, j) (in two dimensions, for example), was stored in q(i,j,m).In order to improve cache usage and to more easily interface with PETSc, a global change was made to the ordering so that the component number comes first; i.e. q(m,i,j).A seemingly minor change like this affects a huge number of lines in the code and cannot easily be automated.The use of version control and regression tests was crucial in the successful completion of the project.compute the waves (or discontinuities) that make up the Riemann solution.In the unsplit algorithm,

Riemann
Clawpack also makes use of transverse Riemann solvers, responsible for computing transport between cells that are only corner (in 2d) or edge (in 3d) adjacent.
For nonlinear systems, the exact solution of the Riemann problem is computationally costly and may involve both discontinuities (shocks and contact waves) and rarefactions.It is almost always preferable to employ inexact Riemann solvers that approximate the solution using discontinuities only, with an appropriate entropy condition.The solvers available in Clawpack are all approximate solvers, although one could easily implement their own exact solver and make it available in the format needed by Clawpack routines.
A common feature in all packages in the Clawpack suite is the use of a standard interface for Fortran Riemann solver routines.This ensures that new solvers or solver improvements developed for one package can immediately be used by all packages.To further facilitate this sharing and to avoid duplication, Riemann solvers are (with rare exceptions) not maintained under the other packages but are collected in a single repository named riemann.Users who develop new solvers for Clawpack are encouraged to submit them to the Riemann repository.
In the Fortran-based packages (Classic, AMRClaw, and GeoClaw) the Riemann solver is selected at compile-time by modifying a problem-specific Makefile.In PyClaw, the Riemann solver to be used is selected at run-time.This is made possible by compiling all of the Riemann solvers (when PyClaw is installed) and generating Python wrappers with f2py.For PyClaw, riemann also provides metadata (such as the number of equations, the number of waves, and the names of the conserved quantities) for each solver so that setup is made more transparent.

ClassicClaw
The classic repository contains code implementing the wave propagation algorithm on a single uniform grid, in much the same form as the original Clawpack 1.0 version of 1994 but with various enhancements added through the years.Following the introduction of Clawpack 4.4 the three-dimensional routines were left out of the Python user interfaces and plotting routines.These have been reintroduced in Clawpack 5. Additionally the OpenMP shared-memory parallelism capabilities have been extended to the three-dimensional code.

AMRClaw
Fortran code in the AMRClaw repository performs block-structured adaptive mesh refinement [4,5] for both Clawpack and GeoClaw applications.The algorithms implemented in AMRClaw are discussed in detail in [7,36], but a short overview is given here to set the stage for a description of recent changes.AMRClaw includes the functionality for: • Coordinating the flagging of points where refinement is needed, with a variety of criteria possible for flagging cells that need refinement from each level to the next finer level (including Richardson extrapolation, gradient testing, or user-specified criteria) 9 , • Organizing the flagged points into efficient grid patches at the next finer level, using the algorithm of [9], • Interpolating the solution to newly created fine grids and initializing auxiliary data (topography, wind velocity, metric data and so on) on these grids, • Averaging fine grid solutions to coarser grids, • Orchestrating the adaptive time stepping (i.e.sub-cycling in time), • Interpolating coarse grid solution to fine grid ghost cells, and 9 See http://www.clawpack.org/flag.html • Maintaining conservation at patch boundaries between resolution levels.
AMRClaw now allows users to specify "regions" in space-time ] in which refinement is forced to be at least at some level L 1 and is allowed to be at most L 2 .This can be useful for constraining refinement, e.g.allowing or ensuring resolution of only a small coastal region in a global tsunami simulation.Previously the user could enforce such conditions by writing a custom flagging routine, but now this is handled in a general manner so that the parameters above can all be specified in the Python problem specification.Multiple regions can be specified, and a simple rule is used to determine the constraints at a grid cell that lies in multiple regions.
Auxiliary arrays are often used in Clawpack to store data that describes the problem and the routine.The routine setaux must then be provided by the user to set these values each time a new grid patch is created.For some applications computing these values can be time-consuming.In Clawpack 5.2, this code was improved to allow reuse of values from previous patches at the same level where possible at each regridding time.This is backward compatible, since no harm is done if previously written routines are used that still compute and overwrite instead of checking a mask.
In Clawpack 5.3 the capability to specify spatially varying boundary conditions was added.For a single grid, it is a simple matter to compute the location of the ghost cells that extend outside the computational domain and set them appropriately.With AMR however, the boundary condition routine can be called for a grid located anywhere in the domain, and may contain fewer or larger numbers of ghost cells.For this reason, the boundary condition routines do not assume a fixed number of ghost cells.
Anisotropic refinement is allowed in both two and three dimensions.This means that the spatial and temporal refinement ratios can be specified independently from one another (as long as the temporal refinement satisfies the CFL condition).In addition, capabilities have been added to automatically select the refinement ratio in time on each level based on the CFL condition.This has only been implemented in GeoClaw.where the wave speed in the shallow water equations depends on the local depth.The finest grids are often located only in shallow coastal regions, so a large refinement ratio in space does not lead to a large refinement ratio in time.
AMRClaw has been parallelized using OpenMP directives using a patch-based decomposition.The main paradigm in structured AMR is a loop over all patches at a given level, where some operation is done on each patch (i.e.taking a time step, finding ghost cells, conservation updates, etc.).This lends itself easily to a parallel for loop construct where each iteration of the loop corresponds to a grid at that level.Dynamic scheduling is used with a chunk size of one, so that one thread is assigned one patch at a time.To help with load balancing, patches at each level are sorted from largest to smallest workload when they are first created, using the total number of cells in the grid as an indicator of work.
Note that this approach causes a memory bulge.Each thread must have its own scratch arrays to save the incoming and outgoing waves and fluxes for future conservation fix-ups.The bulge is directly proportional to the number of threads executing.For stack-based memory allocation per thread, the use of the environment variable OMP STACKSIZE to increase the limit may be necessary.
Figure 1 shows two snapshots of the solution to a three-dimensional shock-bubble interaction problem found in the Clawpack apps repository, illustrating localized phenomena requiring adaptive refinement.In Fig. 2 we show scalability tests and some timings for this example, when run on a 40 core Intel Xeon Haswell machine, using KMP AFFINITY compact.For timing purposes, the only modifications made to the input parameters was to turn off checkpointing and graphics output.The plot on the left shows that most of the wall clock time is in the integration routine (stepgrid), which closely tracks the total time.The second chunk of time is in the regridding, which contains algorithms that are not completely scalable.Very little time is in the filling of ghost cells, mostly from other patches but also includes those at domain boundaries.The efficiency is above 80% until 24 cores, then drops off dramatically.Note that there are only 2 level 1 grids, and an average of 22.8 level 2 grids.Most of the work is on level 3 grids, where there are an average of 138.1 grids over all the level 3 timestep.This is very coarse for large numbers of cores (hence the dropoff in efficiency).At 40 cores, there are less than 4 grids per core, and the grids are very different sizes.The parallelization of AMRClaw and GeoClaw assumes multi-core machines for the target architecture.PyClaw, on the other hand, does not include AMR but uses MPI via PETSc to achieve parallelism on distributed memory machines that scale to tens of thousands of cores (see Section 3.6).
Other frameworks exist, most notably ForestClaw [11], which are being developed in parallel with AMRClaw, that provide scalable AMR calculations on large distributed memory machines.

GeoClaw
The GeoClaw branch of Clawpack was developed to solve the two-dimensional shallow water equations over topography for modeling tsunami generation, propagation, and inundation.The AMRClaw code formed the starting point but it was necessary to make many modifications to support the requirements of this application, as described briefly below.This code originated with the work of George [15,16,17] and was initially called TsunamiClaw.Later it became clear that many other geophysical flow applications have similar requirements and the code was generalized as GeoClaw.
One of the major issues is the treatment of wetting and drying of grid cells at the margins of the flow.The handling of dry states in a Riemann solver is difficult to handle robustly, and has gone through several iterations.GeoClaw must also be well-balanced in order to preserve steady states, in particular the "ocean at rest".To achieve this, the source terms in the momentum equations arising from variations in topography are incorporated into the Riemann solver rather than using a fractional step splitting approach.This is critical for modeling waves that have very small amplitudes relative to the variations in the depth of the ocean.See [35] for a general discussion of such methods and [16,17] for details of the Riemann solver used in GeoClaw.Other features of GeoClaw include the ability to solve the equations in latitude-longitude coordinates on the surface of the sphere, and the incorporation of source terms modeling bottom friction using a Manning formulation.More details about the code and tsunami modeling applications can be found in [6,36].In 2011, a significant effort took place to verify and validate GeoClaw against the US National Tsunami Hazard Mitigation Program (NTHMP) benchmarks [21].NTHMP approval of the code allows GeoClaw to be used in hazard mapping projects that are funded by this program or other federal and state agencies, e.g.[19,20].One such project is illustrated in Fig. 3.
In addition to a variety of tsunami modeling applications, GeoClaw has been used to solve dam break problems in steep terrain [14], storm surge problems [39] (see Fig. 4), and submarine landslides [29].The code also formed the basis for solving the multi-layer shallow water equations for storm surge modeling [37,38], and is currently being extended further to handle debris flow modeling in the packages D-Claw [24,18] (see Fig. 5).
Nearly one quarter of the files in the AMRClaw source library have to be modified for GeoClaw.
There are currently 113 files in the AMRClaw 2D library, of which 26 are replaced by a GeoClawspecific files of the same name in the GeoClaw 2D library.For example, to preserve a flat sea surface when interpolating, it is necessary to interpolate the surface elevation (topography plus water depth) rather than simply interpolating the depth component of the solution vector as would normally be done in AMRClaw.An additional 24 files in the GeoClaw shallow water equations library handle other complications introduced by the need to model tsunamis and storm surge.
Several other substantial improvements in the algorithms implemented in GeoClaw have been made between versions 4.6 and 5.3.0,including: • In depth-averaged flow, the wave speed and therefore the CFL condition depends on the depth.
As a result, flows in shallow water that have been refined spatially may not need to be refined in time.This "variable-time-stepping" was easily added along with the anisotropic capabilities that were added to AMRClaw.
• The ability to specify topography via a set of topo files that may cover overlapping regions at different resolutions has been added.The finite volume method requires cell averages of topography, computed by integrating a piecewise bilinear function constructed from the input topo files over each grid cell.In Clawpack 5.1.0,this was improved to allow an arbitrary number of nested topo grids.When adaptive mesh refinement is used, regridding may take place every few time steps.Improvements were made in 5.2.0 so that topography could be copied rather than always being recomputed in regions where there is an existing old grid.
• The user can now provide multiple dtopo files that specify changes to the initial topography at a series of times.This is used to specify sea-floor motion during a tsunamigenic earthquake, but can also be used to specify submarine landslide motion or a failing dam, for example.
• A number of new Python modules has been developed to assist the user in working with topo and dtopo files.These are documented in the Clawpack documentation and several of them are illustrated with Jupyter notebooks found in the Clawpack Gallery.
• New capabilities were added in 5.0.0 to monitor the maximum of various flow quantities over a specified time range of a simulation.This capability is crucial for many applications where

PyClaw
PyClaw is an object-oriented Python package that provides a convenient way to set up problems and can be used on large distributed-memory parallel machines.For the latter capability, PyClaw relies on PETSc [2].Lower-level code (whatever gets executed repeatedly and needs to be fast) from the earlier Fortran Classic and SharpClaw codes is automatically wrapped at install time using f2py.
Recent applications of PyClaw include studies of laser light trapping by moving refractive index perturbations [42], instabilities of weakly nonlinear detonation waves [13], and effective dispersion of nonlinear waves via diffraction in periodic materials [28].Two of these are depicted in Figure 6.

Librarization and extensibility
Scientific software is easier to use, extend, and integrate with other tools when it is designed as a library [10].Clawpack has always been designed to be extensible, but PyClaw takes this further in several ways.First, it is distributed via a widely-used package management system, pip.Second, the default installation process ("pip install clawpack") provides the user with a fully-compiled code and does not require setting environment variables.Like other Clawpack packages, PyClaw provides several "hooks" for users to plug in custom routines (for instance, to specify boundary conditions).In PyClaw, these routines -including the Riemann solver itself -are selected at run-time, rather than at compile-time.These routines can be written directly in Python, or (if they are performance-critical) in a

Python geometry
PyClaw includes Python classes for describing collections of structured grids and data on them.These classes are also used by the other codes and VisClaw, for post-processing.A mesh in Clawpack always consists of a set of (possibly mapped) tensor-product grids (interval, quadrilateral, or hexahedral), also referred to as patches.At present, PyClaw solvers operate only on a single patch, but the geometry and grids already incorporate multi-patch capabilities for visualization in AMRClaw and GeoClaw.transverse shocks that arise from instabilities; see [13].Right: Dispersion of waves in a layered medium with matched impedance and periodically-varying sound speed; see [28].

PyClaw solvers
PyClaw includes an interface to both the Classic solvers (already described above) and those of Sharp-Claw [26].SharpClaw uses a traditional method-of-lines approach to achieve high-order resolution in space and time.Spatial operators are discretized first, resulting in a system of ODEs that is then solved using Runge-Kutta or linear multistep methods.The spatial derivatives are computed using a weighted essentially non-oscillatory (WENO) reconstruction from cell averages, which suppresses spurious oscillations near discontinuities.The WENO routines in SharpClaw were generated by PyWENO11 , which is a standalone package that generates WENO routines.
The default time stepping routines in SharpClaw are strong stability preserving (SSP) Runge-Kutta methods of order two to four.Some of the methods use extra stages in order to allow more efficient time stepping with larger CFL numbers.Time stepping in SharpClaw has recently been augmented to include linear multistep methods with variable step size.These methods use a time step size selection that ensures the strong stability preserving property, as described in [22].

Parallelism
PyClaw includes a distributed parallel backend that uses PETSc through the Python wrapper petsc4py.
The parallel code uses the same low-level routines without modification.In the high-level routines, only a few hundred lines of Python code deal explicitly with parallel communication, in order to transfer ghost cell information between subdomains and to find the global maximum CFL number in order to adapt the time step size.For instance, the computation shown in the right part of Figure 6 involved more than 120 million degrees of freedom and was run on two racks of the Shaheen I BlueGene/P supercomputer.The code has been demonstrated to scale with better than 90% efficiency in even larger tests on tens of thousands of processors on both the Shaheen I (BlueGene/P) and Shaheen II (Cray XC40) supercomputers at KAUST.A hybrid MPI/OpenMP version is already available in a development branch and will be included in future releases.

VisClaw : Visualizing Clawpack output
A practical way to visualize the results of simulations is essential to any software package for solving PDEs.This is particularly true for simulations making use of adaptive mesh refinement, since Matlab routines, Iplotclaw allows the user to step, interactively, through a time sequence of plots, jump from one frame to another, or interactively explore data from the current time frame.

Tools for visualizing geo-spatial data produced by GeoClaw
The geo-spatial data generated by GeoClaw has particular visualization requirements.Tsunami or storm surge simulations are most useful when the plots showing inundation or flooding levels are overlaid onto background bathymetry or topography.Supplementary one dimensional time series data (e.g.gauge data) numerically interpolated from the simulation at fixed spatial locations are most useful when compared graphically to observational data.Finally, to more thoroughly analyze the computational data, simulation data should be made available in formats that can be easily exported to GIS tools such as ArcGIS 13 or the open source alternative QGIS14 .For exploration of preliminary results or communicating results to non-experts, Google Earth is also helpful.
The latest release of Clawpack includes many specialized VisClaw routines for handling the above issues with plotting geo-spatial data.Topography or bathymetry data that was used in the simulation will be read by the graphing routines, and, using distinct colormaps, both water and land can be viewed on the same plot.Additionally, gauge locations can be added, along with contours of water and land.
One dimensional gauge plots are also created, according to user-customizable routines.In these gauge plotting routines, users can easily include observational data to compare with GeoClaw simulation results.
In addition to HTML and Latex formats available for all Clawpack results, VisClaw will now also produce KML and KMZ files suitable for visualizing results in Google Earth.Using the same matplotlib graphics routines, VisClaw creates PNG files that can be used as GroundOverlay features in a KML file.Other features, such as gauges, borders on AMR grids, and user specified regions can also be shown on Google Earth.All KML and PNG files are compressed into a single KMZ file that can be opened directly in Google Earth or made available on-line.While VisClaw does not have any direct support for ArcGIS or QGIS, the KML files created for Google Earth can be edited for export, along with associated PNG files to these other GIS applications.

Matlab plotting routines
The Matlab plotting tools available in early versions of Clawpack are still included in VisClaw.
While most of the one and two dimensional capabilities available originally in the Matlab suite have been ported to Python and matplotlib, the original Matlab routines are still available in the Matlab suite of plotting tools.Other plotting capabilities, such as two dimensional manifolds embedded in three dimensional space, or three dimensional plots of fully three-dimensional data are only available in the Matlab routines in a way that interfaces directly with Clawpack.More advanced three-dimensional plotting capabilities are planned for future releases of VisClaw.

Conclusions
Clawpack has evolved over the past 20 years from its genesis as a small and focused software package that two core developers could manage without version control.It is now an ecosystem of related projects that share a core philosophy and some common code (notably Riemann solvers and visualization tools), but that are aimed at different user communities and that are developed by overlapping but somewhat distinct groups of developers scattered at many institutions.The adoption of better software engineering practices, in particular the use of Git and GitHub as an open development platform and the use of pull requests to discuss proposed changes, has been instrumental in facilitating the development of many of the new capabilities summarized in this paper.

Future Plans
The Clawpack development team continues to look forward to new ideas and efforts that will allow great accessibility to the project as well as new capabilities that the core development team has not thought of.To this end a number of the broad efforts that are being considered for the next major release of Clawpack include • An increased librarization effort with the Fortran based sub-packages, • An extensible and more accessible interface to the Riemann solvers, • An effort to allow PyClaw and the Clawpack Fortran packages to rely on more of the same code-base, • An increased emphasis on a larger development community, • More support for new frameworks such as ForestClaw, • A refactoring of the visualization tools in VisClaw, along with support for additional backends, particularly for three-dimensional results (e.g.VisIt15 , ParaView16 , or yt17 ).
4.x to the most recent version 5.3.A number of the repositories have seen only minor changes as the bulk of the development is focused on current research interests.There are a number of minor changes not listed here and the interested reader is encouraged to refer to the change logs7 and the individual

Figure 1 :
Figure 1: AMRClaw example demonstrating a shock-bubble interaction in the Euler equations of compressible gas-dynamics at two times, illustrating the need for adaptive refinement to capture localized behavior.The 40 × 10 × 10 grid at Level 1 is refined where needed by factors of 4 and then 2 in this 3-level run.

Figure 2 :
Figure 2: Left is strong scaling results for the AMRClaw example shown in Fig. 1.Right is plot of efficiency based on total computational time.

Figure 3 :
Figure 3: Left: Gray's Harbor showing Westport, WA on southern peninsula.Middle: Simulation of a potential magnitude 9 Cascadia Subduction Zone event, 40 minutes after the earthquake.Right: Design for new Ocosta Elementary School in Westport, based in part on GeoClaw simulations [19].
and call the algorithms of Clawpack.It grew from what was initially a set of data structures and file IO routines that are used by the other Clawpack codes and by VisClaw.These routines were released in an early form in later 4.x versions of Clawpack.Those releases also included a fullyfunctional implementation of the 1D classic algorithm in pure Python.That implementation still exists in PyClaw and is useful for understanding the algorithm.The current release of PyClaw includes access to the classic algorithms as well as the high-order algorithms introduced in SharpClaw [27] (i.e., WENO reconstruction and Runge-Kutta integrators)

Figure 4 :
Figure 4: Top Left: A snapshot of a GeoClaw storm surge simulation of Hurricane Ike at landfall.Top Right: Tide gauge data computed from GeoClaw and adcirc along with observed data at the same location.Bottom: Computational effort and timings for GeoClaw and adcirc.From [39].

Figure 6 :
Figure 6: Left: A two-dimensional detonation wave solution of the reactive Euler equations, showingtransverse shocks that arise from instabilities; see[13].Right: Dispersion of waves in a layered medium with matched impedance and periodically-varying sound speed; see[28].
most available visualization packages do not have tools that conveniently visualize hierarchical AMR data.VisClaw provides support for all of the main Clawpack submodules, including ClassicClaw, AMRClaw, PyClaw and GeoClaw.From the first release in 1994, Clawpack has included tools for visualizing the output of Clawpack and AMRClaw runs.Up until the release of version Clawpack 4.x, these visualization tools consisted primarily of Matlab routines for creating one, two and three dimensional plots including pseudo-color plots, Schlieren plots, contour plots and scatter plots, including radially or spherically symmetric data.Built-in tools were also available for handling one, two and three-dimensional mapped grids.Starting with version 4.x, however, it was recognized that a reliance on proprietary software for visualization prevented a sizable potential user base from making use of the Clawpack software.The one and two dimensional plotting routines were converted from Matlab to matplotlib, a popular open source Python package for producing publication quality graphics for one and two dimensional data[23].With the development of Clawpack Version 5 and above, Python graphics tools have been collected into the VisClaw repository.The VisClaw tools extend the functionality of the Version 4.x Python routines for creating one and two dimensional plots, and adds several new capabilities.Chief among these are the ability to generate output to webpages, where a series of plots can be viewed individually or as an animated time sequence using the Javascript package12 (which was motivated by code in an earlier version of Clawpack).The VisClaw module Iplotclaw provides interactive plotting capabilities from the Python or IPython prompt.Providing much of the same interactive capabilities as the original : A Community-Driven Collection of Approximate Riemann Solvers The methods implemented in Clawpack, and all modern Godunov-type methods for hyperbolic PDEs, are based on the solution of Riemann problems as discussed in Section 1.2.Whereas most existing codes for hyperbolic PDEs use Riemann solvers to compute fluxes, Clawpack Riemann solvers instead 7http://www.clawpack.org/changes.html 8http://docs.scipy.org/doc/numpy-dev/f2py