This library provides a C++ namespace fd
containing a class representing a finite-difference differential operator on a rectangular grid backed by PETSc. To use the class, called fd::RectangularOperator
, simply inherit from it and provide an overload of the Assemble()
function, which assembles the
operator. We provide an example of this with the fd::Laplacian
class, which implements a second-order centered finite difference Laplacian.
You may also want to overload any of the virtual functions in the fd::RectangularOperator
class, particularly ApplyBC
and MatApplyBC
, which apply
boundary conditions to a PETSc Vec
and a Mat
respectively. The default is to apply zero Dirichlet boundary conditions.
fd::RectangularOperator
is a functor in C++11-speak, meaning that it is a class that implements operator()
-- specifically for
fd::RectangularOperator
, we have
virtual PetscErrorCode fd::RectangularOperator::operator()(Vec X, Vec Y) const noexcept;
which applies the operator to a PETSc Vec
X
, storing the result in Y
. Note that PETSc allows you to construct a matrix-free linear operator
in the form of a MatShell
, which you can create in the Assemble()
function if you want to apply your operator matrix-free but still have the features
of a PETSc
Mat
and of the fd::RectangularOperator
class. That way, when we apply the matrix in operator()
, everything works as expected. You could
also overload operator()
itself, but you probably only want to do this if you cannot for whatever reason make a MatShell
to do what you want
(one possible reason could be that your implementation throws C++ exceptions that you don't necessarily want to trap inside your callback and convert
to an appropriate PetscErrorCode
return value, which must be done to guarantee that your exception isn't simply ignored at the C call boundary or worse).
The library also contains utilities for performing tests with the method of manufactured solutions (see fd::TestConvergence<>
in include/Utils.h
),
which is useful to check that you've implemented the discretization correctly in Assemble()
and applied the appropriate boundary conditions. We
provide integration with GNUPlot with the fd::VecToGNUPlot()
function (in include/Utils.h
), and a utility function to
scatter a PETSc Vec
to MPI rank 0 in only 1 function call instead of 4 (useful for plotting a solution that was obtained via multiple-processor solve).
First, download and install PETSc from https://www.mcs.anl.gov/petsc/download/index.html
if you haven't already. See https://www.mcs.anl.gov/petsc/documentation/installation.html
for PETSc installation instructions.
Finally, cd
into petsc_fd
and
make PETSC_DIR=/your/petsc-dir PETSC_ARCH=your-petsc-arch