From 741b1df7c693db1c0970044206d25749801f4291 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 14 Aug 2025 20:51:21 -0700 Subject: [PATCH 01/18] PetscLib: Add setOptionsFromInputFile(TS& ts) Sets options on a TS object using the options in the BOUT.inp [petsc] section. --- include/bout/petsclib.hxx | 7 +++++-- src/sys/petsclib.cxx | 6 ++++++ 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/include/bout/petsclib.hxx b/include/bout/petsclib.hxx index 2008671286..760b18ff91 100644 --- a/include/bout/petsclib.hxx +++ b/include/bout/petsclib.hxx @@ -22,9 +22,9 @@ * so it *must* be included before *any* PETSc header. * ************************************************************************** - * Copyright 2012 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2012 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -113,6 +113,9 @@ public: /// was passed to the constructor. void setOptionsFromInputFile(SNES& snes); + /// Set options for a TS time integrator + void setOptionsFromInputFile(TS& ts); + /*! * Force cleanup. This will call PetscFinalize, printing a warning * if any instances of PetscLib still exist diff --git a/src/sys/petsclib.cxx b/src/sys/petsclib.cxx index f1cf1a9d1b..1830db89d9 100644 --- a/src/sys/petsclib.cxx +++ b/src/sys/petsclib.cxx @@ -119,6 +119,12 @@ void PetscLib::setOptionsFromInputFile(SNES& snes) { BOUT_DO_PETSC(SNESSetFromOptions(snes)); } +void PetscLib::setOptionsFromInputFile(TS& ts) { + BOUT_DO_PETSC(TSSetOptionsPrefix(ts, options_prefix.c_str())); + + BOUT_DO_PETSC(TSSetFromOptions(ts)); +} + void PetscLib::cleanup() { BOUT_OMP_SAFE(critical(PetscLib)) { From 33b984291a799a92466ddc284c2dd78776f94a06 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 14 Aug 2025 20:55:30 -0700 Subject: [PATCH 02/18] PETSc solver: Merge SNES Jacobian coloring with TS interface Uses the TS interface https://petsc.org/release/manual/ts/ allowing many different time integration methods, timestep adaptation methods etc. to be tried. --- src/solver/impls/petsc/petsc.cxx | 1518 +++++++++++++++--------------- src/solver/impls/petsc/petsc.hxx | 102 +- 2 files changed, 789 insertions(+), 831 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index dd63811f0f..22d8594c28 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -2,9 +2,9 @@ * Interface to PETSc solver * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -32,53 +32,230 @@ #include #include +#include +#include +#include +#include +#include -#include +#include +#include -#include // Cell interpolation -#include #include -extern PetscErrorCode solver_f(TS ts, BoutReal t, Vec globalin, Vec globalout, - void* f_data); +#include "petscsnes.h" + +class ColoringStencil { +private: + bool static isInSquare(int const i, int const j, int const n_square) { + return std::abs(i) <= n_square && std::abs(j) <= n_square; + } + bool static isInCross(int const i, int const j, int const n_cross) { + if (i == 0) { + return std::abs(j) <= n_cross; + } + if (j == 0) { + return std::abs(i) <= n_cross; + } + return false; + } + bool static isInTaxi(int const i, int const j, int const n_taxi) { + return std::abs(i) + std::abs(j) <= n_taxi; + } + +public: + auto static getOffsets(int n_square, int n_taxi, int n_cross) { + ASSERT2(n_square >= 0 && n_cross >= 0 && n_taxi >= 0 + && n_square + n_cross + n_taxi > 0); + auto inside = [&](int i, int j) { + return isInSquare(i, j, n_square) || isInTaxi(i, j, n_taxi) + || isInCross(i, j, n_cross); + }; + std::vector> xy_offsets; + auto loop_bound = std::max({n_square, n_taxi, n_cross}); + for (int i = -loop_bound; i <= loop_bound; ++i) { + for (int j = -loop_bound; j <= loop_bound; ++j) { + if (inside(i, j)) { + xy_offsets.emplace_back(i, j); + } + } + } + return xy_offsets; + } +}; + +// PETSc callback function for matrix-free preconditioner +static PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { + // Get the context + void* ctx; + int ierr = PCShellGetContext(pc, &ctx); + CHKERRQ(ierr); + // Run the preconditioner + PetscFunctionReturn(static_cast(ctx)->pre(x, y)); +} + +// PETSc callback function, that evaluates the nonlinear +// function being integrated by TS. +static PetscErrorCode solver_rhs(TS UNUSED(ts), BoutReal t, Vec globalin, Vec globalout, + void* f_data) { + PetscFunctionBegin; + auto* s = static_cast(f_data); + s->rhs(t, globalin, globalout, false); + PetscFunctionReturn(0); +} + +// PETSc callback function for finite differencing. +// May use a linearised form of the RHS function. +static PetscErrorCode solver_rhs_differencing(void* ctx, Vec globalin, Vec globalout) { + PetscFunctionBegin; + auto* s = static_cast(ctx); + s->rhs_differencing(globalin, globalout); + PetscFunctionReturn(0); +} + +// PETSc callback function for coloring. +// May use a linearised form of the RHS function. +static PetscErrorCode solver_rhs_coloring(void* UNUSED(snes), Vec globalin, Vec globalout, + void* ctx) { + PetscFunctionBegin; + auto* s = static_cast(ctx); + s->rhs_differencing(globalin, globalout); + PetscFunctionReturn(0); +} -#if PETSC_VERSION_GE(3, 5, 0) -extern PetscErrorCode solver_rhsjacobian(TS ts, BoutReal t, Vec globalin, Mat J, Mat Jpre, - void* f_data); -extern PetscErrorCode solver_ijacobianfd(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, - void*); +// 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, + void* ctx) { + auto* solver = static_cast(ctx); + solver->shift = shift; + + PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY)); + PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY)); + if (J != Jpre) { + PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY)); + PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY)); + } + + PetscFunctionReturn(0); +} + +// This function is called by the TS object every internal timestep +// It is responsible for triggering interpolation and output at +// the output time interval. +PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, void* ctx) { + PetscFunctionBegin; + + auto* s = static_cast(ctx); + if (t < s->next_output) { + // Not reached output time yet => return + PetscFunctionReturn(0); + } + + PetscErrorCode ierr; + PetscReal tfinal; + static int i = 0; + +#if PETSC_VERSION_GE(3, 8, 0) + ierr = TSGetMaxTime(ts, &tfinal); + CHKERRQ(ierr); #else -extern PetscErrorCode solver_rhsjacobian(TS ts, BoutReal t, Vec globalin, Mat* J, - Mat* Jpre, MatStructure* str, void* f_data); -extern PetscErrorCode solver_ijacobianfd(TS, PetscReal, Vec, Vec, PetscReal, Mat*, Mat*, - MatStructure*, void*); + ierr = TSGetDuration(ts, nullptr, &tfinal); + CHKERRQ(ierr); #endif -extern PetscErrorCode solver_if(TS, BoutReal, Vec, Vec, Vec, void*); + // Duplicate the solution vector X into a work vector + Vec interpolatedX; + ierr = VecDuplicate(X, &interpolatedX); + CHKERRQ(ierr); -/// KSP preconditioner PCShell routines for physics preconditioners -extern PetscErrorCode PhysicsPCApply(PC, Vec x, Vec y); -extern PetscErrorCode PhysicsJacobianApply(Mat J, Vec x, Vec y); -extern PetscErrorCode PhysicsSNESApply(SNES, Vec); + // The internal timestepper may have stepped over multiple output times + while (s->next_output <= t && s->next_output <= tfinal) { + if (s->interpolate) { + ierr = TSInterpolate(ts, s->next_output, interpolatedX); + CHKERRQ(ierr); + } + + // Place the interpolated values into the global variables + const PetscScalar* x; + ierr = VecGetArrayRead(interpolatedX, &x); + CHKERRQ(ierr); + s->load_vars(const_cast(x)); + ierr = VecRestoreArrayRead(interpolatedX, &x); + CHKERRQ(ierr); + + if (s->call_monitors(s->next_output, i++, s->getNumberOutputSteps())) { + PetscFunctionReturn(1); + } + + s->next_output += s->getOutputTimestep(); + } + + // Done with vector, so destroy it + ierr = VecDestroy(&interpolatedX); + CHKERRQ(ierr); + + PetscFunctionReturn(0); +} PetscSolver::PetscSolver(Options* opts) : Solver(opts), diagnose( (*options)["diagnose"].doc("Enable some diagnostic output").withDefault(false)), - adaptive( - (*options)["adaptive"].doc("Use adaptive timestepping").withDefault(false)), - use_precon((*options)["use_precon"] - .doc("Use user-supplied preconditioning function") - .withDefault(false)), - use_jacobian((*options)["use_jacobian"] - .doc("Use user-supplied Jacobian function") - .withDefault(false)), - abstol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)), - reltol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-5)), - adams_moulton((*options)["adams_moulton"] - .doc("Use Adams-Moulton implicit multistep method instead of BDF " - "(requires PETSc to have been built with SUNDIALS)") - .withDefault(false)), + user_precon((*options)["user_precon"] + .doc("Use user-supplied preconditioning function?") + .withDefault(false)), + atol((*options)["atol"].doc("Absolute tolerance").withDefault(1.0e-12)), + rtol((*options)["rtol"].doc("Relative tolerance").withDefault(1.0e-5)), + stol((*options)["stol"] + .doc("Convergence tolerance in terms of the norm of the change in " + "the solution between steps") + .withDefault(1e-8)), + maxnl((*options)["max_nonlinear_iterations"] + .doc("Maximum number of nonlinear iterations per SNES solve") + .withDefault(50)), + maxf((*options)["maxf"] + .doc("Maximum number of function evaluations per SNES solve") + .withDefault(10000)), + maxl((*options)["maxl"].doc("Maximum number of linear iterations").withDefault(20)), + ts_type( + (*options)["ts_type"].doc("PETSc time integrator type").withDefault("beuler")), + adapt_type((*options)["adapt_type"] + .doc("PETSc TSAdaptType timestep adaptation method") + .withDefault("basic")), + snes_type((*options)["snes_type"] + .doc("PETSc nonlinear solver method to use") + .withDefault("newtonls")), + ksp_type((*options)["ksp_type"] + .doc("Linear solver type. By default let PETSc decide (gmres)") + .withDefault("default")), + pc_type( + (*options)["pc_type"] + .doc("Preconditioner type. By default lets PETSc decide (ilu or bjacobi)") + .withDefault("default")), + pc_hypre_type((*options)["pc_hypre_type"] + .doc("hypre preconditioner type: euclid, pilut, parasails, " + "boomeramg, ams, ads") + .withDefault("pilut")), + line_search_type((*options)["line_search_type"] + .doc("Line search type: basic, bt, l2, cp, nleqerr") + .withDefault("default")), + matrix_free((*options)["matrix_free"] + .doc("Use matrix free Jacobian?") + .withDefault(false)), + matrix_free_operator((*options)["matrix_free_operator"] + .doc("Use matrix free Jacobian-vector operator?") + .withDefault(false)), + lag_jacobian((*options)["lag_jacobian"] + .doc("Re-use the Jacobian this number of SNES iterations") + .withDefault(50)), + use_coloring((*options)["use_coloring"] + .doc("Use matrix coloring to calculate Jacobian?") + .withDefault(true)), + kspsetinitialguessnonzero((*options)["kspsetinitialguessnonzero"] + .doc("Set the initial guess to be non-zero") + .withDefault(false)), start_timestep((*options)["start_timestep"] .doc("Initial internal timestep (defaults to output timestep)") .withDefault(getOutputTimestep())), @@ -87,9 +264,9 @@ PetscSolver::PetscSolver(Options* opts) PetscSolver::~PetscSolver() { VecDestroy(&u); - MatDestroy(&J); + MatDestroy(&Jfd); MatDestroy(&Jmf); - MatFDColoringDestroy(&matfdcoloring); + MatFDColoringDestroy(&fdcoloring); TSDestroy(&ts); } @@ -99,74 +276,45 @@ PetscSolver::~PetscSolver() { int PetscSolver::init() { PetscErrorCode ierr; - int neq; - MPI_Comm comm = PETSC_COMM_WORLD; - PetscMPIInt rank; TRACE("Initialising PETSc-dev solver"); - PetscFunctionBegin; - PetscLogEventRegister("PetscSolver::init", PETSC_VIEWER_CLASSID, &init_event); - PetscLogEventRegister("loop_vars", PETSC_VIEWER_CLASSID, &loop_event); - PetscLogEventRegister("solver_f", PETSC_VIEWER_CLASSID, &solver_event); - Solver::init(); - ierr = PetscLogEventBegin(init_event, 0, 0, 0, 0); - CHKERRQ(ierr); - output.write("Initialising PETSc-dev solver\n"); - ierr = bout::globals::mpi->MPI_Comm_rank(comm, &rank); - CHKERRQ(ierr); + int nlocal = getLocalN(); // Number of evolving variables on this processor - PetscInt local_N = getLocalN(); // Number of evolving variables on this processor - - /********** Get total problem size **********/ - 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, BoutComm::get())) { - output_error.write("\tERROR: MPI_Allreduce failed!\n"); - ierr = PetscLogEventEnd(init_event, 0, 0, 0, 0); - CHKERRQ(ierr); - PetscFunctionReturn(1); + throw BoutException("MPI_Allreduce failed!"); } - ierr = VecCreate(BoutComm::get(), &u); - CHKERRQ(ierr); - ierr = VecSetSizes(u, local_N, PETSC_DECIDE); + output_info.write("\t3d fields = {:d}, 2d fields = {:d} neq={:d}, local_N={:d}\n", + n3Dvars(), n2Dvars(), neq, nlocal); + + // Create a vector to contain the state + PetscCall(VecCreate(BoutComm::get(), &u)); + ierr = VecSetSizes(u, nlocal, PETSC_DECIDE); CHKERRQ(ierr); ierr = VecSetFromOptions(u); CHKERRQ(ierr); - ////////// SAVE INITIAL STATE TO PETSc VECTOR /////////// - // Set pointer to data array in vector u. - BoutReal* udata; - + // Save initial state to PETSc Vec + BoutReal* udata; // Pointer to data array in vector u. ierr = VecGetArray(u, &udata); CHKERRQ(ierr); save_vars(udata); ierr = VecRestoreArray(u, &udata); CHKERRQ(ierr); - PetscReal norm; - ierr = VecNorm(u, NORM_1, &norm); - CHKERRQ(ierr); - output_info << "initial |u| = " << norm << "\n"; - - PetscBool J_load; - char load_file[PETSC_MAX_PATH_LEN]; /* jacobian input file name */ - PetscBool J_write = PETSC_FALSE, J_slowfd = PETSC_FALSE; - // Create timestepper ierr = TSCreate(BoutComm::get(), &ts); CHKERRQ(ierr); ierr = TSSetProblemType(ts, TS_NONLINEAR); CHKERRQ(ierr); -#if PETSC_HAS_SUNDIALS - ierr = TSSetType(ts, TSSUNDIALS); - CHKERRQ(ierr); -#else - ierr = TSSetType(ts, TSRK); + ierr = TSSetType(ts, ts_type.c_str()); CHKERRQ(ierr); -#endif ierr = TSSetApplicationContext(ts, this); CHKERRQ(ierr); @@ -175,55 +323,38 @@ int PetscSolver::init() { Vec rhs_vec; ierr = VecDuplicate(u, &rhs_vec); CHKERRQ(ierr); - ierr = TSSetRHSFunction(ts, rhs_vec, solver_f, this); - CHKERRQ(ierr); - ierr = VecDestroy(&rhs_vec); + ierr = TSSetRHSFunction(ts, rhs_vec, solver_rhs, this); CHKERRQ(ierr); // Set up adaptive time-stepping TSAdapt adapt; ierr = TSGetAdapt(ts, &adapt); CHKERRQ(ierr); - if (adaptive) { - ierr = TSAdaptSetType(adapt, TSADAPTBASIC); - CHKERRQ(ierr); - } else { - ierr = TSAdaptSetType(adapt, TSADAPTNONE); - CHKERRQ(ierr); - } - - // Set default absolute/relative tolerances - ierr = TSSetTolerances(ts, abstol, nullptr, reltol, nullptr); + ierr = TSAdaptSetType(adapt, adapt_type.c_str()); CHKERRQ(ierr); -#if PETSC_VERSION_LT(3, 5, 0) - ierr = TSRKSetTolerance(ts, reltol); + // Set default absolute/relative tolerances + // Note: Vector atol and rtol not given + ierr = TSSetTolerances(ts, atol, nullptr, rtol, nullptr); CHKERRQ(ierr); -#endif - + if (ts_type == TSSUNDIALS) { #if PETSC_HAS_SUNDIALS - // Set Sundials tolerances - ierr = TSSundialsSetTolerance(ts, abstol, reltol); - CHKERRQ(ierr); - - // Select Sundials Adams-Moulton or BDF method - if (adams_moulton) { - output.write("\tUsing Adams-Moulton implicit multistep method\n"); - ierr = TSSundialsSetType(ts, SUNDIALS_ADAMS); - CHKERRQ(ierr); - } else { - output.write("\tUsing BDF method\n"); - ierr = TSSundialsSetType(ts, SUNDIALS_BDF); - CHKERRQ(ierr); - } + // The PETSc interface to SUNDIALS' CVODE + TSSundialsSetType(ts, SUNDIALS_BDF); + TSSundialsSetTolerance(ts, atol, rtol); +#else + throw BoutException("PETSc was not built with SUNDIALS. Reconfigure and build PETSc " + "with --download-sundials"); #endif + } - // Initial time and timestep + // Initial time of the simulation state ierr = TSSetTime(ts, simtime); CHKERRQ(ierr); + next_output = simtime; + ierr = TSSetTimeStep(ts, start_timestep); CHKERRQ(ierr); - next_output = simtime; // Total number of steps PetscInt total_steps = mxstep * getNumberOutputSteps(); @@ -245,377 +376,537 @@ int PetscSolver::init() { // Set the current solution ierr = TSSetSolution(ts, u); CHKERRQ(ierr); - - // Create RHSJacobian J - SNES snes, psnes; - KSP ksp, nksp; - PC pc, npc; - PCType pctype; - TSType tstype; - PetscBool pcnone = PETSC_TRUE; - - ierr = TSGetSNES(ts, &snes); - CHKERRQ(ierr); + // Allow TS to step over the final time + // Note: This does not affect intermediate outputs, that + // are always interpolated (in PetscMonitor) ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE); CHKERRQ(ierr); + ierr = TSMonitorSet(ts, PetscMonitor, this, nullptr); -#if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsGetBool(nullptr, nullptr, "-interpolate", &interpolate, nullptr); - CHKERRQ(ierr); -#else - ierr = PetscOptionsGetBool(nullptr, "-interpolate", &interpolate, nullptr); - CHKERRQ(ierr); -#endif + if (ts_type != TSSUNDIALS) { + // Get and configure the SNES nonlinear solver + // Note: SUNDIALS does not use PETSc SNES or KSP + // https://petsc.org/release/manual/ts/#using-sundials-from-petsc - // Check for -output_name to see if user specified a "performance" - // run, if they didn't then use the standard monitor function. TODO: - // use PetscFList -#if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsGetString(nullptr, nullptr, "-output_name", this->output_name, - sizeof this->output_name, &output_flag); - CHKERRQ(ierr); -#else - ierr = PetscOptionsGetString(nullptr, "-output_name", this->output_name, - sizeof this->output_name, &output_flag); - CHKERRQ(ierr); -#endif - - // If the output_name is not specified then use the standard monitor function - if (output_flag != 0U) { - ierr = SNESMonitorSet(snes, PetscSNESMonitor, this, nullptr); + ierr = TSGetSNES(ts, &snes); CHKERRQ(ierr); - } else { - ierr = TSMonitorSet(ts, PetscMonitor, this, nullptr); - CHKERRQ(ierr); - } - ierr = SNESSetTolerances(snes, abstol, reltol, PETSC_DEFAULT, PETSC_DEFAULT, - PETSC_DEFAULT); - CHKERRQ(ierr); - - // Matrix free Jacobian - - if (use_jacobian and hasJacobian()) { - // Use a user-supplied Jacobian function - ierr = MatCreateShell(comm, local_N, local_N, neq, neq, this, &Jmf); - CHKERRQ(ierr); - ierr = MatShellSetOperation(Jmf, MATOP_MULT, - reinterpret_cast(PhysicsJacobianApply)); - CHKERRQ(ierr); - ierr = TSSetIJacobian(ts, Jmf, Jmf, solver_ijacobian, this); - CHKERRQ(ierr); - } else { - // Use finite difference approximation - ierr = MatCreateSNESMF(snes, &Jmf); - CHKERRQ(ierr); - ierr = SNESSetJacobian(snes, Jmf, Jmf, MatMFFDComputeJacobian, this); + ierr = SNESSetType(snes, snes_type.c_str()); CHKERRQ(ierr); - } - - ierr = SNESGetKSP(snes, &ksp); - CHKERRQ(ierr); - ierr = KSPSetTolerances(ksp, reltol, abstol, PETSC_DEFAULT, PETSC_DEFAULT); - CHKERRQ(ierr); - - ierr = KSPGetPC(ksp, &pc); - CHKERRQ(ierr); + // Line search + if (line_search_type != "default") { + SNESLineSearch linesearch; + SNESGetLineSearch(snes, &linesearch); + SNESLineSearchSetType(linesearch, line_search_type.c_str()); + } - if (use_precon and hasPreconditioner()) { + // Set tolerances + // Note: TS should set SNES convergence tolerances + SNESSetTolerances(snes, + PETSC_DETERMINE, // atol + PETSC_DETERMINE, // rtol + PETSC_DETERMINE, // stol + maxnl, maxf); -#if PETSC_VERSION_GE(3, 5, 0) - ierr = SNESGetNPC(snes, &psnes); - CHKERRQ(ierr); -#else - ierr = SNESGetPC(snes, &psnes); - CHKERRQ(ierr); + // Force SNES to take at least one nonlinear iteration. + // This may prevent the solver from getting stuck in false steady state conditions +#if PETSC_VERSION_GE(3, 8, 0) + SNESSetForceIteration(snes, PETSC_TRUE); #endif - ierr = SNESGetKSP(psnes, &nksp); - CHKERRQ(ierr); - ierr = KSPGetPC(nksp, &npc); - CHKERRQ(ierr); - ierr = SNESSetType(psnes, SNESSHELL); - CHKERRQ(ierr); - ierr = SNESShellSetSolve(psnes, PhysicsSNESApply); - CHKERRQ(ierr); - // Use a user-supplied preconditioner + // Re-use Jacobian + // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian + // always updates its reference 'u' vector every nonlinear iteration + SNESSetLagJacobian(snes, lag_jacobian); + SNESSetLagJacobianPersists(snes, PETSC_FALSE); - // Tell PETSc we're using a "shell" preconditioner - ierr = PCSetType(pc, PCSHELL); - CHKERRQ(ierr); - - // Set routine for applying preconditioner - ierr = PCShellSetApply(pc, PhysicsPCApply); - CHKERRQ(ierr); + SNESSetLagPreconditionerPersists(snes, PETSC_FALSE); + SNESSetLagPreconditioner(snes, 1); // Rebuild when Jacobian is rebuilt - // Set context to this solver object - ierr = PCShellSetContext(pc, this); - CHKERRQ(ierr); + // Get and configure the KSP linear solver + SNESGetKSP(snes, &ksp); - // Set name of preconditioner - ierr = PCShellSetName(pc, "PhysicsPreconditioner"); - CHKERRQ(ierr); - - // Need a callback for IJacobian to get shift 'alpha' - ierr = TSSetIJacobian(ts, Jmf, Jmf, solver_ijacobian, this); - CHKERRQ(ierr); + if (ksp_type != "default") { + KSPSetType(ksp, ksp_type.c_str()); + } - // Use right preconditioner - ierr = KSPSetPCSide(ksp, PC_RIGHT); - CHKERRQ(ierr); + if (kspsetinitialguessnonzero) { + // Set the initial guess to be non-zero + KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); + } - } else { - // Default to no preconditioner - ierr = PCSetType(pc, PCNONE); - CHKERRQ(ierr); + KSPSetTolerances(ksp, + PETSC_DEFAULT, // rtol + PETSC_DEFAULT, // abstol + PETSC_DEFAULT, // dtol (divergence tolerance) + maxl); // Maximum number of iterations + + // Set up the Jacobian + // Note: We can have both matrix-free Jacobian-vector operator (Jmf) + // and a finite-difference Jacobian matrix (Jfd) + if (matrix_free or matrix_free_operator) { + /* + PETSc SNES matrix free Jacobian, using a different + operator for differencing. + + See PETSc examples + http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tests/ex7.c.html + and this thread: + http://lists.mcs.anl.gov/pipermail/petsc-users/2014-January/020075.html + + */ + MatCreateSNESMF(snes, &Jmf); + + // Set a function to be called for differencing + // This can be a linearised form of the SNES function + // Note: Unsure if setting this will interfere with TS object + // because the SNES Jacobian depends on the timestep + //MatMFFDSetFunction(Jmf, solver_rhs_differencing, this); + } } - ierr = TSSetFromOptions(ts); - CHKERRQ(ierr); // enable PETSc runtime options - - ierr = PCGetType(pc, &pctype); - CHKERRQ(ierr); - ierr = TSGetType(ts, &tstype); - CHKERRQ(ierr); - output.write("\tTS type {:s}, PC type {:s}\n", tstype, pctype); - ierr = PetscObjectTypeCompare(reinterpret_cast(pc), PCNONE, &pcnone); - CHKERRQ(ierr); - if (pcnone) { - ierr = PetscLogEventEnd(init_event, 0, 0, 0, 0); - CHKERRQ(ierr); - PetscFunctionReturn(0); - } - ierr = PetscObjectTypeCompare(reinterpret_cast(pc), PCSHELL, &pcnone); - CHKERRQ(ierr); - if (pcnone) { - ierr = PetscLogEventEnd(init_event, 0, 0, 0, 0); - CHKERRQ(ierr); - PetscFunctionReturn(0); + // Configure preconditioner + PC pc; + if (ts_type == TSSUNDIALS) { +#if PETSC_HAS_SUNDIALS + TSSundialsGetPC(ts, &pc); +#endif + } else { + // Get PC context from KSP + KSPGetPC(ksp, &pc); } + if (matrix_free) { + // Matrix-free preconditioner - // Create Jacobian matrix to be used by preconditioner - output_info << " Get Jacobian matrix at simtime " << simtime << "\n"; -#if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsGetString(nullptr, nullptr, "-J_load", load_file, - PETSC_MAX_PATH_LEN - 1, &J_load); - CHKERRQ(ierr); -#else - ierr = PetscOptionsGetString(nullptr, "-J_load", load_file, PETSC_MAX_PATH_LEN - 1, - &J_load); - CHKERRQ(ierr); -#endif - if (J_load) { - PetscViewer fd; - output_info << " Load Jmat ...local_N " << local_N << " neq " << neq << "\n"; - ierr = PetscViewerBinaryOpen(comm, load_file, FILE_MODE_READ, &fd); - CHKERRQ(ierr); - ierr = MatCreate(PETSC_COMM_WORLD, &J); - CHKERRQ(ierr); - //ierr = MatSetType(J, MATBAIJ);CHKERRQ(ierr); - ierr = MatSetSizes(J, local_N, local_N, PETSC_DECIDE, PETSC_DECIDE); - CHKERRQ(ierr); - ierr = MatSetFromOptions(J); - CHKERRQ(ierr); - ierr = MatLoad(J, fd); - CHKERRQ(ierr); - ierr = PetscViewerDestroy(&fd); - CHKERRQ(ierr); - } else { // create Jacobian matrix - - /* number of degrees (variables) at each grid point */ - PetscInt dof = n3Dvars(); - - // Maximum allowable size of stencil in x is the number of guard cells - PetscInt stencil_width_estimate = options->operator[]("stencil_width_estimate") - .withDefault(bout::globals::mesh->xstart); - // This is the stencil in each direction (*2) along each dimension - // (*3), plus the point itself. Not sure if this is correct - // though, on several levels: - // 1. Ignores corner points used in e.g. brackets - // 2. Could have different stencil widths in each dimension - // 3. FFTs couple every single point together - PetscInt cols = stencil_width_estimate * 2 * 3 + 1; - PetscInt prealloc; // = cols*dof; - - ierr = MatCreate(comm, &J); - CHKERRQ(ierr); - ierr = MatSetType(J, MATBAIJ); - CHKERRQ(ierr); - ierr = MatSetSizes(J, local_N, local_N, neq, neq); - CHKERRQ(ierr); - ierr = MatSetFromOptions(J); - CHKERRQ(ierr); + if (user_precon) { + output_info.write("\tUsing user-supplied preconditioner\n"); + if (!hasPreconditioner()) { + throw BoutException("Model does not define a preconditioner"); + } + // Set a Shell (matrix-free) preconditioner type + PCSetType(pc, PCSHELL); - // Get nonzero pattern of J - color_none !!! - prealloc = cols * dof * dof; - ierr = MatSeqAIJSetPreallocation(J, prealloc, nullptr); - CHKERRQ(ierr); - ierr = MatMPIAIJSetPreallocation(J, prealloc, nullptr, prealloc, nullptr); - CHKERRQ(ierr); + // Specify the preconditioner function + PCShellSetApply(pc, snesPCapply); - // why nonzeros=295900, allocated nonzeros=2816000/12800000 (*dof*dof), number of mallocs used during MatSetValues calls =256? - prealloc = cols; - ierr = MatSeqBAIJSetPreallocation(J, dof, prealloc, nullptr); - CHKERRQ(ierr); - ierr = MatMPIBAIJSetPreallocation(J, dof, prealloc, nullptr, prealloc, nullptr); - CHKERRQ(ierr); -#if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsHasName(nullptr, nullptr, "-J_slowfd", &J_slowfd); - CHKERRQ(ierr); -#else - ierr = PetscOptionsHasName(nullptr, "-J_slowfd", &J_slowfd); - CHKERRQ(ierr); -#endif - if (J_slowfd != 0U) { // create Jacobian matrix by slow fd - ierr = SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, nullptr); - CHKERRQ(ierr); - output_info << "SNESComputeJacobian J by slow fd...\n"; + // Context used to supply object pointer + PCShellSetContext(pc, this); -#if PETSC_VERSION_GE(3, 5, 0) - ierr = TSComputeRHSJacobian(ts, simtime, u, J, J); + // Set name of preconditioner + ierr = PCShellSetName(pc, "PhysicsPreconditioner"); CHKERRQ(ierr); -#else - MatStructure flg; - ierr = TSComputeRHSJacobian(ts, simtime, u, &J, &J, &flg); - CHKERRQ(ierr); -#endif - output_info << "compute J by slow fd is done.\n"; - } else { // get sparse pattern of the Jacobian - throw BoutException("Sorry, unimplemented way of setting PETSc preconditioner. " - "Either set a preconditioner function with 'setPrecon' " - "in your code, or load a matrix with '-Jload' on the " - "command line, or calculate with finite differences " - "with '-J_slowfd' (on command line)"); + } else { + // Can't use preconditioner because no Jacobian matrix available + PCSetType(pc, PCNONE); } - } - - // Create coloring context of J to be used during time stepping - output_info << " Create coloring ...\n"; + // Callback to get the shift parameter + // Use matrix free for both operator and preconditioner + TSSetIJacobian(ts, Jmf, Jmf, solver_ijacobian, this); - ISColoring iscoloring; -#if PETSC_VERSION_GE(3, 5, 0) - MatColoring coloring; - MatColoringCreate(Jmf, &coloring); - MatColoringSetType(coloring, MATCOLORINGSL); - MatColoringSetFromOptions(coloring); - // Calculate index sets - MatColoringApply(coloring, &iscoloring); - MatColoringDestroy(&coloring); -#else - ierr = MatGetColoring(J, MATCOLORINGSL, &iscoloring); - CHKERRQ(ierr); -#endif - ierr = MatFDColoringCreate(J, iscoloring, &matfdcoloring); - CHKERRQ(ierr); - - ierr = MatFDColoringSetFromOptions(matfdcoloring); - CHKERRQ(ierr); - ierr = ISColoringDestroy(&iscoloring); - CHKERRQ(ierr); - throw BoutException("Coloring is not working"); - ierr = SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring); - CHKERRQ(ierr); - - // Write J in binary for study - see ~petsc/src/mat/examples/tests/ex124.c -#if PETSC_VERSION_GE(3, 7, 0) - ierr = PetscOptionsHasName(nullptr, nullptr, "-J_write", &J_write); - CHKERRQ(ierr); -#else - ierr = PetscOptionsHasName(nullptr, "-J_write", &J_write); - CHKERRQ(ierr); -#endif - if (J_write != 0U) { - PetscViewer viewer = nullptr; - output_info.write("\n[{:d}] Test TSComputeRHSJacobian() ...\n", rank); -#if PETSC_VERSION_GE(3, 5, 0) - ierr = TSComputeRHSJacobian(ts, simtime, u, J, J); - CHKERRQ(ierr); + } else { + // Set PC type from input + if (pc_type != "default") { + PCSetType(pc, pc_type.c_str()); + + if (pc_type == "hypre") { +#if PETSC_HAVE_HYPRE + // Set the type of hypre preconditioner + PCHYPRESetType(pc, pc_hypre_type.c_str()); #else - MatStructure J_structure; - ierr = TSComputeRHSJacobian(ts, simtime, u, &J, &J, &J_structure); - CHKERRQ(ierr); + throw BoutException("PETSc was not configured with Hypre."); #endif + } + } - output.write("[{:d}] TSComputeRHSJacobian is done\n", rank); + // Calculate a Jacobian matrix using finite differences. The finite + // difference Jacobian (Jfd) may be used for both operator and + // preconditioner or, if matrix_free_operator, in only the + // preconditioner. + + if (use_coloring) { + // Use matrix coloring. This greatly reduces the number of + // times the rhs() function needs to be evaluated when + // calculating the Jacobian, by identifying which quantities may + // be simultaneously perturbed. + + // Use global mesh for now + Mesh* mesh = bout::globals::mesh; + + ////////////////////////////////////////////////// + // Get the local indices by starting at 0 + Field3D index = globalIndex(0); + + ////////////////////////////////////////////////// + // Pre-allocate PETSc storage + + output_progress.write("Setting Jacobian matrix sizes\n"); + + const int n2d = f2d.size(); + const int n3d = f3d.size(); + + // Set size of Matrix on each processor to nlocal x nlocal + MatCreate(BoutComm::get(), &Jfd); + MatSetOption(Jfd, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); + MatSetSizes(Jfd, nlocal, nlocal, PETSC_DETERMINE, PETSC_DETERMINE); + MatSetFromOptions(Jfd); + // Determine which row/columns of the matrix are locally owned + int Istart, Iend; + MatGetOwnershipRange(Jfd, &Istart, &Iend); + // Convert local into global indices + // Note: Not in the boundary cells, to keep -1 values + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + index[i] += Istart; + } + // Now communicate to fill guard cells + mesh->communicate(index); + + // Non-zero elements on this processor + std::vector d_nnz; + std::vector o_nnz; + auto n_square = (*options)["stencil:square"] + .doc("Extent of stencil (square)") + .withDefault(0); + auto n_taxi = (*options)["stencil:taxi"] + .doc("Extent of stencil (taxi-cab norm)") + .withDefault(0); + auto n_cross = (*options)["stencil:cross"] + .doc("Extent of stencil (cross)") + .withDefault(0); + // Set n_taxi 2 if nothing else is set + // Probably a better way to do this + if (n_square == 0 && n_taxi == 0 && n_cross == 0) { + output_info.write("Setting solver:stencil:taxi = 2\n"); + n_taxi = 2; + } + + auto const xy_offsets = ColoringStencil::getOffsets(n_square, n_taxi, n_cross); + { + // This is ugly but can't think of a better and robust way to + // count the non-zeros for some arbitrary stencil + // effectively the same loop as the one that sets the non-zeros below + std::vector> d_nnz_map2d(nlocal); + std::vector> o_nnz_map2d(nlocal); + std::vector> d_nnz_map3d(nlocal); + std::vector> o_nnz_map3d(nlocal); + // Loop over every element in 2D to count the *unique* non-zeros + for (int x = mesh->xstart; x <= mesh->xend; x++) { + for (int y = mesh->ystart; y <= mesh->yend; y++) { + + const int ind0 = ROUND(index(x, y, 0)) - Istart; + + // 2D fields + for (int i = 0; i < n2d; i++) { + const PetscInt row = ind0 + i; + // Loop through each point in the stencil + for (const auto& [x_off, y_off] : xy_offsets) { + const int xi = x + x_off; + const int yi = y + y_off; + if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) + || (yi >= mesh->LocalNy)) { + continue; + } + + const int ind2 = ROUND(index(xi, yi, 0)); + if (ind2 < 0) { + continue; // A boundary point + } + + // Depends on all variables on this cell + for (int j = 0; j < n2d; j++) { + const PetscInt col = ind2 + j; + if (col >= Istart && col < Iend) { + d_nnz_map2d[row].insert(col); + } else { + o_nnz_map2d[row].insert(col); + } + } + } + } + // 3D fields + for (int z = 0; z < mesh->LocalNz; z++) { + const int ind = ROUND(index(x, y, z)) - Istart; + + for (int i = 0; i < n3d; i++) { + PetscInt row = ind + i; + if (z == 0) { + row += n2d; + } + + // Depends on 2D fields + for (int j = 0; j < n2d; j++) { + const PetscInt col = ind0 + j; + if (col >= Istart && col < Iend) { + d_nnz_map2d[row].insert(col); + } else { + o_nnz_map2d[row].insert(col); + } + } + + // Star pattern + for (const auto& [x_off, y_off] : xy_offsets) { + const int xi = x + x_off; + const int yi = y + y_off; + + if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) + || (yi >= mesh->LocalNy)) { + continue; + } + + int ind2 = ROUND(index(xi, yi, 0)); + if (ind2 < 0) { + continue; // Boundary point + } + + if (z == 0) { + ind2 += n2d; + } + + // 3D fields on this cell + for (int j = 0; j < n3d; j++) { + const PetscInt col = ind2 + j; + if (col >= Istart && col < Iend) { + d_nnz_map3d[row].insert(col); + } else { + o_nnz_map3d[row].insert(col); + } + } + } + } + } + } + } + + d_nnz.reserve(nlocal); + d_nnz.reserve(nlocal); + + for (int i = 0; i < nlocal; ++i) { + // Assume all elements in the z direction are potentially coupled + d_nnz.emplace_back(d_nnz_map3d[i].size() * mesh->LocalNz + + d_nnz_map2d[i].size()); + o_nnz.emplace_back(o_nnz_map3d[i].size() * mesh->LocalNz + + o_nnz_map2d[i].size()); + } + } + + output_progress.write("Pre-allocating Jacobian\n"); + // Pre-allocate + MatMPIAIJSetPreallocation(Jfd, 0, d_nnz.data(), 0, o_nnz.data()); + MatSeqAIJSetPreallocation(Jfd, 0, d_nnz.data()); + MatSetUp(Jfd); + MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE); + + ////////////////////////////////////////////////// + // Mark non-zero entries + + output_progress.write("Marking non-zero Jacobian entries\n"); + PetscScalar val = 1.0; + for (int x = mesh->xstart; x <= mesh->xend; x++) { + for (int y = mesh->ystart; y <= mesh->yend; y++) { + + const int ind0 = ROUND(index(x, y, 0)); + + // 2D fields + for (int i = 0; i < n2d; i++) { + const PetscInt row = ind0 + i; + + // Loop through each point in the stencil + for (const auto& [x_off, y_off] : xy_offsets) { + const int xi = x + x_off; + const int yi = y + y_off; + if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) + || (yi >= mesh->LocalNy)) { + continue; + } + + int ind2 = ROUND(index(xi, yi, 0)); + if (ind2 < 0) { + continue; // A boundary point + } + + // Depends on all variables on this cell + for (int j = 0; j < n2d; j++) { + PetscInt col = ind2 + j; + ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); + CHKERRQ(ierr); + } + } + } + // 3D fields + for (int z = 0; z < mesh->LocalNz; z++) { + int ind = ROUND(index(x, y, z)); + + for (int i = 0; i < n3d; i++) { + PetscInt row = ind + i; + if (z == 0) { + row += n2d; + } + + // Depends on 2D fields + for (int j = 0; j < n2d; j++) { + PetscInt col = ind0 + j; + ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); + CHKERRQ(ierr); + } + + // Star pattern + for (const auto& [x_off, y_off] : xy_offsets) { + int xi = x + x_off; + int yi = y + y_off; + + if ((xi < 0) || (yi < 0) || (xi >= mesh->LocalNx) + || (yi >= mesh->LocalNy)) { + continue; + } + for (int zi = 0; zi < mesh->LocalNz; ++zi) { + int ind2 = ROUND(index(xi, yi, zi)); + if (ind2 < 0) { + continue; // Boundary point + } + + if (z == 0) { + ind2 += n2d; + } + + // 3D fields on this cell + for (int j = 0; j < n3d; j++) { + PetscInt col = ind2 + j; + ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); + + if (ierr != 0) { + output.write("ERROR: {} {} : ({}, {}) -> ({}, {}) : {} -> {}\n", + row, x, y, xi, yi, ind2, ind2 + n3d - 1); + } + CHKERRQ(ierr); + } + } + } + } + } + } + } + + // Finished marking non-zero entries + + output_progress.write("Assembling Jacobian matrix\n"); + + // Assemble Matrix + MatAssemblyBegin(Jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Jfd, MAT_FINAL_ASSEMBLY); + + { + // Test if the matrix is symmetric + // Values are 0 or 1 so tolerance (1e-5) shouldn't matter + PetscBool symmetric; + ierr = MatIsSymmetric(Jfd, 1e-5, &symmetric); + CHKERRQ(ierr); + if (!symmetric) { + output_warn.write("Jacobian pattern is not symmetric\n"); + } + } + + // The above can miss entries around the X-point branch cut: + // The diagonal terms are complicated because moving in X then Y + // is different from moving in Y then X at the X-point. + // Making sure the colouring matrix is symmetric does not + // necessarily give the correct stencil but may help. + if ((*options)["force_symmetric_coloring"] + .doc("Modifies coloring matrix to force it to be symmetric") + .withDefault(false)) { + Mat Jfd_T; + MatCreateTranspose(Jfd, &Jfd_T); + MatAXPY(Jfd, 1, Jfd_T, DIFFERENT_NONZERO_PATTERN); + } + + output_progress.write("Creating Jacobian coloring\n"); + updateColoring(); - if (J_slowfd != 0U) { - output_info.write("[{:d}] writing J in binary to data/Jrhs_dense.dat...\n", rank); - ierr = PetscViewerBinaryOpen(comm, "data/Jrhs_dense.dat", FILE_MODE_WRITE, &viewer); - CHKERRQ(ierr); } else { - output_info.write("[{:d}] writing J in binary to data/Jrhs_sparse.dat...\n", rank); - ierr = - PetscViewerBinaryOpen(comm, "data/Jrhs_sparse.dat", FILE_MODE_WRITE, &viewer); - CHKERRQ(ierr); + // Brute force calculation + // There is usually no reason to use this, except as a check of + // the coloring calculation. + + MatCreateAIJ( + BoutComm::get(), nlocal, nlocal, // Local sizes + PETSC_DETERMINE, PETSC_DETERMINE, // Global sizes + 3, // Number of nonzero entries in diagonal portion of local submatrix + nullptr, + 0, // Number of nonzeros per row in off-diagonal portion of local submatrix + nullptr, &Jfd); + + if (matrix_free_operator) { + SNESSetJacobian(snes, Jmf, Jfd, SNESComputeJacobianDefault, this); + } else { + SNESSetJacobian(snes, Jfd, Jfd, SNESComputeJacobianDefault, this); + } + + MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); } - ierr = MatView(J, viewer); - CHKERRQ(ierr); - ierr = PetscViewerDestroy(&viewer); - CHKERRQ(ierr); } - ierr = PetscLogEventEnd(init_event, 0, 0, 0, 0); - CHKERRQ(ierr); - - PetscFunctionReturn(0); + // Get runtime options + lib.setOptionsFromInputFile(ts); + + { + // Some reporting + TSType tstype; + TSGetType(ts, &tstype); + output_info.write("TS Type : {}\n", tstype); + TSAdaptType adapttype; + TSAdaptGetType(adapt, &adapttype); + output_info.write("TS Adapt Type : {}\n", adapttype); + SNESType snestype; + SNESGetType(snes, &snestype); + output_info.write("SNES Type : {}\n", snestype); + KSPType ksptype; + KSPGetType(ksp, &ksptype); + if (ksptype) { + output_info.write("KSP Type : {}\n", ksptype); + } + PCType pctype; + PCGetType(pc, &pctype); + if (pctype) { + output_info.write("PC Type : {}\n", pctype); + } + } + return 0; } -/************************************************************************** - * Run - Advance time - **************************************************************************/ - +// Starts the time integrator PetscErrorCode PetscSolver::run() { - // Set when the next call to monitor is desired next_output = simtime + getOutputTimestep(); PetscFunctionBegin; - if (this->output_flag) { - prev_linear_its = 0; - bout_snes_time = bout::globals::mpi->MPI_Wtime(); - } - + // Run the PETSc time integrator + // The PetscMonitor function is responsible for regular outputs CHKERRQ(TSSolve(ts, u)); - // Gawd, everything is a hack - if ((this->output_flag != 0U) and (BoutComm::rank() == 0)) { - Output petsc_info(output_name); - // Don't write to stdout - petsc_info.disable(); - petsc_info.write("SNES Iteration, KSP Iterations, Wall Time, Norm\n"); - for (const auto& info : snes_list) { - petsc_info.write("{:d}, {:d}, {:e}, {:e}\n", info.it, info.linear_its, info.time, - info.norm); - } - } - PetscFunctionReturn(0); } -/************************************************************************** - * RHS function - **************************************************************************/ - -PetscErrorCode PetscSolver::rhs(TS UNUSED(ts), BoutReal t, Vec udata, Vec dudata) { +// Evaluate the user-supplied RHS function +// This method is called from the static PETSc callback functions +PetscErrorCode PetscSolver::rhs(BoutReal t, Vec udata, Vec dudata, bool linear) { TRACE("Running RHS: PetscSolver::rhs({:e})", t); - const BoutReal* udata_array; - BoutReal* dudata_array; + simtime = t; // This will be used in rhs_differencing PetscFunctionBegin; // Load state from PETSc + const BoutReal* udata_array; VecGetArrayRead(udata, &udata_array); load_vars(const_cast(udata_array)); VecRestoreArrayRead(udata, &udata_array); // Call RHS function - run_rhs(t); + run_rhs(t, linear); // Save derivatives to PETSc + BoutReal* dudata_array; VecGetArray(dudata, &dudata_array); save_derivs(dudata_array); VecRestoreArray(dudata, &dudata_array); @@ -623,19 +914,17 @@ PetscErrorCode PetscSolver::rhs(TS UNUSED(ts), BoutReal t, Vec udata, Vec dudata PetscFunctionReturn(0); } -/************************************************************************** - * Preconditioner function - **************************************************************************/ +// This is used for Jacobian-vector product and Jacobian coloring +PetscErrorCode PetscSolver::rhs_differencing(Vec globalin, Vec globalout) { + return rhs(simtime, globalin, globalout, true); +} -PetscErrorCode PetscSolver::pre(PC UNUSED(pc), Vec x, Vec y) { +// Matrix-free preconditioner function +PetscErrorCode PetscSolver::pre(Vec x, Vec y) { TRACE("PetscSolver::pre()"); BoutReal* data; - if (diagnose) { - output << "Preconditioning" << endl; - } - // Load state VecGetArray(state, &data); load_vars(data); @@ -654,355 +943,44 @@ PetscErrorCode PetscSolver::pre(PC UNUSED(pc), Vec x, Vec y) { save_derivs(data); VecRestoreArray(y, &data); - // Petsc's definition of Jacobian differs by a factor from Sundials' + // Petsc's definition of Jacobian differs by a factor from SUNDIALS' + // PETSc solves (scale + J)^-1 + // SUNDIALS solves (I + gamma J)^-1 PetscErrorCode ierr = VecScale(y, shift); CHKERRQ(ierr); return 0; } -/************************************************************************** - * User-supplied Jacobian function J(state) * x = y - **************************************************************************/ - -PetscErrorCode PetscSolver::jac(Vec x, Vec y) { - TRACE("PetscSolver::jac()"); - - BoutReal* data; - - if (diagnose) { - output << "Jacobian evaluation\n"; - } - - // Load state - VecGetArray(state, &data); - load_vars(data); - VecRestoreArray(state, &data); - - // Load vector to be operated on into F_vars - VecGetArray(x, &data); - load_derivs(data); - VecRestoreArray(x, &data); - - // Call the Jacobian function - runJacobian(ts_time); - - // Save the solution from time derivatives - VecGetArray(y, &data); - save_derivs(data); - VecRestoreArray(y, &data); - - // y = a * x - y - PetscErrorCode ierr = VecAXPBY(y, shift, -1.0, x); - CHKERRQ(ierr); - - return 0; -} - -/************************************************************************** - * Static functions which can be used for PETSc callbacks - **************************************************************************/ - -PetscErrorCode solver_f(TS ts, BoutReal t, Vec globalin, Vec globalout, void* f_data) { - PetscFunctionBegin; - auto* s = static_cast(f_data); - PetscLogEventBegin(s->solver_event, 0, 0, 0, 0); - s->rhs(ts, t, globalin, globalout); - PetscLogEventEnd(s->solver_event, 0, 0, 0, 0); - PetscFunctionReturn(0); -} - -/* - FormIFunction = Udot - RHSFunction -*/ -PetscErrorCode solver_if(TS ts, BoutReal t, Vec globalin, Vec globalindot, Vec globalout, - void* f_data) { - PetscErrorCode ierr; - - PetscFunctionBegin; - ierr = solver_f(ts, t, globalin, globalout, f_data); - CHKERRQ(ierr); - - ierr = VecAYPX(globalout, -1.0, globalindot); - CHKERRQ(ierr); // globalout = globalindot + (-1)globalout - PetscFunctionReturn(0); -} - -#if PETSC_VERSION_GE(3, 5, 0) -PetscErrorCode solver_rhsjacobian(TS UNUSED(ts), BoutReal UNUSED(t), Vec UNUSED(globalin), - Mat J, Mat Jpre, void* UNUSED(f_data)) { - PetscErrorCode ierr; - - PetscFunctionBegin; - ierr = MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - ierr = MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - if (J != Jpre) { - ierr = MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - ierr = MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - } - PetscFunctionReturn(0); -} -#else -PetscErrorCode solver_rhsjacobian([[maybe_unused]] TS ts, [[maybe_unused]] BoutReal t, - [[maybe_unused]] Vec globalin, Mat* J, Mat* Jpre, - [[maybe_unused]] MatStructure* str, - [[maybe_unused]] void* f_data) { - PetscErrorCode ierr; - - PetscFunctionBegin; - ierr = MatAssemblyBegin(*Jpre, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - ierr = MatAssemblyEnd(*Jpre, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - if (*J != *Jpre) { - ierr = MatAssemblyBegin(*J, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - ierr = MatAssemblyEnd(*J, MAT_FINAL_ASSEMBLY); - CHKERRQ(ierr); - } - PetscFunctionReturn(0); -} -#endif - -/* - solver_ijacobian - Compute IJacobian = dF/dU + a dF/dUdot - a dummy matrix used for pc=none -*/ -#if PETSC_VERSION_GE(3, 5, 0) -PetscErrorCode solver_ijacobian(TS ts, BoutReal t, Vec globalin, Vec UNUSED(globalindot), - PetscReal a, Mat J, Mat Jpre, void* f_data) { - PetscErrorCode ierr; - - PetscFunctionBegin; - ierr = solver_rhsjacobian(ts, t, globalin, J, Jpre, f_data); - CHKERRQ(ierr); - - ////// Save data for preconditioner - auto* solver = static_cast(f_data); - - if (solver->diagnose) { - output << "Saving state, t = " << t << ", a = " << a << endl; - } - - solver->shift = a; // Save the shift 'a' - solver->state = globalin; // Save system state - solver->ts_time = t; - - PetscFunctionReturn(0); -} -#else -PetscErrorCode solver_ijacobian(TS ts, BoutReal t, Vec globalin, - [[maybe_unused]] Vec globalindot, - [[maybe_unused]] PetscReal a, Mat* J, Mat* Jpre, - MatStructure* str, void* f_data) { - PetscErrorCode ierr; - - PetscFunctionBegin; - ierr = solver_rhsjacobian(ts, t, globalin, J, Jpre, str, (void*)f_data); - CHKERRQ(ierr); - - ////// Save data for preconditioner - PetscSolver* solver = (PetscSolver*)f_data; - - if (solver->diagnose) { - output << "Saving state, t = " << t << ", a = " << a << endl; - } - - solver->shift = a; // Save the shift 'a' - solver->state = globalin; // Save system state - solver->ts_time = t; - - PetscFunctionReturn(0); -} -#endif - -/* - solver_ijacobianfd - Compute IJacobian = dF/dU + a dF/dUdot using finite deference - not implemented yet -*/ -#if PETSC_VERSION_GE(3, 5, 0) -PetscErrorCode solver_ijacobianfd(TS ts, BoutReal t, Vec globalin, - Vec UNUSED(globalindot), PetscReal UNUSED(a), Mat J, - Mat Jpre, void* f_data) { - PetscErrorCode ierr; - - PetscFunctionBegin; - ierr = solver_rhsjacobian(ts, t, globalin, J, Jpre, f_data); - CHKERRQ(ierr); - //*Jpre + a - PetscFunctionReturn(0); -} -#else -PetscErrorCode solver_ijacobianfd(TS ts, BoutReal t, Vec globalin, Vec globalindot, - PetscReal a, Mat* J, Mat* Jpre, MatStructure* str, - void* f_data) { - PetscErrorCode ierr; - - PetscFunctionBegin; - ierr = solver_rhsjacobian(ts, t, globalin, J, Jpre, str, (void*)f_data); - CHKERRQ(ierr); - //*Jpre + a - PetscFunctionReturn(0); -} -#endif -//----------------------------------------- - -PetscErrorCode PhysicsSNESApply(SNES snes, Vec x) { - PetscErrorCode ierr; - Vec F, Fout; - PetscReal fnorm = 0., foutnorm = 0., dot = 0.; - KSP ksp; - PC pc; - Mat A, B; - - PetscFunctionBegin; - ierr = SNESGetJacobian(snes, &A, &B, nullptr, nullptr); - CHKERRQ(ierr); -#if PETSC_VERSION_GE(3, 5, 0) - ierr = SNESComputeJacobian(snes, x, A, B); - CHKERRQ(ierr); -#else - MatStructure diff = DIFFERENT_NONZERO_PATTERN; - ierr = SNESComputeJacobian(snes, x, &A, &B, &diff); - CHKERRQ(ierr); -#endif - ierr = SNESGetKSP(snes, &ksp); - CHKERRQ(ierr); - ierr = KSPGetPC(ksp, &pc); - CHKERRQ(ierr); - ierr = SNESGetFunction(snes, &F, nullptr, nullptr); - CHKERRQ(ierr); - ierr = SNESComputeFunction(snes, x, F); - CHKERRQ(ierr); - ierr = SNESGetSolutionUpdate(snes, &Fout); - CHKERRQ(ierr); - - ierr = PCApply(pc, F, Fout); - CHKERRQ(ierr); - ierr = VecNorm(Fout, NORM_2, &foutnorm); - CHKERRQ(ierr); - ierr = VecAXPY(x, -1., Fout); - CHKERRQ(ierr); - ierr = SNESComputeFunction(snes, x, F); - CHKERRQ(ierr); - ierr = VecNorm(F, NORM_2, &fnorm); - CHKERRQ(ierr); - ierr = VecDot(F, Fout, &dot); - CHKERRQ(ierr); - output_info << " (Debug) function norm: " << fnorm << ", P(f) norm " << foutnorm - << ", F \\cdot Fout " << dot << " "; -#if PETSC_VERSION_GE(3, 5, 0) - Vec func; - ierr = SNESGetFunction(snes, &func, nullptr, nullptr); - CHKERRQ(ierr); - ierr = VecNorm(func, NORM_2, &fnorm); - CHKERRQ(ierr); -#else - ierr = SNESSetFunctionNorm(snes, fnorm); - CHKERRQ(ierr); -#endif - ierr = SNESMonitor(snes, 0, fnorm); - CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -PetscErrorCode PhysicsPCApply(PC pc, Vec x, Vec y) { - int ierr; - - // Get the context - PetscSolver* s; - ierr = PCShellGetContext(pc, reinterpret_cast(&s)); - CHKERRQ(ierr); - - PetscFunctionReturn(s->pre(pc, x, y)); -} - -PetscErrorCode PhysicsJacobianApply(Mat J, Vec x, Vec y) { - // Get the context - PetscSolver* s; - int ierr = MatShellGetContext(J, reinterpret_cast(&s)); - CHKERRQ(ierr); - PetscFunctionReturn(s->jac(x, y)); -} - -PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, void* ctx) { - PetscErrorCode ierr; - auto* s = static_cast(ctx); - PetscReal tfinal, dt; - Vec interpolatedX; - const PetscScalar* x; - static int i = 0; - - PetscFunctionBegin; - ierr = TSGetTimeStep(ts, &dt); - CHKERRQ(ierr); - -#if PETSC_VERSION_GE(3, 8, 0) - ierr = TSGetMaxTime(ts, &tfinal); - CHKERRQ(ierr); -#else - ierr = TSGetDuration(ts, nullptr, &tfinal); - CHKERRQ(ierr); -#endif - - /* Duplicate the solution vector X into a work vector */ - ierr = VecDuplicate(X, &interpolatedX); - CHKERRQ(ierr); - while (s->next_output <= t && s->next_output <= tfinal) { - if (s->interpolate) { - ierr = TSInterpolate(ts, s->next_output, interpolatedX); - CHKERRQ(ierr); - } - - /* Place the interpolated values into the global variables */ - ierr = VecGetArrayRead(interpolatedX, &x); - CHKERRQ(ierr); - s->load_vars(const_cast(x)); - ierr = VecRestoreArrayRead(interpolatedX, &x); - CHKERRQ(ierr); - - if (s->call_monitors(simtime, i++, s->getNumberOutputSteps())) { - PetscFunctionReturn(1); - } - - s->next_output += s->getOutputTimestep(); - simtime = s->next_output; - } - - /* Done with vector, so destroy it */ - ierr = VecDestroy(&interpolatedX); - CHKERRQ(ierr); - - PetscFunctionReturn(0); -} - -PetscErrorCode PetscSNESMonitor(SNES snes, PetscInt its, PetscReal norm, void* ctx) { - PetscErrorCode ierr; - PetscInt linear_its = 0; - BoutReal tmp = .0; - snes_info row; - auto* s = static_cast(ctx); +void PetscSolver::updateColoring() { + // Re-calculate the coloring + MatColoring coloring = NULL; + MatColoringCreate(Jfd, &coloring); + MatColoringSetType(coloring, MATCOLORINGGREEDY); + MatColoringSetFromOptions(coloring); - PetscFunctionBegin; + // Calculate new index sets + ISColoring iscoloring = NULL; + MatColoringApply(coloring, &iscoloring); + MatColoringDestroy(&coloring); - if (!its) { - s->prev_linear_its = 0; + // Replace the old coloring with the new one + MatFDColoringDestroy(&fdcoloring); + MatFDColoringCreate(Jfd, iscoloring, &fdcoloring); + MatFDColoringSetFunction(fdcoloring, bout::cast_MatFDColoringFn(solver_rhs_coloring), + this); + MatFDColoringSetFromOptions(fdcoloring); + MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); + ISColoringDestroy(&iscoloring); + + // Replace the CTX pointer in SNES Jacobian + // Note: TSComputeIJacobianDefaultColor + if (matrix_free_operator) { + // Use matrix-free calculation for operator, finite difference for preconditioner + SNESSetJacobian(snes, Jmf, Jfd, SNESComputeJacobianDefaultColor, fdcoloring); + } else { + SNESSetJacobian(snes, Jfd, Jfd, SNESComputeJacobianDefaultColor, fdcoloring); } - ierr = SNESGetLinearSolveIterations(snes, &linear_its); - CHKERRQ(ierr); - tmp = bout::globals::mpi->MPI_Wtime(); - - row.it = its; - s->prev_linear_its = row.linear_its = linear_its - s->prev_linear_its; - row.time = tmp - s->bout_snes_time; - row.norm = norm; - - s->snes_list.push_back(row); - - PetscFunctionReturn(0); } #endif diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index 686a7c6242..71a09ab61d 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -1,11 +1,10 @@ /************************************************************************** * Interface to PETSc solver - * NOTE: This class needs tidying, generalising to use FieldData interface * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -55,35 +54,6 @@ namespace { RegisterSolver registersolverpetsc("petsc"); } -using BoutReal = PetscScalar; -#define OPT_SIZE 40 - -using rhsfunc = int (*)(BoutReal); - -extern BoutReal simtime; - -/// Monitor function called on every internal timestep -extern PetscErrorCode PetscMonitor(TS, PetscInt, PetscReal, Vec, void* ctx); -/// Monitor function for SNES -extern PetscErrorCode PetscSNESMonitor(SNES, PetscInt, PetscReal, void* ctx); - -/// Compute IJacobian = dF/dU + a dF/dUdot - a dummy matrix used for pc=none -#if PETSC_VERSION_GE(3, 5, 0) -extern PetscErrorCode solver_ijacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, - void*); -#else -extern PetscErrorCode solver_ijacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat*, Mat*, - MatStructure*, void*); -#endif - -/// Data for SNES -struct snes_info { - PetscInt it; - PetscInt linear_its; - PetscReal time; - PetscReal norm; -}; - class PetscSolver : public Solver { public: PetscSolver(Options* opts = nullptr); @@ -95,24 +65,17 @@ public: // These functions used internally (but need to be public) /// Wrapper for the RHS function - PetscErrorCode rhs(TS ts, PetscReal t, Vec globalin, Vec globalout); + PetscErrorCode rhs(PetscReal t, Vec globalin, Vec globalout, bool linear); + /// RHS function for differencing + PetscErrorCode rhs_differencing(Vec globalin, Vec globalout); /// Wrapper for the preconditioner - PetscErrorCode pre(PC pc, Vec x, Vec y); - /// Wrapper for the Jacobian function - PetscErrorCode jac(Vec x, Vec y); + PetscErrorCode pre(Vec x, Vec y); // Call back functions that need to access internal state friend PetscErrorCode PetscMonitor(TS, PetscInt, PetscReal, Vec, void* ctx); - friend PetscErrorCode PetscSNESMonitor(SNES, PetscInt, PetscReal, void* ctx); -#if PETSC_VERSION_GE(3, 5, 0) - friend PetscErrorCode solver_ijacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, - void*); -#else - friend PetscErrorCode solver_ijacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat*, Mat*, - MatStructure*, void*); -#endif - PetscLogEvent solver_event, loop_event, init_event; + friend PetscErrorCode solver_ijacobian(TS, BoutReal, Vec, Vec, PetscReal shift, Mat J, + Mat Jpre, void* ctx); private: BoutReal shift; ///< Shift (alpha) parameter from TS @@ -121,28 +84,45 @@ private: PetscLib lib; ///< Handles initialising, finalising PETSc - Vec u{nullptr}; ///< PETSc solution vector - TS ts{nullptr}; ///< PETSc timestepper object - Mat J{nullptr}; ///< RHS Jacobian - Mat Jmf{nullptr}; - MatFDColoring matfdcoloring{nullptr}; - - bool diagnose; ///< If true, print some information about current stage + 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 + Mat Jmf{nullptr}; ///< Matrix Free Jacobian + Mat Jfd{nullptr}; ///< Finite Difference Jacobian + MatFDColoring fdcoloring{nullptr}; ///< Matrix coloring context BoutReal next_output; ///< When the monitor should be called next PetscBool interpolate{PETSC_TRUE}; ///< Whether to interpolate or not - char output_name[PETSC_MAX_PATH_LEN]; - PetscBool output_flag{PETSC_FALSE}; - PetscInt prev_linear_its; - BoutReal bout_snes_time{0.0}; - std::vector snes_list; + bool diagnose; ///< If true, print some information about current stage + bool user_precon; ///< Use user-supplied preconditioning function? + + BoutReal atol; ///< Absolute tolerance + BoutReal rtol; ///< Relative tolerance + BoutReal stol; ///< Convergence tolerance + + int maxnl; ///< Maximum nonlinear iterations per SNES solve + int maxf; ///< Maximum number of function evaluations allowed in the solver (default: 10000) + int maxl; ///< Maximum linear iterations + + std::string ts_type; ///< PETSc TS time solver type + std::string adapt_type; ///< TSAdaptType timestep adaptation + std::string snes_type; ///< PETSc SNES nonlinear solver type + std::string ksp_type; ///< PETSc KSP linear solver type + std::string pc_type; ///< Preconditioner type + std::string pc_hypre_type; ///< Hypre preconditioner type + std::string line_search_type; ///< Line search type + + bool matrix_free; ///< Use matrix free Jacobian + bool matrix_free_operator; ///< Use matrix free Jacobian in the operator? + int lag_jacobian; ///< Re-use Jacobian + bool use_coloring; ///< Use matrix coloring + void updateColoring(); ///< Updates the coloring using Jfd + + bool kspsetinitialguessnonzero; ///< Set initial guess to non-zero - bool adaptive; ///< Use adaptive timestepping - bool use_precon, use_jacobian; - BoutReal abstol, reltol; - bool adams_moulton; BoutReal start_timestep; int mxstep; }; From 46f7344592dcd9ceab4ac6d94093c95ec7b1aece Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 14 Aug 2025 21:56:16 -0700 Subject: [PATCH 03/18] petsc solver fixes - 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. --- src/solver/impls/petsc/petsc.cxx | 55 +++++++++++++------------------- src/solver/impls/petsc/petsc.hxx | 4 +-- 2 files changed, 23 insertions(+), 36 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 22d8594c28..fce7f1178f 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -104,25 +104,6 @@ static PetscErrorCode solver_rhs(TS UNUSED(ts), BoutReal t, Vec globalin, Vec gl PetscFunctionReturn(0); } -// PETSc callback function for finite differencing. -// May use a linearised form of the RHS function. -static PetscErrorCode solver_rhs_differencing(void* ctx, Vec globalin, Vec globalout) { - PetscFunctionBegin; - auto* s = static_cast(ctx); - s->rhs_differencing(globalin, globalout); - PetscFunctionReturn(0); -} - -// PETSc callback function for coloring. -// May use a linearised form of the RHS function. -static PetscErrorCode solver_rhs_coloring(void* UNUSED(snes), Vec globalin, Vec globalout, - void* ctx) { - PetscFunctionBegin; - auto* s = static_cast(ctx); - s->rhs_differencing(globalin, globalout); - PetscFunctionReturn(0); -} - // 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 @@ -172,9 +153,14 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo // The internal timestepper may have stepped over multiple output times while (s->next_output <= t && s->next_output <= tfinal) { + BoutReal output_time = t; if (s->interpolate) { ierr = TSInterpolate(ts, s->next_output, interpolatedX); - CHKERRQ(ierr); + if (ierr != PETSC_SUCCESS) { + throw BoutException("This PETSc TS does not support interpolation. Use a " + "different method or set solver:interpolate=false"); + } + output_time = s->next_output; } // Place the interpolated values into the global variables @@ -185,11 +171,11 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo ierr = VecRestoreArrayRead(interpolatedX, &x); CHKERRQ(ierr); - if (s->call_monitors(s->next_output, i++, s->getNumberOutputSteps())) { + if (s->call_monitors(output_time, i++, s->getNumberOutputSteps())) { PetscFunctionReturn(1); } - s->next_output += s->getOutputTimestep(); + s->next_output = output_time + s->getOutputTimestep(); } // Done with vector, so destroy it @@ -200,7 +186,9 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo } PetscSolver::PetscSolver(Options* opts) - : Solver(opts), + : Solver(opts), interpolate((*options)["interpolate"] + .doc("Interpolate to regular output times?") + .withDefault(true)), diagnose( (*options)["diagnose"].doc("Enable some diagnostic output").withDefault(false)), user_precon((*options)["user_precon"] @@ -379,7 +367,11 @@ int PetscSolver::init() { // Allow TS to step over the final time // Note: This does not affect intermediate outputs, that // are always interpolated (in PetscMonitor) - ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE); + if (interpolate) { + ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE); + } else { + ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP); + } CHKERRQ(ierr); ierr = TSMonitorSet(ts, PetscMonitor, this, nullptr); @@ -892,7 +884,7 @@ PetscErrorCode PetscSolver::run() { PetscErrorCode PetscSolver::rhs(BoutReal t, Vec udata, Vec dudata, bool linear) { TRACE("Running RHS: PetscSolver::rhs({:e})", t); - simtime = t; // This will be used in rhs_differencing + simtime = t; PetscFunctionBegin; @@ -914,11 +906,6 @@ PetscErrorCode PetscSolver::rhs(BoutReal t, Vec udata, Vec dudata, bool linear) PetscFunctionReturn(0); } -// This is used for Jacobian-vector product and Jacobian coloring -PetscErrorCode PetscSolver::rhs_differencing(Vec globalin, Vec globalout) { - return rhs(simtime, globalin, globalout, true); -} - // Matrix-free preconditioner function PetscErrorCode PetscSolver::pre(Vec x, Vec y) { TRACE("PetscSolver::pre()"); @@ -967,14 +954,16 @@ void PetscSolver::updateColoring() { // Replace the old coloring with the new one MatFDColoringDestroy(&fdcoloring); MatFDColoringCreate(Jfd, iscoloring, &fdcoloring); - MatFDColoringSetFunction(fdcoloring, bout::cast_MatFDColoringFn(solver_rhs_coloring), - this); + // Use the SNES function that is defined by the TS method + // SNESTSFormFunction is defined in PETSc ts.c + // The ctx pointer should be the TS object + MatFDColoringSetFunction(fdcoloring, bout::cast_MatFDColoringFn(SNESTSFormFunction), + ts); MatFDColoringSetFromOptions(fdcoloring); MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); ISColoringDestroy(&iscoloring); // Replace the CTX pointer in SNES Jacobian - // Note: TSComputeIJacobianDefaultColor if (matrix_free_operator) { // Use matrix-free calculation for operator, finite difference for preconditioner SNESSetJacobian(snes, Jmf, Jfd, SNESComputeJacobianDefaultColor, fdcoloring); diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index 71a09ab61d..5108e8cb30 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -66,8 +66,6 @@ public: /// Wrapper for the RHS function PetscErrorCode rhs(PetscReal t, Vec globalin, Vec globalout, bool linear); - /// RHS function for differencing - PetscErrorCode rhs_differencing(Vec globalin, Vec globalout); /// Wrapper for the preconditioner PetscErrorCode pre(Vec x, Vec y); @@ -94,7 +92,7 @@ private: BoutReal next_output; ///< When the monitor should be called next - PetscBool interpolate{PETSC_TRUE}; ///< Whether to interpolate or not + bool interpolate; ///< Interpolate to regular times? bool diagnose; ///< If true, print some information about current stage bool user_precon; ///< Use user-supplied preconditioning function? From e4fe421de08e597ff0db426839d5b394247c1425 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 14 Aug 2025 22:31:31 -0700 Subject: [PATCH 04/18] petsc solver: Recover from SNES failures Uses TSSetMaxSNESFailures, that can be overridden with a petsc option. --- src/solver/impls/petsc/petsc.cxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index fce7f1178f..a6971c7270 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -361,6 +361,9 @@ int PetscSolver::init() { CHKERRQ(ierr); #endif + // Allow TS to recover from SNES failures + PetscCall(TSSetMaxSNESFailures(ts, PETSC_UNLIMITED)); + // Set the current solution ierr = TSSetSolution(ts, u); CHKERRQ(ierr); From 3f67ac9ac22ee0c8b0a13ace14e98720d586d1f0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 14 Aug 2025 22:34:24 -0700 Subject: [PATCH 05/18] petsc: Tidying with PetscCall Replaces CHKERRQ macro call with a wrapper around the PETSc function call. --- src/solver/impls/petsc/petsc.cxx | 59 +++++++++++--------------------- 1 file changed, 20 insertions(+), 39 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index a6971c7270..126aa46532 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -283,48 +283,35 @@ int PetscSolver::init() { // Create a vector to contain the state PetscCall(VecCreate(BoutComm::get(), &u)); - ierr = VecSetSizes(u, nlocal, PETSC_DECIDE); - CHKERRQ(ierr); - ierr = VecSetFromOptions(u); - CHKERRQ(ierr); + PetscCall(VecSetSizes(u, nlocal, PETSC_DECIDE)); + PetscCall(VecSetFromOptions(u)); // Save initial state to PETSc Vec BoutReal* udata; // Pointer to data array in vector u. - ierr = VecGetArray(u, &udata); - CHKERRQ(ierr); + PetscCall(VecGetArray(u, &udata)); save_vars(udata); - ierr = VecRestoreArray(u, &udata); - CHKERRQ(ierr); + PetscCall(VecRestoreArray(u, &udata)); // Create timestepper - ierr = TSCreate(BoutComm::get(), &ts); - CHKERRQ(ierr); - ierr = TSSetProblemType(ts, TS_NONLINEAR); - CHKERRQ(ierr); - ierr = TSSetType(ts, ts_type.c_str()); - CHKERRQ(ierr); - ierr = TSSetApplicationContext(ts, this); - CHKERRQ(ierr); + PetscCall(TSCreate(BoutComm::get(), &ts)); + PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); + PetscCall(TSSetType(ts, ts_type.c_str())); + PetscCall(TSSetApplicationContext(ts, this)); // Set user provided RHSFunction // Need to duplicate the solution vector for the residual Vec rhs_vec; - ierr = VecDuplicate(u, &rhs_vec); - CHKERRQ(ierr); - ierr = TSSetRHSFunction(ts, rhs_vec, solver_rhs, this); - CHKERRQ(ierr); + PetscCall(VecDuplicate(u, &rhs_vec)); + PetscCall(TSSetRHSFunction(ts, rhs_vec, solver_rhs, this)); // Set up adaptive time-stepping TSAdapt adapt; - ierr = TSGetAdapt(ts, &adapt); - CHKERRQ(ierr); - ierr = TSAdaptSetType(adapt, adapt_type.c_str()); - CHKERRQ(ierr); + PetscCall(TSGetAdapt(ts, &adapt)); + PetscCall(TSAdaptSetType(adapt, adapt_type.c_str())); // Set default absolute/relative tolerances // Note: Vector atol and rtol not given - ierr = TSSetTolerances(ts, atol, nullptr, rtol, nullptr); - CHKERRQ(ierr); + PetscCall(TSSetTolerances(ts, atol, nullptr, rtol, nullptr)); if (ts_type == TSSUNDIALS) { #if PETSC_HAS_SUNDIALS // The PETSc interface to SUNDIALS' CVODE @@ -337,12 +324,10 @@ int PetscSolver::init() { } // Initial time of the simulation state - ierr = TSSetTime(ts, simtime); - CHKERRQ(ierr); + PetscCall(TSSetTime(ts, simtime)); next_output = simtime; - ierr = TSSetTimeStep(ts, start_timestep); - CHKERRQ(ierr); + PetscCall(TSSetTimeStep(ts, start_timestep)); // Total number of steps PetscInt total_steps = mxstep * getNumberOutputSteps(); @@ -352,21 +337,17 @@ int PetscSolver::init() { simtime); #if PETSC_VERSION_GE(3, 8, 0) - ierr = TSSetMaxSteps(ts, total_steps); - CHKERRQ(ierr); - ierr = TSSetMaxTime(ts, tfinal); - CHKERRQ(ierr); + PetscCall(TSSetMaxSteps(ts, total_steps)); + PetscCall(TSSetMaxTime(ts, tfinal)); #else - ierr = TSSetDuration(ts, total_steps, tfinal); - CHKERRQ(ierr); + PetscCall(TSSetDuration(ts, total_steps, tfinal)); #endif // Allow TS to recover from SNES failures PetscCall(TSSetMaxSNESFailures(ts, PETSC_UNLIMITED)); // Set the current solution - ierr = TSSetSolution(ts, u); - CHKERRQ(ierr); + PetscCall(TSSetSolution(ts, u)); // Allow TS to step over the final time // Note: This does not affect intermediate outputs, that // are always interpolated (in PetscMonitor) @@ -376,7 +357,7 @@ int PetscSolver::init() { ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP); } CHKERRQ(ierr); - ierr = TSMonitorSet(ts, PetscMonitor, this, nullptr); + PetscCall(TSMonitorSet(ts, PetscMonitor, this, nullptr)); if (ts_type != TSSUNDIALS) { // Get and configure the SNES nonlinear solver From 1aca504e46980c59d3ff77eeac4c2f13253eb90d Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 14 Aug 2025 23:01:53 -0700 Subject: [PATCH 06/18] petsc solver improvements - 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. --- src/solver/impls/petsc/petsc.cxx | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 126aa46532..6154cf39cd 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -346,6 +346,9 @@ int PetscSolver::init() { // Allow TS to recover from SNES failures PetscCall(TSSetMaxSNESFailures(ts, PETSC_UNLIMITED)); + // Recover from step rejections + PetscCall(TSSetMaxStepRejections(ts, PETSC_UNLIMITED)); + // Set the current solution PetscCall(TSSetSolution(ts, u)); // Allow TS to step over the final time @@ -395,9 +398,9 @@ int PetscSolver::init() { // Note: If the 'Amat' Jacobian is matrix free, SNESComputeJacobian // always updates its reference 'u' vector every nonlinear iteration SNESSetLagJacobian(snes, lag_jacobian); - SNESSetLagJacobianPersists(snes, PETSC_FALSE); + SNESSetLagJacobianPersists(snes, PETSC_TRUE); - SNESSetLagPreconditionerPersists(snes, PETSC_FALSE); + SNESSetLagPreconditionerPersists(snes, PETSC_TRUE); SNESSetLagPreconditioner(snes, 1); // Rebuild when Jacobian is rebuilt // Get and configure the KSP linear solver From e9025eed5a8e9f91d52795f232bb3483e9894320 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Thu, 14 Aug 2025 23:23:48 -0700 Subject: [PATCH 07/18] petsc: Add petsc headers and define PETSC_UNLIMITED Introduced in PETSc 3.22 (probably), define PETSC_UNLIMITED if it's undefined. Include more headers for used PETSc functionality. --- src/solver/impls/petsc/petsc.cxx | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 6154cf39cd..7c09b7d2a4 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -43,8 +43,15 @@ #include +#include "petscsys.h" +#include "petscts.h" #include "petscsnes.h" +#ifndef PETSC_UNLIMITED +// Introduced in PETSc 3.22 +#define PETSC_UNLIMITED -3 +#endif + class ColoringStencil { private: bool static isInSquare(int const i, int const j, int const n_square) { From 672172c7c7e72cb524b368d76519830e0954baf5 Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Fri, 15 Aug 2025 06:29:07 +0000 Subject: [PATCH 08/18] Apply clang-format changes --- src/solver/impls/petsc/petsc.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 7c09b7d2a4..051d6d3ee1 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -43,9 +43,9 @@ #include +#include "petscsnes.h" #include "petscsys.h" #include "petscts.h" -#include "petscsnes.h" #ifndef PETSC_UNLIMITED // Introduced in PETSc 3.22 From a0b6cca1ace1b45ce48bb60566ad7336ddc86cee Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 08:19:54 -0700 Subject: [PATCH 09/18] Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> --- src/solver/impls/petsc/petsc.cxx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 051d6d3ee1..55dd948b8a 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -49,7 +49,7 @@ #ifndef PETSC_UNLIMITED // Introduced in PETSc 3.22 -#define PETSC_UNLIMITED -3 +#define PETSC_UNLIMITED (-3) #endif class ColoringStencil { @@ -94,7 +94,7 @@ class ColoringStencil { // PETSc callback function for matrix-free preconditioner static PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { // Get the context - void* ctx; + void* ctx = nullptr; int ierr = PCShellGetContext(pc, &ctx); CHKERRQ(ierr); // Run the preconditioner @@ -178,7 +178,7 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo ierr = VecRestoreArrayRead(interpolatedX, &x); CHKERRQ(ierr); - if (s->call_monitors(output_time, i++, s->getNumberOutputSteps())) { + if (s->call_monitors(output_time, i++, s->getNumberOutputSteps()) != 0) { PetscFunctionReturn(1); } @@ -669,9 +669,9 @@ int PetscSolver::init() { for (int i = 0; i < nlocal; ++i) { // Assume all elements in the z direction are potentially coupled - d_nnz.emplace_back(d_nnz_map3d[i].size() * mesh->LocalNz + d_nnz.emplace_back((d_nnz_map3d[i].size() * mesh->LocalNz) + d_nnz_map2d[i].size()); - o_nnz.emplace_back(o_nnz_map3d[i].size() * mesh->LocalNz + o_nnz.emplace_back((o_nnz_map3d[i].size() * mesh->LocalNz) + o_nnz_map2d[i].size()); } } @@ -687,7 +687,7 @@ int PetscSolver::init() { // Mark non-zero entries output_progress.write("Marking non-zero Jacobian entries\n"); - PetscScalar val = 1.0; + PetscScalar const val = 1.0; for (int x = mesh->xstart; x <= mesh->xend; x++) { for (int y = mesh->ystart; y <= mesh->yend; y++) { @@ -706,14 +706,14 @@ int PetscSolver::init() { continue; } - int ind2 = ROUND(index(xi, yi, 0)); + int const ind2 = ROUND(index(xi, yi, 0)); if (ind2 < 0) { continue; // A boundary point } // Depends on all variables on this cell for (int j = 0; j < n2d; j++) { - PetscInt col = ind2 + j; + PetscInt const col = ind2 + j; ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); CHKERRQ(ierr); } @@ -721,7 +721,7 @@ int PetscSolver::init() { } // 3D fields for (int z = 0; z < mesh->LocalNz; z++) { - int ind = ROUND(index(x, y, z)); + int const ind = ROUND(index(x, y, z)); for (int i = 0; i < n3d; i++) { PetscInt row = ind + i; @@ -731,7 +731,7 @@ int PetscSolver::init() { // Depends on 2D fields for (int j = 0; j < n2d; j++) { - PetscInt col = ind0 + j; + PetscInt const col = ind0 + j; ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); CHKERRQ(ierr); } From dbff5bfc9009d939cf1a34bae5da1be64d3bae4e Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 12:24:12 -0700 Subject: [PATCH 10/18] petsc solver: Runs with SUNDIALS, some tidying 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. --- src/solver/impls/petsc/petsc.cxx | 122 +++++++++++++++---------------- src/solver/impls/petsc/petsc.hxx | 1 + 2 files changed, 59 insertions(+), 64 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 55dd948b8a..df559af5bd 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -25,27 +25,30 @@ #include "bout/build_defines.hxx" -#include "petsc.hxx" - #if BOUT_HAS_PETSC -#include +#include "petsc.hxx" + +#include +#include +#include +#include +#include #include #include #include #include +#include #include #include -#include -#include - -#include - -#include "petscsnes.h" -#include "petscsys.h" -#include "petscts.h" +#include +#include +#include +#include +#include +#include #ifndef PETSC_UNLIMITED // Introduced in PETSc 3.22 @@ -95,8 +98,7 @@ class ColoringStencil { static PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { // Get the context void* ctx = nullptr; - int ierr = PCShellGetContext(pc, &ctx); - CHKERRQ(ierr); + PetscCall(PCShellGetContext(pc, &ctx)); // Run the preconditioner PetscFunctionReturn(static_cast(ctx)->pre(x, y)); } @@ -141,28 +143,24 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo PetscFunctionReturn(0); } - PetscErrorCode ierr; PetscReal tfinal; static int i = 0; #if PETSC_VERSION_GE(3, 8, 0) - ierr = TSGetMaxTime(ts, &tfinal); - CHKERRQ(ierr); + PetscCall(TSGetMaxTime(ts, &tfinal)); #else - ierr = TSGetDuration(ts, nullptr, &tfinal); - CHKERRQ(ierr); + PetscCall(TSGetDuration(ts, nullptr, &tfinal)); #endif // Duplicate the solution vector X into a work vector Vec interpolatedX; - ierr = VecDuplicate(X, &interpolatedX); - CHKERRQ(ierr); + PetscCall(VecDuplicate(X, &interpolatedX)); // The internal timestepper may have stepped over multiple output times while (s->next_output <= t && s->next_output <= tfinal) { BoutReal output_time = t; if (s->interpolate) { - ierr = TSInterpolate(ts, s->next_output, interpolatedX); + int ierr = TSInterpolate(ts, s->next_output, interpolatedX); if (ierr != PETSC_SUCCESS) { throw BoutException("This PETSc TS does not support interpolation. Use a " "different method or set solver:interpolate=false"); @@ -172,11 +170,9 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo // Place the interpolated values into the global variables const PetscScalar* x; - ierr = VecGetArrayRead(interpolatedX, &x); - CHKERRQ(ierr); + PetscCall(VecGetArrayRead(interpolatedX, &x)); s->load_vars(const_cast(x)); - ierr = VecRestoreArrayRead(interpolatedX, &x); - CHKERRQ(ierr); + PetscCall(VecRestoreArrayRead(interpolatedX, &x)); if (s->call_monitors(output_time, i++, s->getNumberOutputSteps()) != 0) { PetscFunctionReturn(1); @@ -186,8 +182,7 @@ PetscErrorCode PetscMonitor(TS ts, PetscInt UNUSED(step), PetscReal t, Vec X, vo } // Done with vector, so destroy it - ierr = VecDestroy(&interpolatedX); - CHKERRQ(ierr); + PetscCall(VecDestroy(&interpolatedX)); PetscFunctionReturn(0); } @@ -214,8 +209,7 @@ PetscSolver::PetscSolver(Options* opts) .doc("Maximum number of function evaluations per SNES solve") .withDefault(10000)), maxl((*options)["maxl"].doc("Maximum number of linear iterations").withDefault(20)), - ts_type( - (*options)["ts_type"].doc("PETSc time integrator type").withDefault("beuler")), + ts_type((*options)["ts_type"].doc("PETSc time integrator type").withDefault("bdf")), adapt_type((*options)["adapt_type"] .doc("PETSc TSAdaptType timestep adaptation method") .withDefault("basic")), @@ -270,7 +264,6 @@ PetscSolver::~PetscSolver() { **************************************************************************/ int PetscSolver::init() { - PetscErrorCode ierr; TRACE("Initialising PETSc-dev solver"); @@ -320,7 +313,7 @@ int PetscSolver::init() { // Note: Vector atol and rtol not given PetscCall(TSSetTolerances(ts, atol, nullptr, rtol, nullptr)); if (ts_type == TSSUNDIALS) { -#if PETSC_HAS_SUNDIALS +#if PETSC_HAVE_SUNDIALS2 // The PETSc interface to SUNDIALS' CVODE TSSundialsSetType(ts, SUNDIALS_BDF); TSSundialsSetTolerance(ts, atol, rtol); @@ -362,11 +355,10 @@ int PetscSolver::init() { // Note: This does not affect intermediate outputs, that // are always interpolated (in PetscMonitor) if (interpolate) { - ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE); + PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_INTERPOLATE)); } else { - ierr = TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP); + PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); } - CHKERRQ(ierr); PetscCall(TSMonitorSet(ts, PetscMonitor, this, nullptr)); if (ts_type != TSSUNDIALS) { @@ -374,11 +366,8 @@ int PetscSolver::init() { // Note: SUNDIALS does not use PETSc SNES or KSP // https://petsc.org/release/manual/ts/#using-sundials-from-petsc - ierr = TSGetSNES(ts, &snes); - CHKERRQ(ierr); - - ierr = SNESSetType(snes, snes_type.c_str()); - CHKERRQ(ierr); + PetscCall(TSGetSNES(ts, &snes)); + PetscCall(SNESSetType(snes, snes_type.c_str())); // Line search if (line_search_type != "default") { @@ -455,7 +444,7 @@ int PetscSolver::init() { // Configure preconditioner PC pc; if (ts_type == TSSUNDIALS) { -#if PETSC_HAS_SUNDIALS +#if PETSC_HAVE_SUNDIALS2 TSSundialsGetPC(ts, &pc); #endif } else { @@ -480,8 +469,7 @@ int PetscSolver::init() { PCShellSetContext(pc, this); // Set name of preconditioner - ierr = PCShellSetName(pc, "PhysicsPreconditioner"); - CHKERRQ(ierr); + PetscCall(PCShellSetName(pc, "PhysicsPreconditioner")); } else { // Can't use preconditioner because no Jacobian matrix available PCSetType(pc, PCNONE); @@ -714,8 +702,7 @@ int PetscSolver::init() { // Depends on all variables on this cell for (int j = 0; j < n2d; j++) { PetscInt const col = ind2 + j; - ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); - CHKERRQ(ierr); + PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES)); } } } @@ -732,8 +719,7 @@ int PetscSolver::init() { // Depends on 2D fields for (int j = 0; j < n2d; j++) { PetscInt const col = ind0 + j; - ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); - CHKERRQ(ierr); + PetscCall(MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES)); } // Star pattern @@ -758,7 +744,7 @@ int PetscSolver::init() { // 3D fields on this cell for (int j = 0; j < n3d; j++) { PetscInt col = ind2 + j; - ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); + int ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); if (ierr != 0) { output.write("ERROR: {} {} : ({}, {}) -> ({}, {}) : {} -> {}\n", @@ -785,8 +771,7 @@ int PetscSolver::init() { // Test if the matrix is symmetric // Values are 0 or 1 so tolerance (1e-5) shouldn't matter PetscBool symmetric; - ierr = MatIsSymmetric(Jfd, 1e-5, &symmetric); - CHKERRQ(ierr); + PetscCall(MatIsSymmetric(Jfd, 1e-5, &symmetric)); if (!symmetric) { output_warn.write("Jacobian pattern is not symmetric\n"); } @@ -842,13 +827,15 @@ int PetscSolver::init() { TSAdaptType adapttype; TSAdaptGetType(adapt, &adapttype); output_info.write("TS Adapt Type : {}\n", adapttype); - SNESType snestype; - SNESGetType(snes, &snestype); - output_info.write("SNES Type : {}\n", snestype); - KSPType ksptype; - KSPGetType(ksp, &ksptype); - if (ksptype) { - output_info.write("KSP Type : {}\n", ksptype); + if (ts_type != TSSUNDIALS) { + SNESType snestype; + SNESGetType(snes, &snestype); + output_info.write("SNES Type : {}\n", snestype); + KSPType ksptype; + KSPGetType(ksp, &ksptype); + if (ksptype) { + output_info.write("KSP Type : {}\n", ksptype); + } } PCType pctype; PCGetType(pc, &pctype); @@ -927,8 +914,7 @@ PetscErrorCode PetscSolver::pre(Vec x, Vec y) { // Petsc's definition of Jacobian differs by a factor from SUNDIALS' // PETSc solves (scale + J)^-1 // SUNDIALS solves (I + gamma J)^-1 - PetscErrorCode ierr = VecScale(y, shift); - CHKERRQ(ierr); + PetscCall(VecScale(y, shift)); return 0; } @@ -951,18 +937,26 @@ void PetscSolver::updateColoring() { // Use the SNES function that is defined by the TS method // SNESTSFormFunction is defined in PETSc ts.c // The ctx pointer should be the TS object - MatFDColoringSetFunction(fdcoloring, bout::cast_MatFDColoringFn(SNESTSFormFunction), - ts); + MatFDColoringSetFunction(fdcoloring, + reinterpret_cast(SNESTSFormFunction), ts); MatFDColoringSetFromOptions(fdcoloring); MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); ISColoringDestroy(&iscoloring); // Replace the CTX pointer in SNES Jacobian - if (matrix_free_operator) { - // Use matrix-free calculation for operator, finite difference for preconditioner - SNESSetJacobian(snes, Jmf, Jfd, SNESComputeJacobianDefaultColor, fdcoloring); + if (ts_type == TSSUNDIALS) { +#if PETSC_HAVE_SUNDIALS2 + // The SUNDIALS interface calls TSGetIJacobian + // https://www.mcs.anl.gov/petsc/petsc-3.14/src/ts/impls/implicit/sundials/sundials.c + TSSetIJacobian(ts, Jfd, Jfd, solver_ijacobian, this); +#endif } else { - SNESSetJacobian(snes, Jfd, Jfd, SNESComputeJacobianDefaultColor, fdcoloring); + if (matrix_free_operator) { + // Use matrix-free calculation for operator, finite difference for preconditioner + SNESSetJacobian(snes, Jmf, Jfd, SNESComputeJacobianDefaultColor, fdcoloring); + } else { + SNESSetJacobian(snes, Jfd, Jfd, SNESComputeJacobianDefaultColor, fdcoloring); + } } } diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index 5108e8cb30..525ca93b40 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -47,6 +47,7 @@ class PetscSolver; #include #include +#include #include From 720dd2ef2f70e61b23cc540e0ecaf4bf8512feef Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 14:44:29 -0700 Subject: [PATCH 11/18] petsc solver: Fix Jacobian coloring with sundials method 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. --- src/solver/impls/petsc/petsc.cxx | 56 +++++++++++++++++++++++++++----- src/solver/impls/petsc/petsc.hxx | 7 ++++ 2 files changed, 54 insertions(+), 9 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index df559af5bd..551552a9b4 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -94,8 +94,9 @@ class ColoringStencil { } }; +namespace { // PETSc callback function for matrix-free preconditioner -static PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { +PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { // Get the context void* ctx = nullptr; PetscCall(PCShellGetContext(pc, &ctx)); @@ -105,14 +106,25 @@ static PetscErrorCode snesPCapply(PC pc, Vec x, Vec y) { // PETSc callback function, that evaluates the nonlinear // function being integrated by TS. -static PetscErrorCode solver_rhs(TS UNUSED(ts), BoutReal t, Vec globalin, Vec globalout, - void* f_data) { +PetscErrorCode solver_rhs(TS UNUSED(ts), BoutReal t, Vec globalin, Vec globalout, + void* f_data) { PetscFunctionBegin; auto* s = static_cast(f_data); s->rhs(t, globalin, globalout, false); PetscFunctionReturn(0); } +// Form function for use with SUNDIALS. +// This is needed because SNESTSFormFunction is not available +// PETSc error: No method snesfunction for TS of type sundials +PetscErrorCode solver_form_function(void* UNUSED(dummy), Vec U, Vec F, void* f_data) { + PetscFunctionBegin; + auto* s = static_cast(f_data); + s->formFunction(U, F); + PetscFunctionReturn(0); +} +} // namespace + // 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 @@ -131,6 +143,14 @@ PetscErrorCode solver_ijacobian(TS, BoutReal, Vec, Vec, PetscReal shift, Mat J, PetscFunctionReturn(0); } +PetscErrorCode solver_ijacobian_color(TS ts, PetscReal t, Vec U, Vec Udot, + PetscReal shift, Mat J, Mat B, void* ctx) { + auto* solver = static_cast(ctx); + solver->shift = shift; + + return TSComputeIJacobianDefaultColor(ts, t, U, Udot, shift, J, B, NULL); +} + // This function is called by the TS object every internal timestep // It is responsible for triggering interpolation and output at // the output time interval. @@ -887,6 +907,15 @@ PetscErrorCode PetscSolver::rhs(BoutReal t, Vec udata, Vec dudata, bool linear) PetscFunctionReturn(0); } +PetscErrorCode PetscSolver::formFunction(Vec U, Vec F) { + PetscFunctionBegin; + PetscCall(rhs(simtime, U, F, true)); // Can be linearised for coloring + + // shift * U - dU/dt + PetscCall(VecAXPBY(F, shift, -1.0, U)); + PetscFunctionReturn(0); +} + // Matrix-free preconditioner function PetscErrorCode PetscSolver::pre(Vec x, Vec y) { TRACE("PetscSolver::pre()"); @@ -934,11 +963,17 @@ void PetscSolver::updateColoring() { // Replace the old coloring with the new one MatFDColoringDestroy(&fdcoloring); MatFDColoringCreate(Jfd, iscoloring, &fdcoloring); - // Use the SNES function that is defined by the TS method - // SNESTSFormFunction is defined in PETSc ts.c - // The ctx pointer should be the TS object - MatFDColoringSetFunction(fdcoloring, - reinterpret_cast(SNESTSFormFunction), ts); + if (ts_type != TSSUNDIALS) { + // Use the SNES function that is defined by the TS method + // SNESTSFormFunction is defined in PETSc ts.c + // The ctx pointer should be the TS object + MatFDColoringSetFunction(fdcoloring, + reinterpret_cast(SNESTSFormFunction), ts); + } else { + // SNESTSFormFunction is not available for SUNDIALS + MatFDColoringSetFunction( + fdcoloring, reinterpret_cast(solver_form_function), this); + } MatFDColoringSetFromOptions(fdcoloring); MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); ISColoringDestroy(&iscoloring); @@ -948,7 +983,10 @@ void PetscSolver::updateColoring() { #if PETSC_HAVE_SUNDIALS2 // The SUNDIALS interface calls TSGetIJacobian // https://www.mcs.anl.gov/petsc/petsc-3.14/src/ts/impls/implicit/sundials/sundials.c - TSSetIJacobian(ts, Jfd, Jfd, solver_ijacobian, this); + // This sets the "TSMatFDColoring" property on the Jacobian, that is used in TSComputeIJacobianDefaultColor + PetscObjectCompose((PetscObject)Jfd, "TSMatFDColoring", (PetscObject)fdcoloring); + // Call a wrapper function that stores the shift in the PetscSolver + TSSetIJacobian(ts, Jfd, Jfd, solver_ijacobian_color, this); #endif } else { if (matrix_free_operator) { diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index 525ca93b40..a9b41333f7 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -67,6 +67,8 @@ public: /// Wrapper for the RHS function PetscErrorCode rhs(PetscReal t, Vec globalin, Vec globalout, bool linear); + /// Residual calculation. + PetscErrorCode formFunction(Vec U, Vec F); /// Wrapper for the preconditioner PetscErrorCode pre(Vec x, Vec y); @@ -76,6 +78,11 @@ public: friend PetscErrorCode solver_ijacobian(TS, BoutReal, Vec, Vec, PetscReal shift, Mat J, Mat Jpre, void* ctx); + // Wrapper around TSComputeIJacobianDefaultColor that saves the shift + // This is used to compute the Jacobian using coloring with SUNDIALS TS method. + friend PetscErrorCode solver_ijacobian_color(TS ts, PetscReal t, Vec U, Vec Udot, + PetscReal shift, Mat J, Mat B, void* ctx); + private: BoutReal shift; ///< Shift (alpha) parameter from TS Vec state; From 1e18bf23b15a93e9ec9a67b921582715bb72a469 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 15:14:05 -0700 Subject: [PATCH 12/18] petsc: Use C-style cast for MatFDColoringFn 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`). --- src/solver/impls/petsc/petsc.cxx | 18 ++++++++++++------ src/solver/impls/petsc/petsc.hxx | 3 +++ src/sys/petsclib.cxx | 1 + 3 files changed, 16 insertions(+), 6 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index 551552a9b4..ce33bc232c 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -44,6 +44,7 @@ #include #include +#include #include #include #include @@ -294,7 +295,8 @@ int PetscSolver::init() { // Get total problem size int neq; if (bout::globals::mpi->MPI_Allreduce(&nlocal, &neq, 1, MPI_INT, MPI_SUM, - BoutComm::get())) { + BoutComm::get()) + != 0) { throw BoutException("MPI_Allreduce failed!"); } @@ -967,12 +969,16 @@ void PetscSolver::updateColoring() { // Use the SNES function that is defined by the TS method // SNESTSFormFunction is defined in PETSc ts.c // The ctx pointer should be the TS object - MatFDColoringSetFunction(fdcoloring, - reinterpret_cast(SNESTSFormFunction), ts); + // + // Note: The cast is horrible but the function signature + // varies between PETSc versions in ways that break + // reinterpret_cast. + MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn)SNESTSFormFunction, ts); } else { - // SNESTSFormFunction is not available for SUNDIALS - MatFDColoringSetFunction( - fdcoloring, reinterpret_cast(solver_form_function), this); + // SNESTSFormFunction is not available for SUNDIALS. + // This solver_form_function needs to know the shift + // (SUNDIALS' gamma) that we capture in solver_ijacobian_color. + MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn)solver_form_function, this); } MatFDColoringSetFromOptions(fdcoloring); MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index a9b41333f7..b7857c7afd 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -47,6 +47,9 @@ class PetscSolver; #include #include +#include +#include +#include #include #include diff --git a/src/sys/petsclib.cxx b/src/sys/petsclib.cxx index 1830db89d9..bae53a8c7a 100644 --- a/src/sys/petsclib.cxx +++ b/src/sys/petsclib.cxx @@ -15,6 +15,7 @@ #include #include #include +#include #include #include From ca355a8a0393fc0b7f56aeb47c9a34ca17642ae6 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 15:44:40 -0700 Subject: [PATCH 13/18] petsc: Cast coloring function to void 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/petsclib.hxx | 1 + src/solver/impls/petsc/petsc.cxx | 25 +++++++++++++++++-------- src/solver/impls/petsc/petsc.hxx | 4 ++-- 3 files changed, 20 insertions(+), 10 deletions(-) diff --git a/include/bout/petsclib.hxx b/include/bout/petsclib.hxx index 760b18ff91..41d7618fd2 100644 --- a/include/bout/petsclib.hxx +++ b/include/bout/petsclib.hxx @@ -60,6 +60,7 @@ class Options; #define PETSC_HAVE_BROKEN_RECURSIVE_MACRO #include // IWYU pragma: export +#include #include #include "bout/boutexception.hxx" diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index ce33bc232c..abd4fbf64b 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -35,19 +35,25 @@ #include #include +#include #include #include #include #include #include #include +#include +#include #include #include +#include +#include #include -#include +#include #include #include +#include #include #include @@ -149,7 +155,7 @@ PetscErrorCode solver_ijacobian_color(TS ts, PetscReal t, Vec U, Vec Udot, auto* solver = static_cast(ctx); solver->shift = shift; - return TSComputeIJacobianDefaultColor(ts, t, U, Udot, shift, J, B, NULL); + return TSComputeIJacobianDefaultColor(ts, t, U, Udot, shift, J, B, nullptr); } // This function is called by the TS object every internal timestep @@ -765,7 +771,7 @@ int PetscSolver::init() { // 3D fields on this cell for (int j = 0; j < n3d; j++) { - PetscInt col = ind2 + j; + PetscInt const col = ind2 + j; int ierr = MatSetValues(Jfd, 1, &row, 1, &col, &val, INSERT_VALUES); if (ierr != 0) { @@ -952,13 +958,13 @@ PetscErrorCode PetscSolver::pre(Vec x, Vec y) { void PetscSolver::updateColoring() { // Re-calculate the coloring - MatColoring coloring = NULL; + MatColoring coloring{nullptr}; MatColoringCreate(Jfd, &coloring); MatColoringSetType(coloring, MATCOLORINGGREEDY); MatColoringSetFromOptions(coloring); // Calculate new index sets - ISColoring iscoloring = NULL; + ISColoring iscoloring{nullptr}; MatColoringApply(coloring, &iscoloring); MatColoringDestroy(&coloring); @@ -972,13 +978,16 @@ void PetscSolver::updateColoring() { // // Note: The cast is horrible but the function signature // varies between PETSc versions in ways that break - // reinterpret_cast. - MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn)SNESTSFormFunction, ts); + // reinterpret_cast and a C cast to MatFDColoringFn. + // This is the cast that PETSc examples use. + MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void))SNESTSFormFunction, + ts); } else { // SNESTSFormFunction is not available for SUNDIALS. // This solver_form_function needs to know the shift // (SUNDIALS' gamma) that we capture in solver_ijacobian_color. - MatFDColoringSetFunction(fdcoloring, (MatFDColoringFn)solver_form_function, this); + MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void))solver_form_function, + this); } MatFDColoringSetFromOptions(fdcoloring); MatFDColoringSetUp(Jfd, iscoloring, fdcoloring); diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index b7857c7afd..78f7d1a1e4 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -48,7 +48,7 @@ class PetscSolver; #include #include -#include +#include #include #include @@ -69,7 +69,7 @@ public: // These functions used internally (but need to be public) /// Wrapper for the RHS function - PetscErrorCode rhs(PetscReal t, Vec globalin, Vec globalout, bool linear); + PetscErrorCode rhs(BoutReal t, Vec udata, Vec dudata, bool linear); /// Residual calculation. PetscErrorCode formFunction(Vec U, Vec F); /// Wrapper for the preconditioner From f9179676563ea24ad7b963ca6a7e589a5dc3164f Mon Sep 17 00:00:00 2001 From: bendudson <219233+bendudson@users.noreply.github.com> Date: Fri, 15 Aug 2025 22:50:57 +0000 Subject: [PATCH 14/18] Apply clang-format changes --- src/solver/impls/petsc/petsc.cxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index abd4fbf64b..c0e3094481 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -980,13 +980,12 @@ void PetscSolver::updateColoring() { // varies between PETSc versions in ways that break // reinterpret_cast and a C cast to MatFDColoringFn. // This is the cast that PETSc examples use. - MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void))SNESTSFormFunction, - ts); + MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, ts); } else { // SNESTSFormFunction is not available for SUNDIALS. // This solver_form_function needs to know the shift // (SUNDIALS' gamma) that we capture in solver_ijacobian_color. - MatFDColoringSetFunction(fdcoloring, (PetscErrorCode (*)(void))solver_form_function, + MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))solver_form_function, this); } MatFDColoringSetFromOptions(fdcoloring); From 3da7deacd60d3595a231255a775fb4b0e9c97b8b Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 16:19:03 -0700 Subject: [PATCH 15/18] Try bout::cast_MatFDColoringFn Getting MatFDColoringSetFunction to work across PETSc versions seems to be unreasonably difficult. --- src/solver/impls/petsc/petsc.cxx | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index c0e3094481..ed48d87312 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -976,16 +976,14 @@ void PetscSolver::updateColoring() { // SNESTSFormFunction is defined in PETSc ts.c // The ctx pointer should be the TS object // - // Note: The cast is horrible but the function signature - // varies between PETSc versions in ways that break - // reinterpret_cast and a C cast to MatFDColoringFn. - // This is the cast that PETSc examples use. - MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, ts); + // Note: The function type casting here is horrible but necessary + MatFDColoringSetFunction(fdcoloring, bout::cast_MatFDColoringFn(SNESTSFormFunction), + ts); } else { // SNESTSFormFunction is not available for SUNDIALS. // This solver_form_function needs to know the shift // (SUNDIALS' gamma) that we capture in solver_ijacobian_color. - MatFDColoringSetFunction(fdcoloring, (PetscErrorCode(*)(void))solver_form_function, + MatFDColoringSetFunction(fdcoloring, bout::cast_MatFDColoringFn(solver_form_function), this); } MatFDColoringSetFromOptions(fdcoloring); From 64dd70e48adea2399c98cd1066bb9f4d7c92f386 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 17:58:11 -0700 Subject: [PATCH 16/18] petsc: More function casting I put on my robe and wizard hat... --- include/bout/petsc_interface.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/bout/petsc_interface.hxx b/include/bout/petsc_interface.hxx index 90dd188d8a..5c4f5ed0e8 100644 --- a/include/bout/petsc_interface.hxx +++ b/include/bout/petsc_interface.hxx @@ -573,7 +573,7 @@ PetscVector operator*(const PetscMatrix& mat, const PetscVector& vec) { namespace bout { template constexpr auto cast_MatFDColoringFn(T func) { - return func; + return reinterpret_cast(func); // NOLINT(*-reinterpret-cast) } } // namespace bout #else From 0413697ab87423a750a42bd109838092891c71b0 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 19:36:29 -0700 Subject: [PATCH 17/18] petsc: Yet more function casting --- include/bout/petsc_interface.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/include/bout/petsc_interface.hxx b/include/bout/petsc_interface.hxx index 5c4f5ed0e8..3ae02f41a6 100644 --- a/include/bout/petsc_interface.hxx +++ b/include/bout/petsc_interface.hxx @@ -6,9 +6,9 @@ * up a linear system. * ************************************************************************** - * Copyright 2019 C. MacMackin + * Copyright 2019 - 2025 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -573,7 +573,7 @@ PetscVector operator*(const PetscMatrix& mat, const PetscVector& vec) { namespace bout { template constexpr auto cast_MatFDColoringFn(T func) { - return reinterpret_cast(func); // NOLINT(*-reinterpret-cast) + return reinterpret_cast(func); // NOLINT(*-reinterpret-cast) } } // namespace bout #else From 56ad5e52dafaf1302a7b268b3d7cce7d122afc02 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Fri, 15 Aug 2025 20:11:30 -0700 Subject: [PATCH 18/18] petsc: Tighten tolerance in test-solver Tests pass with beuler/bdf/cn ts_types if tolerance is sufficiently tight. --- tests/integrated/test-solver/test_solver.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/integrated/test-solver/test_solver.cxx b/tests/integrated/test-solver/test_solver.cxx index 2e4345c8cc..a24db71750 100644 --- a/tests/integrated/test-solver/test_solver.cxx +++ b/tests/integrated/test-solver/test_solver.cxx @@ -100,8 +100,7 @@ int main(int argc, char** argv) { root["imexbdf2"]["adaptive"] = true; root["imexbdf2"]["adaptRtol"] = 1.e-5; - root["petsc"]["nout"] = 10000; - root["petsc"]["output_step"] = end / 10000; + root["petsc"]["rtol"] = 1e-7; root["snes"]["adaptive"] = true;