Petsc matrix operators for FCI by bendudson · Pull Request #3330 · boutproject/BOUT-dev

bendudson · 2026-03-16T20:56:48Z

Reads matrix operators from the mesh file in CSR format, creating PETSc matrices that apply to Field3Ds. Mainly intended to implement FCI parallel operators as an alternative to yup/down communication.

This enables:

Arbitrary parallelisation : PETSc handles the communication pattern, so there is no restriction on the field-line trajectories and processor boundaries. The mapping between mesh global indices and petsc global indices is handled automatically.
Exact conservation of fluxes: Using the Support Operator Method to derive the divergence operator from the transpose of the gradient operator (as in e.g. GRILLIX) ensures that the numerical integral of the divergence of a flux cancels to numerical precision.
Better preconditioning: The LU/MUMPS/STRUMPACK preconditioning used in the snes solver needs to know which cells are connected, so that PETSc can perform coloring and calculate the system Jacobian efficiently. Knowing the forward and backward matrices will enable this to be done.

There are a couple of downsides:

Linearity: The operators are fixed matrices so are linear. Monotonic splines, for example, require nonlinearity.
More communications: Each matrix-vector operation will perform communications. Perhaps there is a way to batch operations into some kind of composite matrix?

Read operators from the mesh in CSR format. PETSc handles communication pattern for MPI decomposition.

Handles the mapping between mesh indices and PETSc global indices.

Simple wrapper functions around the Options functionality, so that Array types can be read from the mesh file.

Reads row, column and weights from Mesh, constructs Petsc Mat objects. Working on composition and operation on Field3D objects.

Map from Field3D, including yup/ydown boundaries, to PETSc Vec, perform matrix-vector product, and map back to a Field3D.

Docstrings to explain PetscMapping, PetscOperator and PetscOperators classes. Extend `BOUT_DO_PETSC` macro to print error message, file and line number.

Alternative to initializer_list constructor that avoids confusion with e.g. `Array{2}` : Is it an array of size 2, or an array of size 1 containing value 2? `std::vector` uses `vector(2)` means size 2; `vector{2}` means value 2. With this commit, both `Array{2}` and `Array(2)` mean size 2; `Array<T>::fromValues({2})` means size 1, value 2. It's more verbose but less prone to typos.

Some care needed with yup/down and yp/ym indices.

`diagonal` creates diagonal matrices from Field3D. `transpose` creates a new matrix that is a transpose of a PetscOperator. This is useful for converting gradients into divergences.

Allows input operators to contain entries in boundary rows. This is important so that the gradient and divergence operators depend on boundary cell values.

github-actions

clang-tidy made some suggestions

github-actions · 2026-03-16T21:01:59Z

+  /// The local index is the offset into the packed local PETSc vector.
+  template <typename Function>
+  void map_local_yup(Function func) const {
+    PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size();


warning: narrowing conversion from 'unsigned int' to signed type 'PetscInt' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]

PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); ^

github-actions · 2026-03-16T21:01:59Z

+  /// The local index is the offset into the packed local PETSc vector.
+  template <typename Function>
+  void map_local_ydown(Function func) const {
+    PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size()


warning: narrowing conversion from 'unsigned int' to signed type 'PetscInt' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]

PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() ^

github-actions · 2026-03-16T21:01:59Z

+/// - applied to a `Field3D`;
+/// - composed with other operators; and
+/// - added or subtracted.
+class PetscOperator {


warning: class 'PetscOperator' defines a destructor but does not define a copy constructor, a copy assignment operator, a move constructor or a move assignment operator [cppcoreguidelines-special-member-functions]

class PetscOperator { ^

github-actions · 2026-03-16T21:02:01Z

+  int mesh_total_cells;
+  if (mesh->get(mesh_total_cells, "total_cells") == 0) {
+    // Check total number of cells
+    if (mesh_total_cells != mapping->size()) {


warning: comparison of integers of different signs: 'int' and 'unsigned int' [clang-diagnostic-sign-compare]

if (mesh_total_cells != mapping->size()) { ^

github-actions · 2026-03-16T21:02:01Z

  IMEXBDF2* s;
-  ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s));
-  CHKERRQ(ierr);
+  PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));


warning: do not use reinterpret_cast [cppcoreguidelines-pro-type-reinterpret-cast]

PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s))); ^

github-actions · 2026-03-16T21:02:01Z

  IMEXBDF2* s;
-  ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s));
-  CHKERRQ(ierr);
+  PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));


warning: multilevel pointer conversion from 'void **' to 'void *', please use explicit cast [bugprone-multi-level-implicit-pointer-conversion]

PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s))); ^

Not available unless configured with PETSc.

github-actions

clang-tidy made some suggestions

github-actions · 2026-03-16T21:59:05Z

+  return Region<Ind3D>(xout_indices);
+}
+
+PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number,


warning: constructor does not initialize these fields: row_start, row_end, mat_mesh_to_petsc, mat_petsc_to_mesh [cppcoreguidelines-pro-type-member-init]

include/bout/petsc_operators.hxx:236:

- PetscInt row_start, row_end; + PetscInt row_start{}, row_end{};

include/bout/petsc_operators.hxx:239:

- Mat mat_mesh_to_petsc; + Mat mat_mesh_to_petsc{};

include/bout/petsc_operators.hxx:242:

- Mat mat_petsc_to_mesh; + Mat mat_petsc_to_mesh{};

github-actions · 2026-03-16T21:59:05Z

+  MatDestroy(&this->mat_petsc_to_mesh);
+}
+
+PetscOperator::PetscOperator(PetscMappingPtr mapping, Array<int> rows, Array<int> cols,


warning: constructor does not initialize these fields: mat_operator [cppcoreguidelines-pro-type-member-init]

include/bout/petsc_operators.hxx:335:

- Mat mat_operator; + Mat mat_operator{};

github-actions · 2026-03-16T21:59:05Z

+  BOUT_DO_PETSC(MatDestroy(&mat));
+}
+
+void PetscOperator::copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec) {


warning: the parameter 'mapping' is copied for each invocation but only used as a const reference; consider making it a const reference [performance-unnecessary-value-param]

Suggested change

void PetscOperator::copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec) {

void PetscOperator::copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec) {

include/bout/petsc_operators.hxx:345:

- static void copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec); + static void copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec);

github-actions · 2026-03-16T21:59:06Z

+  return PetscOperator(this->mapping, mat);
+}
+
+Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const {


warning: method 'meshGetField3D' can be made static [readability-convert-member-functions-to-static]

Suggested change

Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const {

Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) {

include/bout/petsc_operators.hxx:392:

- Field3D meshGetField3D(Mesh* mesh, const std::string& name) const; + static Field3D meshGetField3D(Mesh* mesh, const std::string& name) ;

Check number of cells against global mapping size, not local size.

Handle non-standard CSR inputs: - Negative indices for missing rows - Missing final entry in rows array

github-actions

clang-tidy made some suggestions

Changes PetscMapping::size() to PetscMapping::localSize()

Defined in `bout/petsc_interface.hxx`, `bout::petsc::UniqueVec and `UniqueMat` are `std::unique_ptr` wrappers around PETSc Vec and Mat pointers, with custom deleters. PetscOperator copy constructor and assignment deleted, move assignment set to default. This enables some optimizations to operate in-place on temporaries.

ZedThree · 2026-03-18T11:49:42Z

This looks amazing, thanks @bendudson!

I've not finished going through this yet, but how much is this a replacement for GlobalIndexer and/or PetscMatrix/PetscVector? Or does it work together with them? It looks like this can only do Field3D, whereas GlobalIndexer can do the other fields too?

github-actions

clang-tidy made some suggestions

github-actions · 2026-03-18T14:45:09Z

+struct VecDeleter {
+  void operator()(Vec* vec_ptr) const {
+    VecDestroy(vec_ptr);
+    delete vec_ptr;


warning: deleting a pointer through a type that is not marked 'gsl::owner<>'; consider using a smart pointer instead [cppcoreguidelines-owning-memory]

delete vec_ptr; ^

Additional context

include/bout/petsc_interface.hxx:67: variable declared here

void operator()(Vec* vec_ptr) const { ^

github-actions · 2026-03-18T14:45:09Z

+struct MatDeleter {
+  void operator()(Mat* mat_ptr) const {
+    MatDestroy(mat_ptr);
+    delete mat_ptr;


warning: deleting a pointer through a type that is not marked 'gsl::owner<>'; consider using a smart pointer instead [cppcoreguidelines-owning-memory]

delete mat_ptr; ^

Additional context

include/bout/petsc_interface.hxx:80: variable declared here

void operator()(Mat* mat_ptr) const { ^

Shows how to use the PetscOperators, and compares against the yup/ydown operators. Includes a script to generate the mesh, which currently requires a branch of Zoidberg to run. The code to calculate Grad_par, Div_par and Div_par_Grad_par from forward and backward operators is now in `PetscOperators::getParallel()`. That calculates and returns a collection of `Parallel` operators. This example could be turned into a test: - `Grad_par` should be close between the two methods - Volume integrals of `Div_par` and `Div_par_Grad_par` results should be near zero (flux conserved).

bendudson · 2026-03-18T16:45:53Z

Thanks for going through this @ZedThree ! Sorry there's quite a lot, though most is self-contained. The implementation is independent of GlobalIndexer, PetscVector and PetscMatrix. Perhaps there is a way to reuse that machinery, but the way these operators are created is quite specific: The mesh generator assigns global indices to cells, including y boundary cells, and PetscMapping is calculating the mapping between the mesh global indices and PETSc's global indices that depend on the number of processors / domain decomposition.

I've added an example of how these operators are used (https://github.com/boutproject/BOUT-dev/pull/3330/changes#diff-d24773f587abce8257e893cfcfec54002aafaf875fb9cf93f54076b5c60ec7abR16). There is some amount of setup needed to calculate differential operators from the forward / backward maps, that I've put into a helper function that returns a struct of operators (https://github.com/boutproject/BOUT-dev/pull/3330/changes#diff-62e947100b5e999ca6d6d57a6f1281449adc7e2db9e8b7fc5e9273fea887ad88R302). Not sure this is the right API; suggestions welcome!

ZedThree · 2026-03-18T16:51:26Z

The mesh generator assigns global indices to cells, including y boundary cells, and PetscMapping is calculating the mapping between the mesh global indices and PETSc's global indices that depend on the number of processors / domain decomposition.

I had thought that's essentially what GlobalIndexer was doing, so possibly I am misunderstanding either that or this! Maybe we can discuss at the next dev meeting?

bendudson · 2026-03-18T17:05:43Z

The mesh generator assigns global indices to cells, including y boundary cells, and PetscMapping is calculating the mapping between the mesh global indices and PETSc's global indices that depend on the number of processors / domain decomposition.

I had thought that's essentially what GlobalIndexer was doing, so possibly I am misunderstanding either that or this! Maybe we can discuss at the next dev meeting?

The GlobalIndexer assigns a global index to each cell but doesn't map this to a mesh global index that is independent of the domain decompositions: Changing the number of processors will change the global indices. The key part of PetscMapping is the consistency of cell numbering in the mesh file, row and column indices in the CSR matrices for the operators.

ZedThree · 2026-03-18T17:13:27Z

Ahhh, I see! Thanks!

This is a function rather than a PetscOperator because it requires at least two matrix-vector products. This implementation requires six Mat-Vec products that could probably be optimized.

Trace multiple field lines per cell, to construct smoother maps. A complication is cells that partially intersect a boundary.

Mesh not available in repo, only generation script.

…etsc-operators

github-actions

clang-tidy made some suggestions

dschwoerer · 2026-06-05T09:30:05Z

I noticed we run 56 test in the "static, optimized case", including the test-petsc-operator, while fedora runs 52 + 8 + 1 - and the test-petsc-operator is not among them. It does get however build, just not run.

The Fedora CI run has 3 additional tests disabled, according to the configure stage:

-- Not building test test-laplace-hypre3d, requirement not met: BOUT_HAS_HYPRE
-- Not building test test-laplacexy2-hypre, requirement not met: BOUT_HAS_HYPRE
-- Not building test test-slepc-solver, requirement not met: BOUT_HAS_SLEPC

So on the non-fedora CI should have 64 tests run - so 8 tests are missing.
Taking into account that test-petsc-operator is missing on fedora, that makes at least 9 missing tests.

Any idea why some tests are not run @ZedThree ?

…etsc-operators

github-actions

clang-tidy made some suggestions

github-actions · 2026-06-17T09:40:37Z

  Ind3D _ind() const { return region->bndry_points[pos].index; }
  signed char _boundary_width() const {
-    return region->localmesh->ystart - region->bndry_points[pos].abs_offset;
+    return region->localmesh->ystart - region->bndry_points[pos].abs_offset + 1;


warning: narrowing conversion from 'int' to signed type 'signed char' is implementation-defined [bugprone-narrowing-conversions]

{ ^

github-actions

clang-tidy made some suggestions

github-actions

clang-tidy made some suggestions

ZedThree · 2026-06-19T12:37:42Z

Any idea why some tests are not run @ZedThree ?

Looks like the Fedora build runs make -C build check, but the other builds run cd build && ctest.

Here's the diff between what's run on the latest CI:

diff /tmp/static_optimised /tmp/fedora 
1,2d0
< MMS-bracket
< MMS-derivatives3
4,5d1
< MMS-shiftedmetricinterp
< MMS-upwinding3
11,16d6
< test-boutpp-collect
< test-boutpp-collect-staggered
< test-boutpp-legacy-model
< test-boutpp-mms-ddz
< test-boutpp-print
< test-boutpp-simple-model
37a28
> test-laplace-hypre3d
40a32
> test-laplacexy2-hypre
48a41
> test-slepc-solver
60a54
> test_beuler

So we've got the hypre and slepc tests which aren't running due to missing the dependencies, but also beuler, which only requires petsc, so not sure why that's skipped.

I think the others all require Python, which isn't found on the the static optimised build.

github-actions

clang-tidy made some suggestions

github-actions · 2026-06-24T23:47:34Z

+  auto* coords = mesh->getCoordinates();
+
+  // Parallel spacing in cell space
+  Field3D dl = coords->dy * sqrt(coords->g_22);


warning: no viable conversion from 'std::enable_if_t<is_expr_field2d_v && is_expr_field2d_v<BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>>, BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>') to 'Field3D' [clang-diagnostic-error]

Field3D dl = coords->dy * sqrt(coords->g_22); ^

Additional context

include/bout/field3d.hxx:181: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v && is_expr_field2d_v<BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>>, BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>') to 'Mesh *' for 1st argument

Field3D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, ^

include/bout/field3d.hxx:189: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v && is_expr_field2d_v<BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>>, BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>') to 'const Field3D &' for 1st argument

Field3D(const Field3D& f); ^

include/bout/field3d.hxx:192: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v && is_expr_field2d_v<BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>>, BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>') to 'Field3D &&' for 1st argument

Field3D(Field3D&& f) noexcept { swap(*this, f); } ^

include/bout/field3d.hxx:194: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v && is_expr_field2d_v<BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>>, BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>') to 'const Field2D &' for 1st argument

Field3D(const Field2D& f); ^

include/bout/field3d.hxx:196: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v && is_expr_field2d_v<BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>>, BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::sqrt>, bout::op::Mul>') to 'BoutReal' (aka 'double') for 1st argument

Field3D(BoutReal val, Mesh* localmesh = nullptr); ^

include/bout/field3d.hxx:203: candidate template ignored: could not match 'Field3D' against 'Field2D'

Field3D(const BinaryExpr<Field3D, L, R, Func>& expr) ^

include/bout/fieldops.hxx:401: candidate function

} ^

include/bout/fieldops.hxx:402: candidate function

(), f}; } ^

github-actions · 2026-06-24T23:47:34Z

+  dl.applyParallelBoundary("parallel_neumann_o1");
+
+  // Cell volume
+  Field3D dV = coords->J * coords->dx * coords->dy * coords->dz;


warning: no viable conversion from 'std::enable_if_t<is_expr_field2d_v<BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>> && is_expr_field2d_v, BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>') to 'Field3D' [clang-diagnostic-error]

Field3D dV = coords->J * coords->dx * coords->dy * coords->dz; ^

Additional context

include/bout/field3d.hxx:181: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v<BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>> && is_expr_field2d_v, BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>') to 'Mesh *' for 1st argument

Field3D(Mesh* localmesh = nullptr, CELL_LOC location_in = CELL_CENTRE, ^

include/bout/field3d.hxx:189: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v<BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>> && is_expr_field2d_v, BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>') to 'const Field3D &' for 1st argument

Field3D(const Field3D& f); ^

include/bout/field3d.hxx:192: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v<BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>> && is_expr_field2d_v, BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>') to 'Field3D &&' for 1st argument

Field3D(Field3D&& f) noexcept { swap(*this, f); } ^

include/bout/field3d.hxx:194: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v<BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>> && is_expr_field2d_v, BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>') to 'const Field2D &' for 1st argument

Field3D(const Field2D& f); ^

include/bout/field3d.hxx:196: candidate constructor not viable: no known conversion from 'std::enable_if_t<is_expr_field2d_v<BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>> && is_expr_field2d_v, BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>>' (aka 'BinaryExpr<Field2D, BinaryExpr<Field2D, BinaryExpr<Field2D, Field2D, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>, Field2D, bout::op::Mul>') to 'BoutReal' (aka 'double') for 1st argument

Field3D(BoutReal val, Mesh* localmesh = nullptr); ^

include/bout/field3d.hxx:203: candidate template ignored: could not match 'Field3D' against 'Field2D'

Field3D(const BinaryExpr<Field3D, L, R, Func>& expr) ^

include/bout/fieldops.hxx:401: candidate function

} ^

include/bout/fieldops.hxx:402: candidate function

(), f}; } ^

- Arithmetic expression on FieldMetrics - Check for null source in Mesh::get - Use BOUT_USE_METRIC_3D rather than BOUT_HAS_METRIC_3D

Iterate over elements to preallocate the PETSc matrix before setting element values.

Changes in `next` mean communicating Bxy throws exception (disabled calc parallel slices).

dschwoerer

Looks good, thanks @bendudson!

dschwoerer · 2026-06-30T09:23:25Z

+  /// Convert (x, y, z) to an Ind3D
+  /// Requires mesh size.
+  Ind3D indexAt(int x, int y, int z) const {
+    const int ny = this->getNy();
+    const int nz = this->getNz();
+    return Ind3D{(((x * ny) + y) * nz) + z, ny, nz};
+  }
+


I have also added this in the past in many places. Is the Field the best place for this? Or should it go to to the mesh?

Would it be worth to have it for all fields?

Perhaps in Field, or adding to other fields e.g. Field2D::indexAt should return an Ind2D. It's definitely a useful utility; I think Hermes-3 has at least one like it.

dschwoerer · 2026-06-30T09:44:47Z

+      func(row, i.ym());
+      ++row;


More idiomatic C++:

Suggested change

func(row, i.ym());

++row;

func(row++, i.ym());

If accepted, it also needs to be changed above.

bendudson added 11 commits March 12, 2026 08:50

PetscOperators: FCI operators using PETSc

289ddf2

Read operators from the mesh in CSR format. PETSc handles communication pattern for MPI decomposition.

PetscOperators: Splitting out PetscMapping

c73c29a

Handles the mapping between mesh indices and PETSc global indices.

Mesh::get Arrays of int and BoutReal

724572a

Simple wrapper functions around the Options functionality, so that Array types can be read from the mesh file.

PetscOperators: Implement operators from CSR input

af65816

Reads row, column and weights from Mesh, constructs Petsc Mat objects. Working on composition and operation on Field3D objects.

PetscOperators: Operate on Field3Ds

8a3af0c

Map from Field3D, including yup/ydown boundaries, to PETSc Vec, perform matrix-vector product, and map back to a Field3D.

PetscOperators: Add documentation and error checking

2b71259

Docstrings to explain PetscMapping, PetscOperator and PetscOperators classes. Extend `BOUT_DO_PETSC` macro to print error message, file and line number.

PetscOperators: mapping fixes and unit tests

8ba9503

Some care needed with yup/down and yp/ym indices.

PetscOperator: Add diagonal and transpose methods

16e2768

`diagonal` creates diagonal matrices from Field3D. `transpose` creates a new matrix that is a transpose of a PetscOperator. This is useful for converting gradients into divergences.

PetscOperators: map matrix rows in boundary if present

5182ac4

Allows input operators to contain entries in boundary rows. This is important so that the gradient and divergence operators depend on boundary cell values.

Merge branch 'next' into petsc-operators

1e7aabe

github-actions Bot reviewed Mar 16, 2026

View reviewed changes

PetscOperators: Add BOUT_HAS_PETSC guards

68428b3

Not available unless configured with PETSc.

github-actions Bot reviewed Mar 16, 2026

View reviewed changes

bendudson added 2 commits March 16, 2026 16:46

PetscMapping: Add localSize and globalSize

2c45584

Check number of cells against global mapping size, not local size.

PetscOperator: Handle standard CSR and variations

d859712

Handle non-standard CSR inputs: - Negative indices for missing rows - Missing final entry in rows array

github-actions Bot reviewed Mar 17, 2026

View reviewed changes

bendudson mentioned this pull request Mar 17, 2026

Calculate interpolation weight matrix boutproject/zoidberg#40

Merged

bendudson added 2 commits March 16, 2026 21:41

PetscOperators: Fix unit tests

350eba5

Changes PetscMapping::size() to PetscMapping::localSize()

github-actions Bot reviewed Mar 18, 2026

View reviewed changes

bendudson added 3 commits March 18, 2026 11:13

PetscOperators::Parallel::Div_par_K_Grad_par

d5e1218

This is a function rather than a PetscOperator because it requires at least two matrix-vector products. This implementation requires six Mat-Vec products that could probably be optimized.

WIP: Sub-cell sampling

9a44ff8

Trace multiple field lines per cell, to construct smoother maps. A complication is cells that partially intersect a boundary.

fci-operators example: Don't include mesh file

3992b8c

Mesh not available in repo, only generation script.

Merge branch 'next' of https://github.com/boutproject/BOUT-dev into p…

1ffc592

…etsc-operators

github-actions Bot reviewed Jun 3, 2026

View reviewed changes

Comment thread tests/integrated/test-petsc-operators/test_petsc_operators.cxx

Comment thread tests/integrated/test-petsc-operators/test_petsc_operators.cxx

petsc-operators: Include guards

d4caff6

dschwoerer added 2 commits June 17, 2026 11:01

Merge branch 'next' of https://github.com/boutproject/BOUT-dev into p…

3250fcc

…etsc-operators

Fix off by one error

665c1c0

github-actions Bot reviewed Jun 17, 2026

View reviewed changes

Comment thread src/mesh/petsc_operators.cxx

Add missing headers

8019a94

dschwoerer force-pushed the petsc-operators branch 2 times, most recently from 562e38c to 8019a94 Compare June 18, 2026 07:23

github-actions Bot reviewed Jun 18, 2026

View reviewed changes

Comment thread include/bout/petsc_operators.hxx Outdated

dschwoerer added 3 commits June 18, 2026 09:34

Apply clang-tidy improvements

325d03b

FCI and 2D metrics don't mix well

e9d5656

New zoidberg has a check for inner/outer convention

c4a788d

dschwoerer mentioned this pull request Jun 24, 2026

Fix off by one error #3406

Merged

Merge branch 'next' into petsc-operators

cf8cec0

github-actions Bot reviewed Jun 24, 2026

View reviewed changes

bendudson added 5 commits June 29, 2026 12:04

Merge branch 'next' into petsc-operators

260c916

Fixes and tidying

83f1a56

- Arithmetic expression on FieldMetrics - Check for null source in Mesh::get - Use BOUT_USE_METRIC_3D rather than BOUT_HAS_METRIC_3D

PetscOperator::assembleFromCSR preallocate

ca092e8

Iterate over elements to preallocate the PETSc matrix before setting element values.

test-petsc-operators: Fixes

b5dc10b

Changes in `next` mean communicating Bxy throws exception (disabled calc parallel slices).

Merge branch 'next' into petsc-operators

3d0f5d2

bendudson requested review from ZedThree and dschwoerer June 30, 2026 02:59

Manual page on PetscOperator

9c53d77

dschwoerer approved these changes Jun 30, 2026

View reviewed changes

bendudson merged commit 2384555 into next Jun 30, 2026
23 checks passed

bendudson deleted the petsc-operators branch June 30, 2026 15:42

	void PetscOperator::copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec) {
	void PetscOperator::copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec) {

	Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const {
	Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) {

Uh oh!

Conversation

bendudson commented Mar 16, 2026

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 16, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

ZedThree commented Mar 18, 2026

Uh oh!

github-actions Bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 18, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions Bot Mar 18, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

bendudson commented Mar 18, 2026

Uh oh!

ZedThree commented Mar 18, 2026

Uh oh!

bendudson commented Mar 18, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

bendudson commented Mar 18, 2026 •

edited

Loading