Merge Jacobian coloring into PETSc solver - #3161
Conversation
Sets options on a TS object using the options in the BOUT.inp [petsc] section.
Uses the TS interface https://petsc.org/release/manual/ts/ allowing many different time integration methods, timestep adaptation methods etc. to be tried.
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
| // Compute IJacobian = dF/dU + a dF/dUdot | ||
| // This is a dummy matrix that saves the shift. | ||
| // The shift is later used in the matrix-free preconditioner | ||
| PetscErrorCode solver_ijacobian(TS, BoutReal, Vec, Vec, PetscReal shift, Mat J, Mat Jpre, |
There was a problem hiding this comment.
warning: all parameters should be named in a function [readability-named-parameter]
| PetscErrorCode solver_ijacobian(TS, BoutReal, Vec, Vec, PetscReal shift, Mat J, Mat Jpre, | |
| PetscErrorCode solver_ijacobian(TS /*unused*/, BoutReal /*unused*/, Vec /*unused*/, Vec /*unused*/, PetscReal shift, Mat J, Mat Jpre, |
| // Compute IJacobian = dF/dU + a dF/dUdot | ||
| // This is a dummy matrix that saves the shift. | ||
| // The shift is later used in the matrix-free preconditioner | ||
| PetscErrorCode solver_ijacobian(TS, BoutReal, Vec, Vec, PetscReal shift, Mat J, Mat Jpre, |
There was a problem hiding this comment.
warning: no header providing "PetscReal" is directly included [misc-include-cleaner]
PetscErrorCode solver_ijacobian(TS, BoutReal, Vec, Vec, PetscReal shift, Mat J, Mat Jpre,
^- Read `interpolate` option, allowing solution interpolation to be disabled for TS methods that don't support it (e.g. `pseudo`). - Fix coloring method: The TS function should be used for coloring, not the RHS function directly. This is because SNES is solving a modified function that depends on the time step and TS method.
Uses TSSetMaxSNESFailures, that can be overridden with a petsc option.
Replaces CHKERRQ macro call with a wrapper around the PETSc function call.
| ierr = VecGetArray(u, &udata); | ||
| CHKERRQ(ierr); | ||
| // Save initial state to PETSc Vec | ||
| BoutReal* udata; // Pointer to data array in vector u. |
There was a problem hiding this comment.
warning: variable 'udata' is not initialized [cppcoreguidelines-init-variables]
| BoutReal* udata; // Pointer to data array in vector u. | |
| BoutReal* udata = nullptr; // Pointer to data array in vector u. |
- Allow Jacobian and preconditioner to persist across SNES solves. This is useful because many steps can typically be taken using the same preconditioner. - Allow solver to recover from time step rejections by setting TSSetMaxStepRejections. This doesn't seem to work as expected with the `pseudo` TS method.
| CHKERRQ(ierr); | ||
| #endif | ||
| // Recover from step rejections | ||
| PetscCall(TSSetMaxStepRejections(ts, PETSC_UNLIMITED)); |
There was a problem hiding this comment.
warning: use of undeclared identifier 'PETSC_UNLIMITED' [clang-diagnostic-error]
PetscCall(TSSetMaxStepRejections(ts, PETSC_UNLIMITED));
^| } else { | ||
| ierr = TSMonitorSet(ts, PetscMonitor, this, nullptr); | ||
| CHKERRQ(ierr); | ||
| ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP); |
There was a problem hiding this comment.
warning: no header providing "TS_EXACTFINALTIME_MATCHSTEP" is directly included [misc-include-cleaner]
ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
^| ierr = SNESSetTolerances(snes, abstol, reltol, PETSC_DEFAULT, PETSC_DEFAULT, | ||
| PETSC_DEFAULT); | ||
| CHKERRQ(ierr); | ||
| PetscCall(TSMonitorSet(ts, PetscMonitor, this, nullptr)); |
There was a problem hiding this comment.
warning: no header providing "TSMonitorSet" is directly included [misc-include-cleaner]
PetscCall(TSMonitorSet(ts, PetscMonitor, this, nullptr));
^| } else { | ||
| // Use finite difference approximation | ||
| ierr = MatCreateSNESMF(snes, &Jmf); | ||
| ierr = TSGetSNES(ts, &snes); |
There was a problem hiding this comment.
warning: no header providing "TSGetSNES" is directly included [misc-include-cleaner]
ierr = TSGetSNES(ts, &snes);
^| CHKERRQ(ierr); | ||
| // Line search | ||
| if (line_search_type != "default") { | ||
| SNESLineSearch linesearch; |
There was a problem hiding this comment.
warning: variable 'linesearch' is not initialized [cppcoreguidelines-init-variables]
| SNESLineSearch linesearch; | |
| SNESLineSearch linesearch = nullptr; |
| PCSetType(pc, pc_type.c_str()); | ||
|
|
||
| if (pc_type == "hypre") { | ||
| #if PETSC_HAVE_HYPRE |
There was a problem hiding this comment.
warning: no header providing "PETSC_HAVE_HYPRE" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.cxx:26:
+ #include <petscconf.h>| if (pc_type == "hypre") { | ||
| #if PETSC_HAVE_HYPRE | ||
| // Set the type of hypre preconditioner | ||
| PCHYPRESetType(pc, pc_hypre_type.c_str()); |
There was a problem hiding this comment.
warning: no header providing "PCHYPRESetType" is directly included [misc-include-cleaner]
PCHYPRESetType(pc, pc_hypre_type.c_str());
^|
|
||
| ////////////////////////////////////////////////// | ||
| // Get the local indices by starting at 0 | ||
| Field3D index = globalIndex(0); |
There was a problem hiding this comment.
warning: no header providing "Field3D" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.cxx:26:
+ #include "bout/field3d.hxx"|
|
||
| output_progress.write("Setting Jacobian matrix sizes\n"); | ||
|
|
||
| const int n2d = f2d.size(); |
There was a problem hiding this comment.
warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]
const int n2d = f2d.size();
^| output_progress.write("Setting Jacobian matrix sizes\n"); | ||
|
|
||
| const int n2d = f2d.size(); | ||
| const int n3d = f3d.size(); |
There was a problem hiding this comment.
warning: narrowing conversion from 'size_type' (aka 'unsigned long') to signed type 'int' is implementation-defined [bugprone-narrowing-conversions]
const int n3d = f3d.size();
^Introduced in PETSc 3.22 (probably), define PETSC_UNLIMITED if it's undefined. Include more headers for used PETSc functionality.
| if (bout::globals::mpi->MPI_Allreduce(&local_N, &neq, 1, MPI_INT, MPI_SUM, | ||
| // Get total problem size | ||
| int neq; | ||
| if (bout::globals::mpi->MPI_Allreduce(&nlocal, &neq, 1, MPI_INT, MPI_SUM, |
There was a problem hiding this comment.
warning: implicit conversion 'int' -> 'bool' [readability-implicit-bool-conversion]
src/solver/impls/petsc/petsc.cxx:283:
- BoutComm::get())) {
+ BoutComm::get()) != 0) {Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
|
|
||
| // Star pattern | ||
| for (const auto& [x_off, y_off] : xy_offsets) { | ||
| int xi = x + x_off; |
There was a problem hiding this comment.
warning: variable 'xi' of type 'int' can be declared 'const' [misc-const-correctness]
| int xi = x + x_off; | |
| int const xi = x + x_off; |
| // Star pattern | ||
| for (const auto& [x_off, y_off] : xy_offsets) { | ||
| int xi = x + x_off; | ||
| int yi = y + y_off; |
There was a problem hiding this comment.
warning: variable 'yi' of type 'int' can be declared 'const' [misc-const-correctness]
| int yi = y + y_off; | |
| int const yi = y + y_off; |
| { | ||
| // Test if the matrix is symmetric | ||
| // Values are 0 or 1 so tolerance (1e-5) shouldn't matter | ||
| PetscBool symmetric; |
There was a problem hiding this comment.
warning: no header providing "PetscBool" is directly included [misc-include-cleaner]
PetscBool symmetric;
^| { | ||
| // Test if the matrix is symmetric | ||
| // Values are 0 or 1 so tolerance (1e-5) shouldn't matter | ||
| PetscBool symmetric; |
There was a problem hiding this comment.
warning: variable 'symmetric' is not initialized [cppcoreguidelines-init-variables]
PetscBool symmetric;
^| if (ksptype) { | ||
| output_info.write("KSP Type : {}\n", ksptype); | ||
| } | ||
| PCType pctype; |
There was a problem hiding this comment.
warning: no header providing "PCType" is directly included [misc-include-cleaner]
PCType pctype;
^| if (ksptype) { | ||
| output_info.write("KSP Type : {}\n", ksptype); | ||
| } | ||
| PCType pctype; |
There was a problem hiding this comment.
warning: variable 'pctype' is not initialized [cppcoreguidelines-init-variables]
| PCType pctype; | |
| PCType pctype = nullptr; |
| output_info.write("KSP Type : {}\n", ksptype); | ||
| } | ||
| PCType pctype; | ||
| PCGetType(pc, &pctype); |
There was a problem hiding this comment.
warning: no header providing "PCGetType" is directly included [misc-include-cleaner]
PCGetType(pc, &pctype);
^| } | ||
| PCType pctype; | ||
| PCGetType(pc, &pctype); | ||
| if (pctype) { |
There was a problem hiding this comment.
warning: implicit conversion 'PCType' (aka 'const char *') -> 'bool' [readability-implicit-bool-conversion]
| if (pctype) { | |
| if (pctype != nullptr) { |
| PetscFunctionBegin; | ||
|
|
||
| // Load state from PETSc | ||
| const BoutReal* udata_array; |
There was a problem hiding this comment.
warning: variable 'udata_array' is not initialized [cppcoreguidelines-init-variables]
| const BoutReal* udata_array; | |
| const BoutReal* udata_array = nullptr; |
Runs SUNDIALS through PETSc, but preconditioning slows down the simulation. Something wrong with the Jacobian calculation perhaps. Some tidying, adding header files etc., to address Clang Tidy comments.
| static PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { | ||
| // Get the context | ||
| void* ctx = nullptr; | ||
| PetscCall(PCShellGetContext(pc, &ctx)); |
There was a problem hiding this comment.
warning: multilevel pointer conversion from 'void **' to 'void *', please use explicit cast [bugprone-multi-level-implicit-pointer-conversion]
PetscCall(PCShellGetContext(pc, &ctx));
^| while (s->next_output <= t && s->next_output <= tfinal) { | ||
| BoutReal output_time = t; | ||
| if (s->interpolate) { | ||
| int ierr = TSInterpolate(ts, s->next_output, interpolatedX); |
There was a problem hiding this comment.
warning: variable 'ierr' of type 'int' can be declared 'const' [misc-const-correctness]
| int ierr = TSInterpolate(ts, s->next_output, interpolatedX); | |
| int const ierr = TSInterpolate(ts, s->next_output, interpolatedX); |
| // Re-calculate the coloring | ||
| MatColoring coloring = NULL; | ||
| MatColoringCreate(Jfd, &coloring); | ||
| MatColoringSetType(coloring, MATCOLORINGGREEDY); |
There was a problem hiding this comment.
warning: no header providing "MATCOLORINGGREEDY" is directly included [misc-include-cleaner]
MatColoringSetType(coloring, MATCOLORINGGREEDY);
^| // Re-calculate the coloring | ||
| MatColoring coloring = NULL; | ||
| MatColoringCreate(Jfd, &coloring); | ||
| MatColoringSetType(coloring, MATCOLORINGGREEDY); |
There was a problem hiding this comment.
warning: no header providing "MatColoringSetType" is directly included [misc-include-cleaner]
MatColoringSetType(coloring, MATCOLORINGGREEDY);
^| MatColoring coloring = NULL; | ||
| MatColoringCreate(Jfd, &coloring); | ||
| MatColoringSetType(coloring, MATCOLORINGGREEDY); | ||
| MatColoringSetFromOptions(coloring); |
There was a problem hiding this comment.
warning: no header providing "MatColoringSetFromOptions" is directly included [misc-include-cleaner]
MatColoringSetFromOptions(coloring);
^This is fiddly because TS doesn't create a SNES object when TSType is TSSUNDIALS. The nonlinear function being solved in sundials is similar to the PETSc convention, but scaled by the "gamma" / shift factor.
bendudson
changed the title
| Vec u{nullptr}; ///< PETSc solution vector | ||
| TS ts{nullptr}; ///< PETSc timestepper object | ||
| SNES snes{nullptr}; ///< PETSc nonlinear solver object | ||
| KSP ksp{nullptr}; ///< PETSc linear solver |
There was a problem hiding this comment.
warning: no header providing "KSP" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.hxx:30:
+ #include <petscksp.h>| KSP ksp{nullptr}; ///< PETSc linear solver | ||
| Mat Jmf{nullptr}; ///< Matrix Free Jacobian | ||
| Mat Jfd{nullptr}; ///< Finite Difference Jacobian | ||
| MatFDColoring fdcoloring{nullptr}; ///< Matrix coloring context |
There was a problem hiding this comment.
warning: no header providing "MatFDColoring" is directly included [misc-include-cleaner]
MatFDColoring fdcoloring{nullptr}; ///< Matrix coloring context
^The function signature seems to vary between PETSc versions in ways that break `reinterpret_cast` (and `bout::cast_MatFDColoringFn` in `include/bout/petsc_interface.hxx`).
This regrettable necessity follows the PETSc examples. Function pointers are cast to `PetscErrorCode (*)(void)`, which means that any future change to the signature will not be caught but will give undefined behavior. May future maintainers forgive me. Some drive-by tidying to appease Clang-tidy.
| #include <bout/output.hxx> | ||
| namespace { | ||
| // PETSc callback function for matrix-free preconditioner | ||
| PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { |
There was a problem hiding this comment.
warning: no header providing "PC" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.cxx:26:
+ #include <petscpctypes.h>|
|
||
| extern PetscErrorCode solver_f(TS ts, BoutReal t, Vec globalin, Vec globalout, | ||
| void* f_data); | ||
| PetscReal tfinal; |
There was a problem hiding this comment.
warning: variable 'tfinal' is not initialized [cppcoreguidelines-init-variables]
src/solver/impls/petsc/petsc.cxx:48:
- #include <petsc.h>
+ #include <math.h>
+ #include <petsc.h>| PetscReal tfinal; | |
| PetscReal tfinal = NAN; |
| MatFDColoringSetFromOptions(fdcoloring); | ||
| MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); | ||
| ISColoringDestroy(&iscoloring); | ||
|
|
There was a problem hiding this comment.
warning: no header providing "ISColoringDestroy" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.cxx:26:
+ #include <petscis.h>|
|
||
| /// Wrapper for the RHS function | ||
| PetscErrorCode rhs(TS ts, PetscReal t, Vec globalin, Vec globalout); | ||
| PetscErrorCode rhs(BoutReal t, Vec udata, Vec dudata, bool linear); |
There was a problem hiding this comment.
warning: no header providing "BoutReal" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.hxx:28:
- #include "bout/build_defines.hxx"
+ #include "bout/bout_types.hxx"
+ #include "bout/build_defines.hxx"|
|
||
| /// Wrapper for the RHS function | ||
| PetscErrorCode rhs(TS ts, PetscReal t, Vec globalin, Vec globalout); | ||
| PetscErrorCode rhs(BoutReal t, Vec udata, Vec dudata, bool linear); |
There was a problem hiding this comment.
warning: no header providing "PetscErrorCode" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.hxx:30:
+ #include <petscsystypes.h>|
|
||
| /// Wrapper for the RHS function | ||
| PetscErrorCode rhs(TS ts, PetscReal t, Vec globalin, Vec globalout); | ||
| PetscErrorCode rhs(BoutReal t, Vec udata, Vec dudata, bool linear); |
There was a problem hiding this comment.
warning: no header providing "Vec" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.hxx:30:
+ #include <petscvec.h>Getting MatFDColoringSetFunction to work across PETSc versions seems to be unreasonably difficult.
I put on my robe and wizard hat...
Tests pass with beuler/bdf/cn ts_types if tolerance is sufficiently tight.
commented
Nov 3, 2025
|
Future work might be to pull out a common colouring helper function? |
Takes the coloring code from the
snes/beulersolver, and merges it into thepetscsolver. This uses the PETScTSinterface, enabling use of many time-integration algorithms, including sundials/cvode.PETSc TS manual: https://petsc.org/release/manual/ts/
Use with e.g.
solver:type=petsc solver:ts_type=beulerOther
ts_typechoices includebdf,pseudo(https://petsc.org/release/manualpages/TS/TSPSEUDO/)See list of methods here: https://petsc.org/release/overview/integrator_table/#integrator-table
Other solver settings include:
adapt_type, that sets the adaptive timestep method https://petsc.org/release/manualpages/TS/TSAdaptType/snes_type, the nonlinear solver methodksp_type, the linear solver methodpc_type, the preconditioner methoduse e.g.
petsc:ts_monitor,petsc:snes_monitorto print diagnostics;petsc:helpto print the available PETSc options.The combination of
bdfmethod and STRUMPACK direct solver seems to be quite effective:That solves in 2 1/2 minutes a 2D transport problem that took CVODE (without preconditioning) 2 1/2 hours.
If PETSc is compiled with sundials (configure flag
--download-sundials) then the Jacobian coloring and PETSc preconditioners can be used with CVODE (2.5.0):This solved the same 2D transport problem in 1 1/2 minutes.
There are many (many!) settings to try and tune.