Convert `LaplacePetsc` to use `PetscMatrix` by ZedThree · Pull Request #3293 · boutproject/BOUT-dev

ZedThree · 2026-03-02T11:34:18Z

As well as making this implementation a lot cleaner, this is necessary for the Z parallelisation so that we can use the GlobalIndexer, otherwise it's really horrible to work out the domain decomposition (as the grid is Z-first and the processors are X-first, and the matrix may have different number of rows on each processor due to the X boundaries).

I had to remove the TRACE macro from the PetscMatrix as this gets called in a hot loop ((nx * nz)^2 times!) and the string formatting really adds up -- yet more reason to remove MsgStack?

In debug builds, this has a bit of a weird effect on performance:

next:

Timer name	Total time (s)	% of top	Hits	Mean time/hit (s)
petscsolve	4.36366336	1.00	6	0.727277226
petscsetup	1.50196716	0.34	6	0.250327859

this PR:

Timer name	Total time (s)	% of top	Hits	Mean time/hit (s)
petscsolve	3.76327226	1.00	6	0.627212044
petscsetup	2.25778754	0.60	6	0.376297923

The solve time goes down, but the setup time goes up.

The string formatting done in a hot loop ends up being _very_ expensive, 10x slowdown

github-actions

clang-tidy made some suggestions

github-actions · 2026-03-02T11:38:26Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+// TODO: handle fourth order
+auto set_stencil(const Mesh& localmesh, bool fourth_order) {
+  OperatorStencil<IndPerp> stencil;
+  IndexOffset<IndPerp> zero;


warning: variable 'zero' of type 'IndexOffset' (aka 'IndexOffset<SpecificInd<IND_TYPE::IND_PERP>>') can be declared 'const' [misc-const-correctness]

Suggested change

IndexOffset<IndPerp> zero;

IndexOffset<IndPerp> const zero;

github-actions · 2026-03-02T11:38:26Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+    });
+  }
+
+  std::vector offsetsVec(offsets.begin(), offsets.end());


warning: variable 'offsetsVec' of type 'std::vector<IndexOffset<SpecificInd<IND_TYPE::IND_PERP>>, std::allocator<IndexOffset<SpecificInd<IND_TYPE::IND_PERP>>>>' (aka 'vector<IndexOffset<SpecificInd<IND_TYPE::IND_PERP>>, std::allocator<IndexOffset<SpecificInd<IND_TYPE::IND_PERP>>>>') can be declared 'const' [misc-const-correctness]

Suggested change

std::vector offsetsVec(offsets.begin(), offsets.end());

std::vector const offsetsVec(offsets.begin(), offsets.end());

github-actions · 2026-03-02T11:38:26Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+    }
+  }
+
+  stencil.add([](IndPerp UNUSED(ind)) -> bool { return true; }, {zero});


warning: no header providing "UNUSED" is directly included [misc-include-cleaner]

src/invert/laplace/impls/petsc/petsc_laplace.cxx:27:

+ #include "bout/unused.hxx"

github-actions · 2026-03-02T11:38:26Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

      fourth_order(
          (*opts)["fourth_order"].doc("Use fourth order stencil").withDefault(false)),
