.. _sec-parallel-operators-petsc-fci: 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 :ref:`sec-diffops` 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 :ref:`sec-diffops` instead. Requirements ------------ Build requirements ~~~~~~~~~~~~~~~~~~ This feature requires: - PETSc support in BOUT++ - 3D metrics enabled in BOUT++ See :ref:`sec-PETSc-install` for PETSc installation and configuration details. Mesh and grid requirements ~~~~~~~~~~~~~~~~~~~~~~~~~~ These operators are designed for the FCI transform described in :ref:`sec-fci`. 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 :ref:`sec-parallel-bc-fci`. A typical sequence is: .. code-block:: C++ 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: .. code-block:: C++ #include Construct the operator factory from the mesh, then build the grouped parallel operators: .. code-block:: C++ 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: .. code-block:: C++ PetscForwardOperator forward( forward_leg_mapping, cell_mapping, meshGetArray("forward_rows"), meshGetArray("forward_columns"), meshGetArray("forward_weights")); PetscBackwardOperator backward( backward_leg_mapping, cell_mapping, meshGetArray("backward_rows"), meshGetArray("backward_columns"), meshGetArray("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()`` / ``meshGetArray()`` 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: .. code-block:: C++ 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: .. code-block:: C++ 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 :ref:`sec-diffops` 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 :ref:`sec-diffops` 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 -------- - :ref:`sec-diffops` - :ref:`sec-fci` - :ref:`sec-parallel-bc-fci` - :ref:`sec-petsc`