Algorithmic Approach
Our Algorithmic Approach is based on finite difference discretizations of PDEs on locally
structured grids, combined with techniques for introducing adaptivity
and irregular geometry, plus methods for representing kinetic problems
using particles.
Our approach corresponds to a line of investigation that has been the
subject of research on the part of the PIs in this proposal for over a
decade. It has been highly successful in a variety of applications areas, and
represents a good match to the requirements of the applications described
above.
High-resolution finite difference methods on structured grids
In grid-based methods for numerical PDEs, it is necessary to make a
fundamental choice of discretization technologies. We have chosen to use
locally structured grids, i.e. ones that are based on defining discrete
unknowns on a rectangular
discretization of the spatial independent variables.
This is in contrast to
unstructured grids, in which the discretized
dependent variables are defined at
nodes in space corresponding to a graph with more general
connectivity.
We will focus on
finite difference methods on such grids, although much of the
lower-level infrastructure could be used equally well for finite element
discretizations.
For the time-dependent problems described here, we believe
that there are significant advantages to the structured grid approach
(although we also note that unstructured finite element methods
remain the methods of choice for problems such as computing
eigenmodes in complex geometries). In the
applications described here, structured grid methods have been used
extensively, so there is a body of experience in those communities with
numerical methods based on that approach.
It is well-understood how to construct stable and arbitrarily high order
accurate discretizations to PDEs on such grids.
Regularity in space leads to regular data layouts, making it easier
to optimize for various types of data locality.
The regular geometric structure of the spatial distribution of unknowns
leads to fast efficient iterative solvers for elliptic and
parabolic problems based on multigrid iteration.
Adaptive gridding based on local refinement of structured grids,
discussed below, is
actually a more mature technology for time-dependent problems in which
one must regrid often: comparably efficient grid generation for
unstructured grid methods is still an open research question.
Coupling to particle methods is inexpensive and well-understood.
Locating the particles on the grid is trivial,
and interpolation between particles and
fields has an extensive
body of practical experience and mathematical analysis.
Volume-of-Fluid Methods for Irregular Boundaries
There are two closely-related requirements for representing irregular
boundaries. To represent engineering geometries,
we need to compute solutions in domains whose
whose shape (and trajectory, if they are moving) are known. In
addition, in the gas jet problem for plasma-laser accelerators, and
in detailed modeling of the initial breakup of a fuel jet
for combustion, we need to be able to represent free boundaries,
for which PDE's are being solved on both sides of a boundary
whose location and shape is changing as a function of time and
is computed as part of of the solution. In both cases, we will use a
volume-of-fluid representation of the discretized solution near the
boundary. In this approach, the
surface is represented by its intersection with an underlying
rectangular grid. This leads to a natural, finite-volume discretization
of the PDE on irregular control volumes adjacent to the boundary.
In the case where the surface is a domain boundary, this approach is
variously referred to as a Cartesian grid, embedded boundary, or cut
cell method (in this proposal, we will use the term embedded
boundary). One of the principal advantages of the embedded boundary
method is that the problem of generating the description of the geometry
on the grid starting from surface tesselations produced, for example,
by a CAD system
has been completely solved. In
addition, there has been considerable progress on the development of
discretizations for this type of geometry representation. In this approach,
the primary dependent variables are centered at the centers of the
rectangular grid cells, even if those cells correspond to irregular
control volumes. Because the values are centered at regular points on
the grid, it is routine to construct finite difference
approximations to various partial differential operators on the
irregular control volumes that are consistent, conservative,
and stable/well-conditioned.
In the case in which the PDE is being solved on both
sides of a free boundary, the method is referred to as a volume-of-fluid
method.
We use a similar approach for the case of free boundaries.
Solutions are double-valued in each cell,
with both values centered at the Cartesian cell centers.
The construction of difference approximations to various operators
follows along the lines of that in the embedded boundary case, except
that the appropriate well-posed boundary conditions are applied
along the free boundary.
In using these ideas to represent two-phase flow at low Mach numbers, we
have a variety of options, from tracking only the density jump and representing all of the other dependent
variables as single-valued, to making all quantities double-valued, as was done for tracking hyperbolic fronts.
Adaptive Mesh Refinement
Our adaptive mesh methods are based on the block-structured adaptive
mesh refinement (AMR) algorithms of Berger and Oliger.
In this approach, the regions to be refined
are organized into rectangular patches of several hundred to several
thousand grid cells per patch. Thus one is able to use the
high-resolution rectangular grid methods described above to advance the
solution in time. Furthermore, the overhead in managing the irregular
data is amortized over a relatively large amounts or floating point work
on regular arrays. For time-dependent problems, refinement is performed
in time as well as in space. Each level of spatial refinement uses
its own stable time step, with the time steps at a level
constrained to be integer multiples of the time steps at all finer
levels. AMR is a mature technology for problems without geometry,
with a variety of implementations
and applications for various nonlinear combinations of elliptic, parabolic
and hyperbolic PDE's. In particular, the collection of
applications here address most of the major algorithmic issues in
developing adaptive algorithms for the applications described above,
such as
adaptive multigrid solvers for Poisson's equation, the representation of
non-ideal effects in MHD, and the coupling of particle methods to AMR field
solvers. AMR has also been successfully coupled to the
volume-of-fluid algorithms described above.
Particle Methods for Kinetic Problems
To solve the kinetic equations as they arise in
plasma physics, particle-beam modeling, and spray modeling,
particle-in-cell methods are most often used.
In this approach, the distribution
function of the particles is represented by a collection of
macro-particles, while the spatial field variables are represented as
functions on a discrete grid in physical space, with the equations for
those fields discretized using finite differences or finite elements.
The coupling between the particles and the fields is implemented using
interpolation. Quantities such as charges and currents that are
required as sources in the field
equations are interpolated onto the grid from the particle
representation, and field variables required for the particle
dynamics are interpolated to the particle locations from the grid.
Within this approach, there are a range of design variations, based on
tradeoffs between accuracy and computational cost. In an electrostatic
PIC code, one can interpolate charge onto the grid,
solve Poisson's equation for the electric field potential,
and interpolate a difference
approximation to the electric field onto the particles.
One could also interpolate
the dipole field onto the grid, solve three Poisson's equations
for the electric field components, and interpolate the field directly
from the grid. The latter is more accurate, at the cost of increasing the
number of Poisson solutions.
In addition to pursuing the classical PIC approaches to particle methods for
kinetics problems, we
will also provide support for some of the methods that have been
developed in the mathematics community.
There are a variety of fast and accurate algorithms for evaluating
electrostatic potentials
that may be of use in this setting, such as
the fast multipole method (FMM) and the method of
local corrections (MLC). For
both of these methods, the nonlocal component of the
potential induced by a collection of
N particles can be computed in no more than O(N log N)
operations in a way that
essentially eliminates the self-force problem and decouples the
effect of the short-range smoothing of the potential from the
grid spacing (in fact, FMM is a grid-free method).
In addition, the MLC has been applied to the problem of
interpolating forces onto the grid for general particle / mesh
methods, and in particular can be applied
to the case of Vlasov-Maxwell, as well as to the particle methods
arising in the simulation of sprays in combustion.