-      lib(opts) {
+      indexer(std::make_shared<GlobalIndexer<FieldPerp>>(


warning: no header providing "std::make_shared" is directly included [misc-include-cleaner]

src/invert/laplace/impls/petsc/petsc_laplace.cxx:27:

+ #include <memory>

github-actions · 2026-03-02T11:38:26Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

-      }
-    }
+    // NOTE: Only A0 is the A from setCoefA ()
+    const BoutReal A0 = A[i];


warning: 'A0' is confusable with 'AO' [misc-confusable-identifiers]

const BoutReal A0 = A[i]; ^

Additional context

/usr/lib/petscdir/petsc3.22/x86_64-linux-gnu-real/include/petscao.h:19: other declaration found here

typedef struct _p_AO *AO; ^

github-actions · 2026-03-02T11:38:27Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+
+    BoutReal dx = coords->dx[i];
+    BoutReal dx2 = SQ(dx);
+    BoutReal dz = coords->dz[i];


warning: variable 'dz' of type 'BoutReal' (aka 'double') can be declared 'const' [misc-const-correctness]

Suggested change

BoutReal dz = coords->dz[i];

BoutReal const dz = coords->dz[i];

github-actions · 2026-03-02T11:38:27Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+    BoutReal dx = coords->dx[i];
+    BoutReal dx2 = SQ(dx);
+    BoutReal dz = coords->dz[i];
+    BoutReal dz2 = SQ(dz);


warning: variable 'dz2' of type 'BoutReal' (aka 'double') can be declared 'const' [misc-const-correctness]

Suggested change

BoutReal dz2 = SQ(dz);

BoutReal const dz2 = SQ(dz);

github-actions · 2026-03-02T11:38:27Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+    BoutReal dx2 = SQ(dx);
+    BoutReal dz = coords->dz[i];
+    BoutReal dz2 = SQ(dz);
+    BoutReal dxdz = dx * dz;


warning: variable 'dxdz' of type 'BoutReal' (aka 'double') can be declared 'const' [misc-const-correctness]

Suggested change

BoutReal dxdz = dx * dz;

BoutReal const dxdz = dx * dz;

github-actions · 2026-03-02T11:38:27Z

tests/integrated/test-petsc_laplace/test_petsc_laplace.cxx

+
+  // Need this here to ensure PETSc isn't finalised until after the global mesh,
+  // otherwise we get problems from `MPI_Comm_free` on the X communicator
+  PetscLib lib{};


warning: variable 'lib' of type 'PetscLib' can be declared 'const' [misc-const-correctness]

Suggested change

PetscLib lib{};

PetscLib const lib{};

github-actions · 2026-03-02T11:38:27Z

tests/integrated/test-petsc_laplace_MAST-grid/test_petsc_laplace_MAST_grid.cxx


  BoutInitialise(argc, argv);
+
+  PetscLib lib{};


warning: variable 'lib' of type 'PetscLib' can be declared 'const' [misc-const-correctness]

Suggested change

PetscLib lib{};

PetscLib const lib{};

ZedThree · 2026-03-02T11:50:51Z

Hmm, the 3D build is failing with

[0]PETSC ERROR: Unable to find requested PC type mumps

even though I explicitly check that PETSc was compiled with mumps. I guess it's fine to always use sor

dschwoerer · 2026-03-06T13:31:56Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+      operator2D(i, i) = (25.0 / 12.0) * factor;
+      operator2D(i, i.xm(1)) = -4.0 * factor;
+      operator2D(i, i.xm(2)) = 3.0 * factor;
+      operator2D(i, i.xm(3)) = (-4.0 / 3.0) * factor;
+      operator2D(i, i.xm(4)) = (1.0 / 4.0) * factor;


There should be two points where this is used, but on one of them we should be able to use also i.xp() ?
Would that not be better?

dschwoerer · 2026-03-06T13:32:34Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

+    operator2D(l, l.xpp().zm()) = A3 / (18.0 * dxdz);
+    operator2D(l, l.xpp()) = (-1.0 / 12.0) * ((A1 / dx2) + (A4 / dx));
+    operator2D(l, l.xpp().zp()) = -1.0 * A3 / (18.0 * dxdz);
+    operator2D(l, l.xpp().zpp()) = A3 / (144.0 * dxdz);


We could pull out the A3/dxdz and similar.

I guess MMS tests are checking whether all the factors are correct?

-1.0 * A3 / (18.0 * dxdz) could be simplified to

auto const A3fac = A3 / dxdz` A3fac / -18.0

Not sure the compiler will do that, as that might only be allowed with -Ofast and similar.

dschwoerer · 2026-03-06T13:38:05Z

src/invert/laplace/impls/petsc/petsc_laplace.cxx

  PetscFunctionReturn(laplace->precon(x, y)); // NOLINT
 }
+
+// TODO: handle fourth order


That seems to be done now?

ZedThree and others added 6 commits March 2, 2026 10:21

Convert LaplacePetsc to use PetscMatrix

d60e541

Remove TRACE macro from PetscMatrix

2305c08

The string formatting done in a hot loop ends up being _very_ expensive, 10x slowdown

Include all headers in petsc interface

e22b39e

Use preconditioner for petsc MAST grids test

06ea972

[bot] Apply format changes

c07af62

[bot] Add last format changes commit to ignore file

57bb02d

github-actions bot reviewed Mar 2, 2026

View reviewed changes

dschwoerer reviewed Mar 6, 2026

View reviewed changes

	std::vector offsetsVec(offsets.begin(), offsets.end());
	std::vector const offsetsVec(offsets.begin(), offsets.end());

	BoutReal dz = coords->dz[i];
	BoutReal const dz = coords->dz[i];

Conversation

ZedThree commented Mar 2, 2026

Uh oh!

github-actions bot left a comment

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

github-actions bot Mar 2, 2026

Choose a reason for hiding this comment

Uh oh!

ZedThree commented Mar 2, 2026

Uh oh!

dschwoerer Mar 6, 2026

Choose a reason for hiding this comment

Uh oh!

dschwoerer Mar 6, 2026

Choose a reason for hiding this comment

Uh oh!

dschwoerer Mar 6, 2026

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants