Parallel operators with PETSc (FCI)#

BOUT++ provides a PETSc-based interface for constructing parallel operators on FCI grids from mesh metadata. These operators are intended for cases where the standard field-based parallel operators described in Differential operators are not the best fit, for example when using FCI and wanting a conservative support-operator discretisation built from the traced field-line geometry.

This interface constructs sparse matrix PETSc operators that act on Field3D data and the associated FCI parallel-boundary values. The main entry point is the PetscOperators class in bout/petsc_operators.hxx.

When to use these operators#

Use the PETSc FCI operators when all of the following are true:

  • The mesh uses the FCI parallel transform, mesh:paralleltransform:type = fci

  • BOUT++ is built with PETSc support

  • BOUT++ is built with 3D metrics enabled

  • The grid file includes the FCI operator weights generated during grid construction.

If you are using the identity or shifted-metric parallel transforms, or if you simply need the standard field operators such as Grad_par or Div_par_K_Grad_par, use the routines documented in Differential operators instead.

Requirements#

Build requirements#

This feature requires:

  • PETSc support in BOUT++

  • 3D metrics enabled in BOUT++

See PETSc for PETSc installation and configuration details.

Mesh and grid requirements#

These operators are designed for the FCI transform described in FCI method. The grid file must contain the additional numbering, interpolation, and weight arrays needed to construct the PETSc matrices. These are typically generated by the Zoidberg-based FCI grid tooling.

To construct PetscOperators, the mesh metadata currently required are:

  • Cell-space numbering fields:

    • cell_number

    • forward_cell_number

    • backward_cell_number

  • Cell-space size information:

    • total_cells

  • Leg-numbering fields:

    • forward_leg_interior_number

    • forward_leg_boundary_number

    • backward_leg_interior_number

    • backward_leg_boundary_number

  • Leg-space size information:

    • n_forward_legs

    • n_backward_legs

  • CSR data for the forward and backward interpolation operators:

    • forward_rows

    • forward_columns

    • forward_weights

    • backward_rows

    • backward_columns

    • backward_weights

  • Leg-weight arrays:

    • forward_leg_weights

    • backward_leg_weights

If these metadata arrays are not present, constructing PetscOperators will fail during setup.

Boundary conditions and communication#

The PETSc FCI operators use the same parallel boundary data as the rest of the FCI workflow. In particular, parallel boundary conditions must be applied to the parallel slices after communication, as described in FCI boundary conditions.

A typical sequence is:

f.applyBoundary();
mesh->communicate(f);
f.applyParallelBoundary("parallel_neumann_o1");

The exact boundary operator should match the operator you are constructing or comparing against. For example, diffusion operators usually use Neumann parallel-boundary values, while other comparisons may use Dirichlet values.

Constructing the operators#

Include the header:

#include <bout/petsc_operators.hxx>

Construct the operator factory from the mesh, then build the grouped parallel operators:

PetscOperators ops(mesh);
auto parallel = ops.getParallel();

The returned parallel object contains:

  • Cell-to-cell operators such as Grad_par, Div_par, and Div_par_Grad_par

  • Half-step operators on the forward and backward FCI legs

  • Interpolation, injection, and restriction operators between cell space and leg space

  • The cell volume field dV used in the support-operator construction

These operators can then be applied directly to Field3D values.

The PetscOperator class#

The sparse operator type is a templated PetscOperator class. Operators map between a row space and a column space, which are tracked in the type system. In practice, the most common custom operator is likely to be a PetscCellOperator, which maps from cell space to cell space.

Available aliases include:

  • PetscCellOperator for cell-to-cell operators

  • PetscForwardOperator and PetscBackwardOperator for cell-to-leg operators

  • PetscForwardToCellOperator and PetscBackwardToCellOperator for leg-to-cell operators

Creating a PetscOperator#

There are two main ways to create a PetscOperator:

  • From CSR arrays stored in the mesh file

  • From an already-assembled PETSc Mat

The CSR constructor takes:

  • An output-space mapping

  • An input-space mapping

  • A CSR row-pointer array rows

  • A CSR column-index array cols

  • A CSR weight array weights

This is how the built-in forward and backward FCI interpolation operators are constructed:

PetscForwardOperator forward(
    forward_leg_mapping, cell_mapping,
    meshGetArray<int>("forward_rows"),
    meshGetArray<int>("forward_columns"),
    meshGetArray<BoutReal>("forward_weights"));

PetscBackwardOperator backward(
    backward_leg_mapping, cell_mapping,
    meshGetArray<int>("backward_rows"),
    meshGetArray<int>("backward_columns"),
    meshGetArray<BoutReal>("backward_weights"));

To implement a new operator, one straightforward route is:

  1. Generate the CSR arrays for the operator you want

  2. Write those arrays into the mesh or grid metadata

  3. Read them back with meshGetArray<int>() / meshGetArray<BoutReal>()

  4. Construct the corresponding PetscOperator

For many custom differential operators, PetscCellOperator is the natural starting point because both the input and output live in the cell space.

Operator application and algebra#

Once constructed, PetscOperator supports:

  • Application to compatible PETSc vectors

  • Application to Field3D when the input space is cell space

  • transpose() to swap input and output spaces

  • Operator addition: A + B

  • Operator subtraction: A - B

  • Scalar multiplication: A * alpha

  • Operator composition: A * B when the inner and outer spaces match

  • Diagonal operator construction with PetscOperator::diagonal(...) for same-space operators

For example:

auto Grad_plus = Inv_dl_plus * (Forward - Inject_plus);
auto Div_minus = Neg_inv_dV * Grad_plus.transpose() * DV_leg_plus;
auto Grad_par = ((Restrict_minus * Grad_plus)
                 + (Restrict_plus * Grad_minus)) * 0.5;

This operator algebra is one of the main advantages of the PETSc-based representation: once the primitive interpolation and injection operators exist, new support-operator constructions can be built by combining them directly.

Example#

The basic usage pattern is:

PetscOperators ops(mesh);
auto parallel = ops.getParallel();

Field3D f = ...;
mesh->communicate(f);
f.applyParallelBoundary("parallel_neumann_o1");

Field3D grad_f = parallel.Grad_par(f);
Field3D div_grad_f = parallel.Div_par_Grad_par(f);

Field3D K = 1.0;
mesh->communicate(K);
K.applyParallelBoundary("parallel_neumann_o1");

Field3D diffusive = parallel.Div_par_K_Grad_par(K, f);

The integrated test in tests/integrated/test-petsc-operators provides a more complete example comparing these PETSc operators with the existing field-based parallel operators on an FCI grid.

Relation to the standard parallel operators#

The standard routines in Differential operators operate directly on fields and are the default choice in most models. The PETSc FCI operators documented here are an additional interface for constructing sparse matrix operators from FCI geometry. The main advantage is that sparse matrices can be manipulated to e.g. calculate divergence operators from gradient operators (the Support Operator Method).

In practice:

  • Use Differential operators for the usual BOUT++ differential-operator workflow

  • Use this page when you need the PETSc/FCI operator construction explicitly

Performance and limitations#

Operator setup is substantially more expensive than applying an already-built operator, because the PETSc matrices must be assembled from the FCI metadata. Construct the operators once, then reuse them.

These operators also require additional memory for the PETSc matrices and vectors, so they are best suited to cases where that setup cost is justified by repeated operator application or by the need for the support-operator construction.

Current limitations include:

  • The interface is intended for FCI grids, not the identity or shifted-metric transforms

  • PETSc and 3D metrics are required

  • The grid metadata must be generated ahead of time

See also#