diff --git a/CMakeLists.txt b/CMakeLists.txt index 3caf4d8caf..573af32389 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -179,6 +179,7 @@ set(BOUT_SOURCES ./include/bout/paralleltransform.hxx ./include/bout/petsc_interface.hxx ./include/bout/petsc_operators.hxx + ./include/bout/petsc_jacobian.hxx ./include/bout/petsclib.hxx ./include/bout/physicsmodel.hxx ./include/bout/rajalib.hxx @@ -313,6 +314,7 @@ set(BOUT_SOURCES ./src/mesh/parallel_boundary_op.cxx ./src/mesh/parallel_boundary_region.cxx ./src/mesh/petsc_operators.cxx + ./src/mesh/petsc_jacobian.cxx ./src/mesh/surfaceiter.cxx ./src/mesh/tokamak_coordinates.cxx ./src/physics/gyro_average.cxx diff --git a/include/bout/petsc_jacobian.hxx b/include/bout/petsc_jacobian.hxx new file mode 100644 index 0000000000..120279c683 --- /dev/null +++ b/include/bout/petsc_jacobian.hxx @@ -0,0 +1,54 @@ +#pragma once + +#ifndef BOUT_PETSC_JACOBIAN_H +#define BOUT_PETSC_JACOBIAN_H + +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include + +/// Insert the nonzero pattern of @p sub into the variable block +/// (@p out_var, @p in_var) of the Jacobian @p Jfd. +/// +/// @p Jfd is a square matrix of size (n_evolving * nvars) where nvars is +/// inferred as Jfd_global_size / sub_global_size. Each nonzero (r, c) in +/// @p sub produces an entry at (r * nvars + out_var, c * nvars + in_var) +/// in @p Jfd. +/// +/// Limitation: this helper assumes that @p Jfd uses a uniform per-cell +/// interleaving of the form global_row = cell * nvars + var over the +/// evolving cells only. It is therefore not valid for solver layouts that mix +/// Field2D and Field3D variables, include evolving boundary variables, or use +/// any other non-uniform cell-to-row mapping. Callers must check that the +/// solver Jacobian ordering matches this assumption before using this helper. +/// +/// @p Jfd must already be preallocated. Entries are inserted with +/// INSERT_VALUES; MatAssemblyBegin/End must be called by the caller after +/// all insertions are complete. +/// +/// @param Jfd The Jacobian matrix to populate. Must be preallocated. +/// @param sub Evolving-cell submatrix providing the nonzero pattern. +/// @param out_var Row variable index in [0, nvars). +/// @param in_var Column variable index in [0, nvars). +void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var); + +/// @brief Insert the nonzero pattern of @p sub into every variable block of +/// @p Jfd. +/// +/// Equivalent to calling addOperatorSparsity(Jfd, sub, out_var, in_var) for +/// every combination of @p out_var and @p in_var in [0, nvars), where +/// @c nvars is inferred as @c Jfd_global / @c sub_global. +/// +/// The same layout restriction applies as for the block-specific overload: +/// @p Jfd must use the uniform cell * nvars + var ordering over +/// evolving cells only. +/// +/// @param Jfd The Jacobian matrix to populate. Must be preallocated. +/// @param sub Evolving-cell submatrix providing the nonzero pattern. +void addOperatorSparsity(Mat Jfd, Mat sub); + +#endif // BOUT_HAS_PETSC + +#endif // BOUT_PETSC_JACOBIAN_H diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx index 5a7df5b088..83069b9f60 100644 --- a/include/bout/petsc_operators.hxx +++ b/include/bout/petsc_operators.hxx @@ -290,6 +290,37 @@ public: } } + /// @brief Create a PETSc IS containing the global PETSc indices of all + /// locally owned evolving interior cells, in the order that + /// mapOwnedInteriorCells visits them. + /// + /// The returned IS selects the evolving subset of the full cell space C, + /// excluding inner/outer X-boundary cells and forward/backward parallel + /// boundary virtual cells. It is the correct IS to pass to + /// MatCreateSubMatrix when restricting a PetscCellOperator to the degrees + /// of freedom that the SNES solver actually integrates. + /// + /// The caller owns the returned IS and must call ISDestroy when finished. + /// + /// @returns A PETSc IS of local size equal to the number of locally owned + /// evolving cells, holding global PETSc row indices. + IS makeEvolvingIS() const; + + /// @brief Extract the evolving-cell submatrix from a cell-to-cell operator. + /// + /// Restricts @p op to the rows and columns that correspond to evolving + /// interior cells, discarding any rows or columns that belong to inner/outer + /// X-boundary cells or forward/backward parallel boundary virtual cells. + /// + /// The returned Mat is an independent copy (MAT_INITIAL_MATRIX): subsequent + /// modifications to @p op do not affect it. The caller owns the returned + /// Mat and must call MatDestroy when finished. + /// + /// @param op A cell-to-cell operator whose row and column space is the full + /// cell space C managed by this mapping. + /// @returns A square Mat of global size n_evolving × n_evolving. + Mat extractEvolvingSubmatrix(const PetscOperator& op) const; + friend PetscOperator makeNeumannOperator(const PetscCellMappingPtr& mapping, BoundaryDirection direction); @@ -678,6 +709,16 @@ public: /// low-level callers that need direct PETSc access. Mat raw() const { return *this->mat_operator; } + /// @brief Return the shared mapping for the operator's output (row) space. + const std::shared_ptr& getOutMapping() const { + return out_mapping; + } + + /// @brief Return the shared mapping for the operator's input (column) space. + const std::shared_ptr& getInMapping() const { + return in_mapping; + } + /// @brief Construct a diagonal operator from a vector of diagonal entries. /// /// Only available when @p OutSpaceTag == @p InSpaceTag (square, same-space @@ -1077,6 +1118,6 @@ PetscCellOperator makeNeumannOperator(const PetscCellMappingPtr& mapping, #warning PETSc not available. No PetscOperators. -#endif +#endif // BOUT_HAS_PETSC #endif // BOUT_PETSC_OPERATORS diff --git a/include/bout/solver.hxx b/include/bout/solver.hxx index 09ede6a32b..bfaf73b7d8 100644 --- a/include/bout/solver.hxx +++ b/include/bout/solver.hxx @@ -1,20 +1,10 @@ /************************************************************************** * Base class for all solvers. Specifies required interface functions * - * Changelog: - * - * 2009-08 Ben Dudson, Sean Farley - * * Major overhaul, and changed API. Trying to make consistent - * interface to PETSc and SUNDIALS solvers - * - * 2013-08 Ben Dudson - * * Added OO-style API, to allow multiple physics models to coexist - * For now both APIs are supported - * ************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2026 BOUT++ contributors * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -45,7 +35,6 @@ #include "bout/monitor.hxx" #include "bout/options.hxx" #include "bout/region.hxx" -#include "bout/unused.hxx" #include #include @@ -72,6 +61,10 @@ using TimestepMonitorFunc = int (*)(Solver* solver, BoutReal simtime, BoutReal l #include "bout/field2d.hxx" #include "bout/field3d.hxx" #include "bout/generic_factory.hxx" +#if BOUT_HAS_PETSC +#include "bout/petsc_interface.hxx" +#include "bout/petsc_operators.hxx" +#endif #include "bout/vector2d.hxx" #include "bout/vector3d.hxx" @@ -81,6 +74,8 @@ using TimestepMonitorFunc = int (*)(Solver* solver, BoutReal simtime, BoutReal l #include #include +#include +#include using SolverType = std::string; constexpr auto SOLVERCVODE = "cvode"; @@ -210,6 +205,30 @@ using RegisterUnavailableSolver = SolverFactory::RegisterUnavailableInFactory; */ class Solver { public: + /// Variable reference handle + class VarRef { + public: + static constexpr int AllValue = -1; + static constexpr int InvalidValue = -2; + + constexpr VarRef() = default; + + static constexpr VarRef All() { return VarRef(AllValue); } + static constexpr VarRef Invalid() { return VarRef(InvalidValue); } + + constexpr bool isAll() const { return value == AllValue; } + constexpr bool isInvalid() const { return value == InvalidValue; } + constexpr bool isConcrete() const { return value >= 0; } + constexpr bool isValid() const { return isAll() || isConcrete(); } + constexpr int index() const { return value; } + + private: + explicit constexpr VarRef(int value) : value(value) {} + + int value{InvalidValue}; + friend class Solver; + }; + Solver(Options* opts = nullptr); virtual ~Solver() = default; @@ -265,6 +284,21 @@ public: virtual void add(Vector3D& v, const std::string& name, const std::string& description = ""); + /// Get the solver-variable reference for a scalar field/component by name. + VarRef getVarRef(std::string_view name) const; + +#if BOUT_HAS_PETSC + /// Register a Jacobian sparsity contribution for all variable blocks. + bool addJacobianPattern(const PetscCellOperator& op); + + /// Register a Jacobian sparsity contribution for one or more variable blocks. + /// + /// Either variable reference may be VarRef::All(), which expands when the + /// Jacobian matrix is created during solver initialisation. + virtual bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var); +#endif + /// Returns true if constraints available virtual bool constraints() { return has_constraints; } @@ -394,7 +428,7 @@ protected: bool covariant{false}; /// For vectors bool evolve_bndry{false}; /// Are the boundary regions being evolved? std::string name; /// Name of the variable - std::string description{""}; /// Description of what the variable is + std::string description; /// Description of what the variable is }; /// A structure for iterating over fields @@ -688,6 +722,26 @@ private: /// Physics model being evolved PhysicsModel* model{nullptr}; +protected: +#if BOUT_HAS_PETSC + struct DeferredJacobianPattern { + bout::petsc::UniqueMat submatrix{new Mat{nullptr}}; + VarRef out_var; + VarRef in_var; + }; + + /// Queue a Jacobian-pattern contribution for PETSc-preconditioner-based solvers. + bool queueJacobianPattern(const PetscCellOperator& op, VarRef out_var, VarRef in_var); + + /// Insert any queued Jacobian-pattern contributions into an existing Jacobian matrix. + void applyQueuedJacobianPatterns(Mat Jfd) const; + + /// Check whether the current solver variable layout is compatible with + /// addOperatorSparsity(). + bool canApplyQueuedJacobianPatterns() const; +#endif + +private: /// Should non-split physics models be treated as diffusive? bool is_nonsplit_model_diffusive{true}; @@ -704,6 +758,10 @@ private: /// List of timestep monitor functions std::list timestep_monitors; +#if BOUT_HAS_PETSC + std::vector deferred_jacobian_patterns; +#endif + /// Should be run before user RHS is called void pre_rhs(BoutReal t); /// Should be run after user RHS is called diff --git a/manual/sphinx/index.rst b/manual/sphinx/index.rst index 1d4aaeab10..b63d5c1236 100644 --- a/manual/sphinx/index.rst +++ b/manual/sphinx/index.rst @@ -67,6 +67,7 @@ The documentation is divided into the following sections: :name: bout-interfaces user_docs/time_integration + user_docs/preconditioning user_docs/parallel-transforms user_docs/differential_operators user_docs/parallel_operators @@ -92,7 +93,6 @@ The documentation is divided into the following sections: :name: coordinates user_docs/coordinates - user_docs/preconditioning user_docs/BOUT_Gradperp_op .. toctree:: diff --git a/manual/sphinx/user_docs/parallel_operators.rst b/manual/sphinx/user_docs/parallel_operators.rst index b44e60ad70..23f4c73417 100644 --- a/manual/sphinx/user_docs/parallel_operators.rst +++ b/manual/sphinx/user_docs/parallel_operators.rst @@ -292,6 +292,86 @@ Current limitations include: - PETSc and 3D metrics are required - The grid metadata must be generated ahead of time +Using operator stencils in solver Jacobians +------------------------------------------- + +The ``PetscCellOperator`` objects can also be used to augment the Jacobian +pattern of PETSc-backed implicit solvers. This is useful when the default +solver stencil is too small for the operators used in the model, but the +coupling still has a sparse cell-to-cell structure. + +At a high level: + +- Build or reuse a ``PetscCellOperator`` +- Get solver variable references with ``solver->getVarRef(...)`` +- Register the operator with ``solver->addJacobianPattern(...)`` before + ``solver->init()`` + +For example: + +.. code-block:: C++ + + #if BOUT_HAS_PETSC + PetscOperators ops(mesh); + auto parallel = ops.getParallel(); + + auto n = solver->getVarRef("n"); + auto T = solver->getVarRef("T"); + + // Apply the operator stencil to every Jacobian block + solver->addJacobianPattern(parallel.Div_par_Grad_par); + + // Apply the same cell stencil only to dF_n / dT + solver->addJacobianPattern(parallel.Grad_par, n, T); + #endif + +The one-argument overload is equivalent to: + +.. code-block:: C++ + + solver->addJacobianPattern(op, Solver::VarRef::All(), Solver::VarRef::All()); + +The ``VarRef::All()`` sentinel expands when the solver Jacobian is created: + +- ``(All, All)`` inserts the stencil into every variable block +- ``(out, All)`` inserts it into one output-variable row block against all inputs +- ``(All, in)`` inserts it into all output-variable row blocks against one input +- ``(out, in)`` inserts it into a single block + +This interface is only available when BOUT++ is built with PETSc. Even then, +``addJacobianPattern(...)`` may return ``false`` if the chosen solver does not +use the PETSc preconditioner/Jacobian path. See :ref:`sec-time-integration` and +:ref:`sec-preconditioning` for the solver-side behavior. + + +How solver registration works +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +``PetscOperator`` tracks its row and column mappings, and exposes them through +the operator API. For Jacobian-pattern registration, BOUT++ uses the operator's +cell mapping to restrict the full cell-space matrix to the evolving subset of +cells that actually appear in the solver state vector. + +The resulting evolving-cell submatrix is then inserted into one or more +Jacobian variable blocks using the helper documented in +``bout/petsc_jacobian.hxx``. This means the operator contributes only sparsity +information here; the Jacobian entries are still computed later by the solver's +finite-difference machinery. + +Current limitations +~~~~~~~~~~~~~~~~~~~ + +This solver-Jacobian path currently assumes the Jacobian uses a uniform +per-cell interleaving of the evolving variables. In practice that means it is +currently intended for the same cases documented in +``bout/petsc_jacobian.hxx``: + +- Evolving ``Field3D`` variables only +- No evolving boundary cells in the solver state +- No solver layouts with mixed or non-uniform cell-to-row mappings + +If those assumptions do not hold, use the standard solver coloring stencil or a +different Jacobian/preconditioning strategy instead. See also -------- @@ -299,4 +379,6 @@ See also - :ref:`sec-diffops` - :ref:`sec-fci` - :ref:`sec-parallel-bc-fci` +- :ref:`sec-time-integration` +- :ref:`sec-preconditioning` - :ref:`sec-petsc` diff --git a/manual/sphinx/user_docs/preconditioning.rst b/manual/sphinx/user_docs/preconditioning.rst index 79d8e050ae..78a25814bd 100644 --- a/manual/sphinx/user_docs/preconditioning.rst +++ b/manual/sphinx/user_docs/preconditioning.rst @@ -1,351 +1,284 @@ .. default-role:: math -====================== -BOUT++ preconditioning -====================== +.. _sec-preconditioning: -:Author: B.Dudson, University of York +===================================== +BOUT++ preconditioning and Jacobians +===================================== -Introduction -============ +Implicit solvers repeatedly solve linearised systems of the form +`(I - \gamma J) x = b`. BOUT++ supports several ways to help those solves: -This manual describes some of the ways BOUT++ could (and in some cases -does) support preconditioning, Jacobian calculations and other methods -to speed up simulations. This manual assumes that you’re familiar with -how BOUT++ works internally. +- a user-supplied preconditioner registered with ``solver->setPrecon(...)`` +- a user-supplied Jacobian-vector product registered with + ``solver->setJacobian(...)`` +- PETSc finite-difference Jacobians, usually accelerated with coloring +- extra Jacobian sparsity supplied from ``PetscCellOperator`` stencils with + ``solver->addJacobianPattern(...)`` -Some notation: The ODE being solved is of the form +This page starts with the user-facing workflow and then moves into the current +implementation details and limitations. -.. math:: {\frac{\partial {\mathbf{f}}}{\partial t}} = {\mathbf{F}}\left({\mathbf{f}}\right) -Here the state vector `f = \left(f_0, f_1, f_2, \ldots\right)^T` -is a vector containing the evolving (3D) variables -`f_i\left(x,y,z\right)`. +Choosing an approach +-------------------- -The Jacobian of this system is then +The most useful question is usually not "which feature exists?" but "which +approximation is easiest to provide for this model?". -.. math:: {\mathbb{J}}= {\frac{\partial {\mathbf{F}}}{\partial {\mathbf{f}}}} +- Use matrix-free methods when assembling a Jacobian is too expensive or the + coupling is not well represented by a sparse local stencil. +- Add a user preconditioner when you know an approximate inverse for the stiff + part of the physics. +- Use PETSc finite-difference Jacobians when the coupling is sparse and local + enough that a structured sparsity pattern is a good approximation. +- Add ``PetscCellOperator`` sparsity when the solver's default stencil is too + small, but the coupling is still naturally represented by a sparse + cell-to-cell operator. -The order of the elements in the vector `{\mathbf{f}}` -is determined in the solver code and SUNDIALS, so here just assume that -there exists a map `\mathbb{I}` between a global index `k` -and (variable, position) i.e. `\left(i,x,y,z\right)` +The solver options that choose between these modes are described in +:ref:`sec-time-integration`. PETSc option prefixes are described in +:ref:`sec-petsc`. -.. math:: \mathbf{I} : \left(i,x,y,z\right) \mapsto k -and its inverse +User-supplied preconditioners +----------------------------- -.. math:: \mathbf{I}^{-1} : k \mapsto \left(i,x,y,z\right) +A user preconditioner is a callback that approximately solves +`(I - \gamma J) x = b` for the solver's current linearisation. The callback +signature is: -Some problem-specific operations which can be used to speed up the -timestepping +.. code-block:: C++ -#. Jacobian-vector multiply: Given a vector, multiply it by - `{\mathbb{J}}` + int precon(BoutReal t, BoutReal gamma, BoutReal delta); -#. Preconditioner multiply: Given a vector, multiply by an approximate - inverse of `\mathbb{M} = \mathbb{I} - \gamma\mathbb{J}` +Register it during ``PhysicsModel::init``: -#. Calculate the stencils i.e. non-zero elements in - `{\mathbb{J}}` +.. code-block:: C++ -#. Calculate the non-zero elements of `{\mathbb{J}}` + int init(bool restarting) override { + solver->setPrecon(precon); + return 0; + } -Physics problems -================ +Inside the callback: -Some interesting physics problems of increasing difficulty +- the current state is stored in the usual evolving variables +- the vector to be preconditioned is stored in the corresponding + ``ddt(variable)`` fields +- the callback should overwrite those ``ddt(...)`` fields with the + preconditioned result -Resistive drift-interchange instability ---------------------------------------- +For CVODE, enable a user preconditioner with for example: -A “simple” test problem of 2 fields, which results in non-trivial -turbulent results. Supports resistive drift wave and interchange -instabilities. +.. code-block:: cfg -.. math:: - - \begin{aligned} - {\frac{\partial N_i}{\partial t}} + {{\mathbf{v}}_E}\cdot\nabla N_i &=& 0 \\ - {\frac{\partial \omega}{\partial t}} + {{\mathbf{v}}_E}\cdot\nabla\omega &=& 2\omega_{ci}{\mathbf{b}}\times\kappa\cdot\nabla P + N_iZ_i e\frac{4\pi V_A^2}{c^2}\nabla_{||}j_{||} \\ - \nabla_\perp^2\omega / N_i &=& \phi \\ - 0.51\nu_{ei}j_{||} &=& \frac{e}{m_e}\partial_{||}\phi + \frac{T_e}{N_i m_e}\partial_{||} N_i\end{aligned} - -Reduced 3-field MHD -------------------- - -This is a 3-field system of pressure `P`, magnetic flux -`\psi` and vorticity `U`: + [solver] + type = cvode + cvode_precon_method = user -.. math:: - - {\mathbf{f}} = \left(\begin{array}{c} - P \\ - \psi \\ - U - \end{array}\right) +The ``examples/test-precon`` case shows the standard pattern. For the simple +wave system .. math:: - \begin{aligned} - {\frac{\partial \psi}{\partial t}} &=& -\frac{1}{B_0}\nabla_{||}\phi \\ - &=& -\frac{1}{B_0}\left[{\mathbf{b}}_0 - \left({\mathbf{b}}_0\times\nabla\psi\right)\right]\cdot\nabla\phi \\ - &=& -\frac{1}{B_0}{\mathbf{b}}_0\cdot\nabla\phi - \frac{1}{B_0}\left({\mathbf{b}}_0\times\nabla\phi\right)\cdot\nabla\psi \\ - \Rightarrow \frac{d \psi}{dt} &=& -\frac{1}{B_0}{\mathbf{b}}_0\cdot\nabla \phi\end{aligned} + \frac{\partial u}{\partial t} = \partial_{||} v + \qquad + \frac{\partial v}{\partial t} = \partial_{||} u -The coupled set of equations to be solved are therefore +one useful approximation factors `(I - \gamma J)^{-1}` into inexpensive pieces: .. math:: - \begin{aligned} - \frac{1}{B_0}\nabla_\perp^2\phi &=& U \\ - \left({\frac{\partial }{\partial t}} + {\mathbf{v}}_E\cdot\nabla\right)\psi &=& -\frac{1}{B_0}{\mathbf{b}}_0\cdot\nabla\phi \\ - \left({\frac{\partial }{\partial t}} + {\mathbf{v}}_E\cdot\nabla\right)P &=& 0 \\ - \left({\frac{\partial }{\partial t}} + {\mathbf{v}}_E\cdot\nabla\right)U &=& \frac{1}{\rho}B_0^2\left[{\mathbf{b}}_0 - \left({\mathbf{b}}_0\times\nabla\psi\right)\right]\cdot\left(\frac{J_{||0}}{B_0} - \frac{1}{\mu_0}\nabla_\perp^2\psi\right) \nonumber \\ - &+& \frac{1}{\rho}{\mathbf{b}}_0\times{\mathbf{\kappa}}_0\cdot\nabla P \\ - {\mathbf{v}}_E &=& \frac{1}{B_0}{\mathbf{b}}_0\times\nabla\phi\end{aligned} + \left(\begin{array}{cc} + 1 & -\gamma \partial_{||} \\ + -\gamma \partial_{||} & 1 + \end{array}\right)^{-1} + = + \left(\begin{array}{cc} + 1 & \gamma \partial_{||} \\ + 0 & 1 + \end{array}\right) + \left(\begin{array}{cc} + 1 & 0 \\ + 0 & (1 - \gamma^2 \partial_{||}^2)^{-1} + \end{array}\right) + \left(\begin{array}{cc} + 1 & 0 \\ + \gamma \partial_{||} & 1 + \end{array}\right) -The Jacobian of this system is therefore: +In BOUT++ that becomes a sequence of field operations and an ``InvertPar`` solve: -.. math:: - :label: eq:mhdjacobian +.. code-block:: C++ - \mathbb{J} = - \left[ \begin{array}{c|c|c} - \color{blue}{-{\mathbf{v}}_E\cdot\nabla} & 0 & \left[{\mathbf{b}}_0\times\nabla\left(P_0 + \color{blue}{P}\right)\cdot\nabla\right]\nabla_\perp^{-2} \\ - \hline - 0 & \color{blue}{-{\mathbf{v}}_E\cdot\nabla} & \left({\mathbf{b}}_0\cdot\nabla\right)\nabla_\perp^{-2} \\ - \hline - 2{\mathbf{b}}_0\times{\mathbf{\kappa}}_0\cdot\nabla& -\frac{B_0^2}{\mu_0\rho}\left({\mathbf{b}}_0 \color{blue}{-{\mathbf{b}}_0\times\nabla\psi}\right)\cdot\nabla\nabla_\perp^2& \color{blue}{-{\mathbf{v}}_E\cdot\nabla} \\ - & + \frac{B_0^2}{\rho}\left[{\mathbf{b}}_0\times\nabla\left(\frac{J_{||0}}{B_0}\right)\right]\cdot\nabla & \\ - & + \color{blue}{\frac{B_0^2}{\mu_0\rho}\nabla\left(\nabla_\perp^2\psi\right)\cdot\left({\mathbf{b}}_0\times\nabla\right)} & - \end{array}\right] + int precon(BoutReal t, BoutReal gamma, BoutReal delta) { + mesh->communicate(ddt(u)); + ddt(v) = gamma * Grad_par(ddt(u)) + ddt(v); -Where the blue terms are only included in nonlinear simulations. + inv->setCoefB(-SQ(gamma)); + ddt(v) = inv->solve(ddt(v)); -This Jacobian has large dense blocks because of the Laplacian inversion -terms (involving `\nabla_\perp^{-2}` which couples together all -points in an X-Z plane. The way to make `{\mathbb{J}}` -sparse is to solve `\phi` as a constraint (using e.g. the IDA -solver) which moves the Laplacian inversion to the preconditioner. + mesh->communicate(ddt(v)); + ddt(u) = ddt(u) + gamma * Grad_par(ddt(v)); -Solving `\phi` as a constraint ------------------------------------- + ddt(u).applyBoundary("dirichlet"); + ddt(v).applyBoundary("dirichlet"); + return 0; + } -The evolving state vector becomes +This is the main design pattern for analytic preconditioners: keep only the +stiff physics, approximate the rest, and return something that is cheap but +captures the dominant fast scales. -.. math:: - {\mathbf{f}} = \left(\begin{array}{c} - P \\ - \psi \\ - U \\ - \phi - \end{array}\right) +Jacobian-vector products +------------------------ -UEDGE equations ---------------- +If the model can apply the Jacobian to a vector directly, register that +operation with ``solver->setJacobian(...)``. This is primarily useful for +matrix-free methods, where the solver needs Jacobian-vector products but does +not assemble a full sparse Jacobian matrix. -The UEDGE benchmark is a 4-field model with the following equations: +For problems where a sparse matrix is more useful than a matrix-free operator, +PETSc-backed solvers can instead assemble the Jacobian by finite differences. -.. math:: - \begin{aligned} - {\frac{\partial N_i}{\partial t}} + {V_{||}}\partial_{||}N_i &=& -N_i\nabla_{||}{V_{||}}+\nabla_\psi\left(D_\perp \partial_\psi N_i\right) \\ - {\frac{\partial \left(N_i{V_{||}}\right)}{\partial t}} + {V_{||}}\partial_{||}\left(N_i{V_{||}}\right) &=& -\partial_{||}P + \nabla_\psi\left(N_i\mu_\perp\partial_\psi{V_{||}}\right) \\ - \frac{3}{2}{\frac{\partial }{\partial t}}\left(N_iT_e\right) &=& \nabla_{||}\left(\kappa_e\partial_{||}T_e\right) + \nabla_\psi\left(N_i\chi_\perp\partial_\perp T_e\right) \\ - \frac{3}{2}{\frac{\partial }{\partial t}}\left(N_iT_i\right) &=& \nabla_{||}\left(\kappa_i\partial_{||}T_i\right) + \nabla_\psi\left(N_i\chi_\perp\partial_\perp T_i\right)\end{aligned} +PETSc finite-difference Jacobians +--------------------------------- -This set of equations is good in that there is no inversion needed, and -so the Jacobian is sparse everywhere. The state vector is +When a PETSc-backed implicit solver runs with ``matrix_free = false``, BOUT++ +can assemble an explicit sparse Jacobian by finite differences. In most cases +``use_coloring = true`` is the right choice: -.. math:: +.. code-block:: cfg - {\mathbf{f}} = \left(\begin{array}{c} - N_i \\ - {V_{||}}\\ - T_e \\ - T_i \\ - \end{array}\right) + [solver] + matrix_free = false + use_coloring = true + lag_jacobian = 5 -The Jacobian is: + stencil:taxi = 2 + stencil:square = 0 + stencil:cross = 0 -.. math:: +The default pattern comes from a solver-side stencil and is then passed to PETSc +coloring, which groups independent columns so the Jacobian can be built with far +fewer RHS evaluations than a brute-force finite-difference calculation. - \mathbb{J} = - \left( \begin{array}{c|c|c|c} - -{V_{||}}\partial_{||} - \nabla_{||}{V_{||}}+ \nabla_\psi D_\perp\partial_\psi & -\partial_{||}N_i - N_i\nabla_{||} & 0 & 0 \\ - -\frac{1}{N_i}{\frac{\partial {V_{||}}}{\partial t}} - \frac{1}{N_i}{V_{||}}{\mathbb{J}}_{N_iN_i} & & & - \end{array}\right) +The most important current assumption is that the Jacobian is sparse and local +in the solver ordering. If the RHS includes long-range couplings such as Fourier +transforms, matrix inversions, or other nonlocal operators, the default stencil +can miss real dependencies and the implicit solve may fail to converge. In those +cases the usual alternatives are: -If instead the state vector is +- stay matrix-free +- reformulate the problem so the long-range solve becomes a constraint +- supply a better sparsity pattern when one exists -.. math:: - - {\mathbf{f}} = \left(\begin{array}{c} - N_i \\ - N_i{V_{||}}\\ - N_iT_e \\ - N_iT_i \\ - \end{array}\right) - -then the Jacobian is - -.. Result is missing! - -2-fluid turbulence ------------------- - -Jacobian-vector multiply -======================== - -This is currently implemented into the CVODE (SUNDIALS) solver. - -Preconditioner-vector multiply -============================== +The main tuning options are: -.. _reduced-3-field-mhd-1: +- ``use_coloring = false`` to fall back to a brute-force finite-difference + Jacobian +- ``stencil:taxi``, ``stencil:square``, and ``stencil:cross`` to widen the + default solver stencil +- ``force_symmetric_coloring = true`` to make the coloring pattern symmetric +- ``lag_jacobian`` to reuse the assembled Jacobian across successive nonlinear + iterations -Reduced 3-field MHD -------------------- -The matrix `\mathbb{M}` to be inverted can therefore be written +Augmenting Jacobian structure with PetscOperator +------------------------------------------------ -.. math:: +If the coupling you want to expose to the solver is already available as a +``PetscCellOperator``, you can register its sparsity pattern directly with the +solver before ``solver->init()``. - \mathbb{M} = - \left[ \begin{array}{ccc} - \mathbb{D} & 0 & \mathbb{U}_P \\ - 0 & \mathbb{D} & \mathbb{U}_\psi \\ - \mathbb{L}_P & \mathbb{L}_\psi & \mathbb{D} - \end{array}\right] +For example: -where +.. code-block:: C++ -.. math:: \mathbb{D} = \mathbb{I} \color{blue}{+ \gamma{\mathbf{v}}_E\cdot\nabla} + #if BOUT_HAS_PETSC + PetscOperators ops(mesh); + auto parallel = ops.getParallel(); -For small flow velocities, the inverse of `\mathbb{D}` can be -approximated using the Binomial theorem: + auto n = solver->getVarRef("n"); + auto T = solver->getVarRef("T"); -.. math:: - :label: eq:dapprox + // Insert the operator stencil into every Jacobian block + solver->addJacobianPattern(parallel.Div_par_Grad_par); - \mathbb{D}^{-1} \simeq \mathbb{I} \color{blue}{- \gamma{\mathbf{v}}_E\cdot\nabla} + // Insert the operator stencil only into dF_n / dT + solver->addJacobianPattern(parallel.Grad_par, n, T); + #endif -Following [chacon-2008]_, [chacon-2002]_, `\mathbb{M}` can be -re-written as +The one-argument form is shorthand for: -.. math:: +.. code-block:: C++ - \mathbb{M} = - \left[ \begin{array}{cc} - \mathbb{E} & \mathbb{U} \\ - \mathbb{L} & \mathbb{D} - \end{array}\right] \qquad \mathbb{E} = - \left[ \begin{array}{cc} - \mathbb{D} & 0 \\ - 0 & \mathbb{D} - \end{array}\right] \qquad \mathbb{U} = - \left(\begin{array}{c} - \mathbb{U}_P \\ - \mathbb{U}_\psi - \end{array}\right) \qquad \mathbb{L} = \left(\mathbb{L}_P \quad \mathbb{L}_\psi\right) - -The Schur factorization of `\mathbb{M}` yields ([chacon-2008]_) + solver->addJacobianPattern(op, Solver::VarRef::All(), Solver::VarRef::All()); -.. math:: +``Solver::VarRef::All()`` expands when the solver Jacobian is created: - \mathbb{M}^{-1} = - \left[ \begin{array}{cc} - \mathbb{E} & \mathbb{U} \\ - \mathbb{L} & \mathbb{D} - \end{array}\right]^{-1} = - \left[ \begin{array}{cc} - \mathbb{I} & -\mathbb{E}^{-1}\mathbb{U} \\ - 0 & \mathbb{I} - \end{array}\right] - \left[ \begin{array}{cc} - \mathbb{E}^{-1} & 0 \\ - 0 & \mathbb{P}_{Schur}^{-1} - \end{array}\right] - \left[ \begin{array}{cc} - \mathbb{I} & 0 \\ - -\mathbb{L}\mathbb{E}^{-1} & \mathbb{I} - \end{array}\right] - -Where -`\mathbb{P}_{Schur} = \mathbb{D} - \mathbb{L}\mathbb{E}^{-1}\mathbb{U}` -is the Schur complement. Note that this inversion is exact so far. Since -`\mathbb{E}` is block-diagonal, and `\mathbb{D}` can be -easily approximated using equation :eq:`eq:dapprox`, this -simplifies the problem to inverting `\mathbb{P}_{Schur}`, which is -much smaller than `\mathbb{M}`. - -A possible approximation to `\mathbb{P}_{Schur}` is to neglect: - -- All drive terms - - - the curvature term `\mathbb{L}_P` - - - the `J_{||0}` term in `\mathbb{L}_\psi` - -- All nonlinear terms (blue terms in equation :eq:`eq:mhdjacobian`), - including perpendicular terms (so `\mathbb{D} = \mathbb{I}`) - -This gives +- ``(All, All)`` inserts the stencil into every variable block +- ``(out, All)`` inserts it into one output-variable row block against all inputs +- ``(All, in)`` inserts it into all output-variable row blocks against one input +- ``(out, in)`` inserts it into a single block -.. math:: +This interface is guarded by ``BOUT_HAS_PETSC``. Even in PETSc-enabled builds, +``addJacobianPattern(...)`` may return ``false`` if the chosen solver does not +use the PETSc preconditioner/Jacobian path. - \begin{aligned} - \mathbb{P}_{Schur} &\simeq& \mathbb{I} + \gamma^2 \frac{B_0^2}{\mu_0\rho}\left({\mathbf{b}}_0\cdot\nabla\right)\nabla_\perp^2\left({\mathbf{b}}_0\cdot\nabla\right)\nabla_\perp^{-2} \nonumber \\ - &\simeq& \mathbb{I} + \gamma^2 V_A^2 \left({\mathbf{b}}_0\cdot\nabla\right)^2\end{aligned} -Where the commutation of parallel and perpendicular derivatives is also -an approximation. This remaining term is just the shear Alfvén wave -propagating along field-lines, the fastest wave supported by these -equations. +How operator-supplied sparsity is applied +----------------------------------------- -Stencils -======== +The operator itself lives on the full PETSc cell space, including boundary and +virtual-cell entries used by the operator construction. The solver Jacobian does +not: it contains only the evolving degrees of freedom. -Jacobian calculation -==================== +To bridge those two layouts, BOUT++ uses the operator's stored cell mapping to +extract the evolving-cell submatrix first. That restricted submatrix is then +inserted into one or more solver Jacobian blocks using the helper documented in +``bout/petsc_jacobian.hxx``. -The (sparse) Jacobian matrix elements can be calculated automatically -from the physics code by keeping track of the (linearised) operations -going through the RHS function. +The important consequence is that ``addJacobianPattern(...)`` contributes only +the nonzero structure. The Jacobian entries themselves are still computed later +by the solver's finite-difference machinery. -For each point, keep the value (as usual), plus the non-zero elements in -that row of `{\mathbb{J}}` and the constant: result = -Ax + b Keep track of elements using product rule. +At the moment the sequence is: -:: +1. Model or component code calls ``solver->addJacobianPattern(...)`` during + setup. +2. The solver stores the restricted operator submatrix and the requested output + and input variable references. +3. During ``solver->init()``, a supported PETSc-backed solver creates its + default Jacobian pattern. +4. The queued operator patterns are replayed into that Jacobian. +5. PETSc coloring and finite-difference Jacobian assembly use the augmented + sparsity pattern. - class Field3D { - data[ngx][ngy][ngz]; // The data as now - - int JacIndex; // Variable index in Jacobian - SparseMatrix *jac; // Set of rows for indices (JacIndex,*,*,*) - }; -JacIndex is set by the solver, so for the system +Current limitations +------------------- -.. math:: +The helper used to insert operator sparsity into the solver Jacobian currently +assumes a uniform per-cell variable ordering of the form +``global_row = cell * nvars + var`` over evolving cells only. That means this +path is currently intended for: - {\mathbf{f}} = \left(\begin{array}{c} - P \\ - \psi \\ - U - \end{array}\right) +- evolving ``Field3D`` variables +- no evolving boundary cells +- solver layouts that use a uniform per-cell interleaving of variables -``P.JacIndex = 0``, ``psi.JacIndex = 1`` and ``U.JacIndex = 2``. All -other fields are given ``JacIndex = -1``. +It is not currently valid for mixed ``Field2D`` and ``Field3D`` solver layouts, +for layouts that include evolving boundary variables, or for any other +non-uniform cell-to-row mapping. -SparseMatrix stores the non-zero Jacobian components for the set of rows -corresponding to this variable. Evolving variables do not have an -associated ``SparseMatrix`` object, but any fields which result from -operations on evolving fields will have one. -.. [chacon-2008] L. Chacón, An optimal, parallel, fully implicit Newton-Krylov solver for three-dimensional viscoresistive magnetohydrodynamics, POP, 2008, 15, 056103 +See also +-------- -.. [chacon-2002] L. Chacón, D.A. Knoll, and J.M. Finn, An Implicit, Nonlinear Reduced Resistive MHD Solver, JCP, 2002, 178, 15-36 +- :ref:`sec-time-integration` +- :ref:`sec-parallel-operators-petsc-fci` +- :ref:`sec-petsc` diff --git a/manual/sphinx/user_docs/time_integration.rst b/manual/sphinx/user_docs/time_integration.rst index 097b62eeda..58ba7c4694 100644 --- a/manual/sphinx/user_docs/time_integration.rst +++ b/manual/sphinx/user_docs/time_integration.rst @@ -721,6 +721,14 @@ The default and recommended approach for most problems: The coloring algorithm exploits the sparse structure of the Jacobian to reduce the number of function evaluations needed for finite differencing. +If the default solver stencil is too narrow for part of the physics, PETSc +cell-space operators can add extra sparsity information before coloring is set +up. Register those contributions during model setup with +``solver->addJacobianPattern(...)``. This is particularly useful for +``PetscCellOperator`` stencils assembled from the FCI support-operator tools in +:ref:`sec-parallel-operators-petsc-fci`. See :ref:`sec-preconditioning` for the +full solver-side behavior and current limitations. + Jacobian coloring stencil ^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -997,263 +1005,42 @@ Finally, delete the model and solver when finished:: **Note:** If an ODE needs to be solved multiple times, at the moment it is recommended to delete the solver, and create a new one each time. -.. _sec-preconditioning: - Preconditioning --------------- -At every time step, an implicit scheme such as BDF has to solve a -non-linear problem to find the next solution. This is usually done using -Newton’s method, each step of which involves solving a linear (matrix) -problem. For :math:`N` evolving variables is an :math:`N\times N` matrix -and so can be very large. By default matrix-free methods are used, in -which the Jacobian :math:`\mathcal{J}` is approximated by finite -differences (see next subsection), and so this matrix never needs to be -explicitly calculated. Finding a solution to this matrix can still be -difficult, particularly as :math:`\delta t` gets large compared with -some time-scales in the system (i.e. a stiff problem). - -A preconditioner is a function which quickly finds an approximate -solution to this matrix, speeding up convergence to a solution. A -preconditioner does not need to include all the terms in the problem -being solved, as the preconditioner only affects the convergence rate -and not the final solution. A good preconditioner can therefore -concentrate on solving the parts of the problem with the fastest -time-scales. - -A simple example [1]_ is a coupled wave equation, solved in the -``test-precon`` example code: - -.. math:: - - \frac{\partial u}{\partial t} = \partial_{||}v \qquad \frac{\partial - v}{\partial t} = \partial_{||} u - -First, calculate the Jacobian of this set of equations by taking -partial derivatives of the time-derivatives with respect to each of the -evolving variables - -.. math:: - - \mathcal{J} = (\begin{array}{cc} - \frac{\partial}{\partial u}\frac{\partial u}{\partial t} & - \frac{\partial}{\partial v}\frac{\partial u}{\partial t}\\ - \frac{\partial}{\partial u}\frac{\partial v}{\partial t} & - \frac{\partial}{\partial v}\frac{\partial v}{\partial t} - \end{array} - ) = (\begin{array}{cc} - 0 & \partial_{||} \\ - \partial_{||} & 0 - \end{array} - ) - -In this case :math:`\frac{\partial u}{\partial t}` doesn’t depend on -:math:`u` nor :math:`\frac{\partial v}{\partial t}` on :math:`v`, so the -diagonal is empty. Since the equations are linear, the Jacobian doesn’t -depend on :math:`u` or :math:`v` and so - -.. math:: - - \frac{\partial}{\partial t}(\begin{array}{c} u \\ - v \end{array}) = \mathcal{J} (\begin{array}{c} u \\ - v \end{array} ) - -In general for non-linear functions :math:`\mathcal{J}` gives the -change in time-derivatives in response to changes in the state variables -:math:`u` and :math:`v`. - -In implicit time stepping, the preconditioner needs to solve an equation - -.. math:: - - \mathcal{I} - \gamma \mathcal{J} - -where :math:`\mathcal{I}` is the identity matrix, and :math:`\gamma` -depends on the time step and method (e.g. :math:`\gamma = \delta t` for -backwards Euler method). For the simple wave equation problem, this is - -.. math:: - - \mathcal{I} - \gamma \mathcal{J} = (\begin{array}{cc} - 1 & -\gamma\partial_{||} \\ - -\gamma\partial_{||} & 1 - \end{array} - ) - -This matrix can be block inverted using Schur factorisation [2]_ - -.. math:: +Implicit solvers repeatedly solve linearised systems of the form +:math:`(I - \gamma J)\,x = b`. Preconditioning improves the convergence of +those inner solves by applying a cheap approximate inverse, while +matrix-free or finite-difference Jacobian information controls how the solver +models the couplings in :math:`J`. - (\begin{array}{cc} - {\mathbf{E}} & {\mathbf{U}} \\ - {\mathbf{L}} & {\mathbf{D}} - \end{array})^{-1} - = (\begin{array}{cc} - {\mathbf{I}} & -{\mathbf{E}}^{-1}{\mathbf{U}} \\ - 0 & {\mathbf{I}} - \end{array} - )(\begin{array}{cc} - {\mathbf{E}}^{-1} & 0 \\ - 0 & {\mathbf{P}}_{Schur}^{-1} - \end{array} - )(\begin{array}{cc} - {\mathbf{I}} & 0 \\ - -{\mathbf{L}}{\mathbf{E}}^{-1} & {\mathbf{I}} - \end{array} - ) - -where -:math:`{\mathbf{P}}_{Schur} = {\mathbf{D}} - {\mathbf{L}}{\mathbf{E}}^{-1}{\mathbf{U}}` -Using this, the wave problem becomes: +For most models, the relevant entry points are: -.. math:: - :label: precon - - (\begin{array}{cc} 1 & -\gamma\partial_{||} \\ - -\gamma\partial_{||} & 1 \end{array})^{-1} = (\begin{array}{cc} 1 & \gamma\partial_{||}\\ - 0 & 1 \end{array} )(\begin{array}{cc} 1 & 0 \\ - 0 & (1 -\gamma^2\partial^2_{||})^{-1} \end{array} )(\begin{array}{cc} 1 & 0\\ - \gamma\partial_{||} & 1 \end{array} ) - -The preconditioner is implemented by defining a function of the form - -:: - - int precon(BoutReal t, BoutReal gamma, BoutReal delta) { - ... - } - -which takes as input the current time, the :math:`\gamma` factor -appearing above, and :math:`\delta` which is only important for -constrained problems (not discussed here... yet). The current state of -the system is stored in the state variables (here ``u`` and ``v`` ), -whilst the vector to be preconditioned is stored in the time derivatives -(here ``ddt(u)`` and ``ddt(v)`` ). At the end of the preconditioner the -result should be in the time derivatives. A preconditioner which is just -the identity matrix and so does nothing is therefore:: - - int precon(BoutReal t, BoutReal gamma, BoutReal delta) { - } - -To implement the preconditioner in equation :eq:`precon`, first apply the -rightmost matrix to the given vector: - -.. math:: +- ``solver->setPrecon(...)`` for a user-supplied preconditioner +- ``solver->setJacobian(...)`` for a Jacobian-vector product +- ``matrix_free = false`` and ``use_coloring = true`` for PETSc-backed + finite-difference Jacobians +- ``solver->addJacobianPattern(...)`` when a ``PetscCellOperator`` provides a + better sparsity pattern than the default solver stencil - (\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) = (\begin{array}{cc} - 1 & 0 \\ - \gamma\partial_{||} & 1 - \end{array} - )(\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) +The detailed discussion, worked example, and implementation notes are in +:ref:`sec-preconditioning`. This section keeps the focus on solver options and +setup. -:: - - int precon(BoutReal t, BoutReal gamma, BoutReal delta) { - mesh->communicate(ddt(u)); - //ddt(u) = ddt(u); - ddt(v) = gamma*Grad_par(ddt(u)) + ddt(v); - -note that since the preconditioner is linear, it doesn’t depend on -:math:`u` or :math:`v`. As in the RHS function, since we are taking a -differential of ``ddt(u)``, it first needs to be communicated to -exchange guard cell values. - -The second matrix - -.. math:: - - (\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) \rightarrow (\begin{array}{cc} - 1 & 0 \\ - 0 & (1 - \gamma^2\partial^2_{||})^{-1} - \end{array} - )(\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) - -doesn’t alter :math:`u`, but solves a parabolic equation in the -parallel direction. There is a solver class to do this called -`InvertPar` which solves the equation :math:`(A + B\partial_{||}^2)x = -b` where :math:`A` and :math:`B` are `Field2D` or constants [3]_. In -`PhysicsModel::init` we create one of these solvers:: - - InvertPar *inv; // Parallel inversion class - int init(bool restarting) { - ... - inv = InvertPar::Create(); - inv->setCoefA(1.0); - ... - } - -In the preconditioner we then use this solver to update :math:`v`:: - - inv->setCoefB(-SQ(gamma)); - ddt(v) = inv->solve(ddt(v)); - -which solves -:math:`ddt(v) \rightarrow (1 - \gamma^2\partial_{||}^2)^{-1} ddt(v)`. -The final matrix just updates :math:`u` using this new solution for -:math:`v` - -.. math:: - - (\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) \rightarrow (\begin{array}{cc} - 1 & \gamma\partial_{||} \\ - 0 & 1 - \end{array} - )(\begin{array}{c} - \texttt{ddt(u)} \\ - \texttt{ddt(v)} - \end{array} - ) - -:: - - mesh->communicate(ddt(v)); - ddt(u) = ddt(u) + gamma*Grad_par(ddt(v)); - -Finally, boundary conditions need to be imposed, which should be -consistent with the conditions used in the RHS:: - - ddt(u).applyBoundary("dirichlet"); - ddt(v).applyBoundary("dirichlet"); - -To use the preconditioner, pass the function to the solver in -`PhysicsModel::init`:: - - int init(bool restarting) { - solver->setPrecon(precon); - ... - } - -then in the ``BOUT.inp`` settings file switch on the preconditioner - -.. code-block:: bash +Jacobian-vector function +------------------------ - [solver] - type = cvode # Need CVODE or PETSc - cvode_precon_method = user # Use user-supplied preconditioner - rightprec = false # Use Right preconditioner (default left) +If the physics model can apply the Jacobian to a vector directly, it can provide +that routine through ``solver->setJacobian(...)`` during ``PhysicsModel::init``. +This is mainly useful in matrix-free workflows, where the solver needs +Jacobian-vector products but does not assemble a full sparse matrix. -Jacobian function ------------------ +When a PETSc-backed implicit solver is configured with +``matrix_free = false``, BOUT++ can instead assemble an explicit sparse +Jacobian by finite differences. In that case, the coloring options described +below determine how the sparsity pattern is used, and +``solver->addJacobianPattern(...)`` can add extra block structure from +``PetscCellOperator`` stencils before the colored finite differencing is set up. DAE constraint equations ------------------------ @@ -1458,73 +1245,80 @@ here:\ https://computation.llnl.gov/casc/sundials/support/notes.html). This may in some cases be less efficient. -Implementation internals ------------------------- +Solver interface +---------------- -The solver is the interface between BOUT++ and the time-integration -code such as SUNDIALS. All solvers implement the `Solver` -class interface (see ``src/solver/generic_solver.hxx``). +The solver is the interface between BOUT++ and the time-integration library. +All time integrators implement the ``Solver`` interface, and most user-facing +setup happens during ``PhysicsModel::init`` before ``solver->init()`` is called. -First all the fields which are to be evolved need to be added to the -solver. These are always done in pairs, the first specifying the field, -and the second the time-derivative:: +The main setup methods are: - void add(Field2D &v, Field2D &F_v, const char* name); +.. code-block:: C++ -This is normally called in the `PhysicsModel::init` initialisation routine. -Some solvers (e.g. IDA) can support constraints, which need to be added -in the same way as evolving fields:: + void add(Field2D& v, const std::string& name, + const std::string& description = ""); + void add(Field3D& v, const std::string& name, + const std::string& description = ""); + void add(Vector2D& v, const std::string& name, + const std::string& description = ""); + void add(Vector3D& v, const std::string& name, + const std::string& description = ""); - bool constraints(); - void constraint(Field2D &v, Field2D &C_v, const char* name); + bool constraints(); + void constraint(Field3D& v, Field3D& residual, std::string name); -The ``constraints()`` function tests whether or not the current solver -supports constraints. The format of ``constraint(...)`` is the same as -``add``, except that now the solver will attempt to make ``C_v`` zero. -If ``constraint`` is called when the solver doesn’t support them then an -error should occur. + typedef int (*PhysicsPrecon)(BoutReal t, BoutReal gamma, BoutReal delta); + void setPrecon(PhysicsPrecon f); -If the physics model implements a preconditioner or Jacobian-vector -multiplication routine, these can be passed to the solver during -initialisation:: + typedef int (*Jacobian)(BoutReal t); + void setJacobian(Jacobian j); - typedef int (*PhysicsPrecon)(BoutReal t, BoutReal gamma, BoutReal delta); - void setPrecon(PhysicsPrecon f); // Specify a preconditioner - typedef int (*Jacobian)(BoutReal t); - void setJacobian(Jacobian j); // Specify a Jacobian +#if BOUT_HAS_PETSC + Solver::VarRef getVarRef(std::string_view name) const; + bool addJacobianPattern(const PetscCellOperator& op); + bool addJacobianPattern(const PetscCellOperator& op, Solver::VarRef out_var, + Solver::VarRef in_var); +#endif -If the solver doesn’t support these functions then the calls will just -be ignored. +``add(...)`` registers evolving variables. ``constraint(...)`` registers DAE +constraints for solvers that support them. ``setPrecon(...)`` and +``setJacobian(...)`` register optional callbacks used by implicit solvers. -Once the problem to be solved has been specified, the solver can be -initialised using:: +``getVarRef(...)`` resolves the solver's internal variable numbering from the +name used when the variable or component was added. Those references can then be +used with ``addJacobianPattern(...)`` to target one Jacobian block, one row or +column of blocks, or all blocks via ``Solver::VarRef::All()``. - int init(); +For example: -which returns an error code (0 on success). This is currently called in -:doc:`bout++.cxx<../_breathe_autogen/file/bout_09_09_8cxx>`:: +.. code-block:: C++ - if (solver.init()) { - output.write("Failed to initialise solver. Aborting\n"); - return(1); - } + #if BOUT_HAS_PETSC + auto n = solver->getVarRef("n"); + auto T = solver->getVarRef("T"); -which passes the (physics module) RHS function `PhysicsModel::rhs` to the -solver along with the number and size of the output steps. + solver->addJacobianPattern(op); + solver->addJacobianPattern(op, n, T); + #endif -:: +These Jacobian-pattern registrations are designed to be called during model or +component setup, before the solver knows the final matrix size. BOUT++ therefore +queues them and applies them later when a supported PETSc-backed solver creates +its Jacobian structure. + +Behavior depends on the solver: - typedef int (*MonitorFunc)(BoutReal simtime, int iter, int NOUT); - int run(MonitorFunc f); +- ``setPrecon(...)`` and ``setJacobian(...)`` are ignored by solvers that do + not use those callbacks +- ``addJacobianPattern(...)`` is only available in PETSc-enabled builds +- even with PETSc, ``addJacobianPattern(...)`` may return ``false`` if the + chosen solver does not use the PETSc preconditioner/Jacobian path -.. [1] - Taken from a talk by L.Chacon available here - https://bout2011.llnl.gov/pdf/talks/Chacon_bout2011.pdf +Once the problem has been specified, the solver is initialised with:: + + int init(); -.. [2] - See paper https://arxiv.org/abs/1209.2054 for an application to - 2-fluid equations +and then run in the usual way:: -.. [3] This `InvertPar` class can handle cases with closed - field-lines and twist-shift boundary conditions for tokamak - simulations + int run(); diff --git a/src/mesh/petsc_jacobian.cxx b/src/mesh/petsc_jacobian.cxx new file mode 100644 index 0000000000..457c667de5 --- /dev/null +++ b/src/mesh/petsc_jacobian.cxx @@ -0,0 +1,65 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "bout/assert.hxx" +#include "bout/petsc_jacobian.hxx" +#include "bout/petsclib.hxx" + +#include + +void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) { + // Infer nvars from global sizes + PetscInt jfd_global{0}, sub_global{0}; + BOUT_DO_PETSC(MatGetSize(Jfd, &jfd_global, nullptr)); + BOUT_DO_PETSC(MatGetSize(sub, &sub_global, nullptr)); + + ASSERT1(sub_global > 0); + ASSERT1(jfd_global % sub_global == 0); + + const PetscInt nvars = jfd_global / sub_global; + + ASSERT1(out_var >= 0 && out_var < static_cast(nvars)); + ASSERT1(in_var >= 0 && in_var < static_cast(nvars)); + + // Iterate over locally owned rows of sub and insert into Jfd + PetscInt rstart{0}, rend{0}; + BOUT_DO_PETSC(MatGetOwnershipRange(sub, &rstart, &rend)); + + const PetscScalar one = 1.0; + + for (PetscInt sub_row = rstart; sub_row < rend; ++sub_row) { + PetscInt ncols{0}; + const PetscInt* sub_cols{nullptr}; + const PetscScalar* vals{nullptr}; + BOUT_DO_PETSC(MatGetRow(sub, sub_row, &ncols, &sub_cols, &vals)); + + const PetscInt jfd_row = (sub_row * nvars) + out_var; + + for (PetscInt k = 0; k < ncols; ++k) { + const PetscInt jfd_col = (sub_cols[k] * nvars) + in_var; + BOUT_DO_PETSC(MatSetValues(Jfd, 1, &jfd_row, 1, &jfd_col, &one, INSERT_VALUES)); + } + + BOUT_DO_PETSC(MatRestoreRow(sub, sub_row, &ncols, &sub_cols, &vals)); + } +} + +void addOperatorSparsity(Mat Jfd, Mat sub) { + PetscInt jfd_global{0}, sub_global{0}; + MatGetSize(Jfd, &jfd_global, nullptr); + MatGetSize(sub, &sub_global, nullptr); + + ASSERT1(sub_global > 0); + ASSERT1(jfd_global % sub_global == 0); + + const int nvars = static_cast(jfd_global / sub_global); + + for (int out_var = 0; out_var < nvars; ++out_var) { + for (int in_var = 0; in_var < nvars; ++in_var) { + addOperatorSparsity(Jfd, sub, out_var, in_var); + } + } +} + +#endif // BOUT_HAS_PETSC diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx index 573044a0a1..581e02eb24 100644 --- a/src/mesh/petsc_operators.cxx +++ b/src/mesh/petsc_operators.cxx @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -131,6 +132,33 @@ PetscCellMapping::PetscCellMapping(const Field3D& cell_number, local_indices); } +IS PetscCellMapping::makeEvolvingIS() const { + // Collect global PETSc indices in mapOwnedInteriorCells order. + // Reserve the known count up front to avoid reallocation. + std::vector indices; + indices.reserve(static_cast(evolving_region.size())); + + mapOwnedInteriorCells( + [&](PetscInt row, const Ind3D& /*i*/, int /*stored*/) { indices.push_back(row); }); + + IS is; + BOUT_DO_PETSC(ISCreateGeneral(BoutComm::get(), static_cast(indices.size()), + indices.data(), PETSC_COPY_VALUES, &is)); + return is; +} + +Mat PetscCellMapping::extractEvolvingSubmatrix( + const PetscOperator& op) const { + IS is = makeEvolvingIS(); + + Mat sub; + BOUT_DO_PETSC(MatCreateSubMatrix(op.raw(), is, is, MAT_INITIAL_MATRIX, &sub)); + + BOUT_DO_PETSC(ISDestroy(&is)); + + return sub; +} + PetscLegMapping::PetscLegMapping(int total_legs, std::vector local_leg_indices) { std::sort(local_leg_indices.begin(), local_leg_indices.end()); local_leg_indices.erase(std::unique(local_leg_indices.begin(), local_leg_indices.end()), diff --git a/src/solver/impls/cvode/cvode.cxx b/src/solver/impls/cvode/cvode.cxx index 2cdeeed9c3..2edac8e292 100644 --- a/src/solver/impls/cvode/cvode.cxx +++ b/src/solver/impls/cvode/cvode.cxx @@ -491,6 +491,7 @@ int CvodeSolver::init() { Field3D index = globalIndex(0); PetscCall(petsc_preconditioner.createJacobianPattern( index, *options, local_N, n2Dvars(), n3Dvars(), BoutComm::get())); + applyQueuedJacobianPatterns(petsc_preconditioner.jacobian()); PetscCall( petsc_preconditioner.updateColoring(CvodeSolver::petscFormFunction, this)); PetscCall(MatFDColoringSetF(petsc_preconditioner.coloring(), petsc_f)); diff --git a/src/solver/impls/cvode/cvode.hxx b/src/solver/impls/cvode/cvode.hxx index c242e4128a..bdf601cfd2 100644 --- a/src/solver/impls/cvode/cvode.hxx +++ b/src/solver/impls/cvode/cvode.hxx @@ -44,10 +44,10 @@ RegisterUnavailableSolver #include "../../sundials_nvector_interface.hxx" #include "bout/bout_types.hxx" -#include "bout/region.hxx" #include "bout/sundials_backports.hxx" #if BOUT_HAS_PETSC +#include "bout/petsc_operators.hxx" #include "bout/petsc_preconditioner.hxx" #include "bout/petsclib.hxx" @@ -79,12 +79,19 @@ class CvodeSolver : public Solver { public: explicit CvodeSolver(Options* opts = nullptr); ~CvodeSolver() override; + using Solver::addJacobianPattern; BoutReal getCurrentTimestep() override { return hcur; } int init() override; int run() override; BoutReal run(BoutReal tout); +#if BOUT_HAS_PETSC + bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) override { + return queueJacobianPattern(op, out_var, in_var); + } +#endif void resetInternalFields() override; diff --git a/src/solver/impls/petsc/petsc.cxx b/src/solver/impls/petsc/petsc.cxx index d123b56508..e6fbe95249 100644 --- a/src/solver/impls/petsc/petsc.cxx +++ b/src/solver/impls/petsc/petsc.cxx @@ -519,6 +519,7 @@ int PetscSolver::init() { Field3D index = globalIndex(0); PetscCall(petsc_preconditioner.createJacobianPattern( index, *options, nlocal, n2Dvars(), n3Dvars(), BoutComm::get())); + applyQueuedJacobianPatterns(petsc_preconditioner.jacobian()); output_progress.write("Creating Jacobian coloring\n"); updateColoring(); } else { diff --git a/src/solver/impls/petsc/petsc.hxx b/src/solver/impls/petsc/petsc.hxx index ab8e5bf9e0..fea8a03411 100644 --- a/src/solver/impls/petsc/petsc.hxx +++ b/src/solver/impls/petsc/petsc.hxx @@ -2,7 +2,7 @@ * Interface to PETSc solver * ************************************************************************** - * Copyright 2010 - 2025 BOUT++ contributors + * Copyright 2010 - 2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov * @@ -27,6 +27,7 @@ #define BOUT_PETSC_SOLVER_H #include "bout/build_defines.hxx" +#include "bout/petsc_operators.hxx" #include "bout/solver.hxx" #if not BOUT_HAS_PETSC @@ -40,6 +41,7 @@ RegisterUnavailableSolver class PetscSolver; +#include #include #include #include @@ -53,7 +55,7 @@ class PetscSolver; #include #include -#include +#include namespace { RegisterSolver registersolverpetsc("petsc"); @@ -63,9 +65,14 @@ class PetscSolver : public Solver { public: PetscSolver(Options* opts = nullptr); ~PetscSolver(); + using Solver::addJacobianPattern; int init() override; int run() override; + bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) override { + return queueJacobianPattern(op, out_var, in_var); + } // These functions used internally (but need to be public) diff --git a/src/solver/impls/snes/snes.cxx b/src/solver/impls/snes/snes.cxx index 07d16652f2..415eb77a09 100644 --- a/src/solver/impls/snes/snes.cxx +++ b/src/solver/impls/snes/snes.cxx @@ -78,6 +78,7 @@ PetscErrorCode SNESSolver::FDJinitialise() { Field3D index = globalIndex(0); PetscCall(petsc_preconditioner.createJacobianPattern( index, *options, nlocal, n2Dvars(), n3Dvars(), BoutComm::get())); + applyQueuedJacobianPatterns(petsc_preconditioner.jacobian()); output_progress.write("Creating Jacobian coloring\n"); PetscCall(petsc_preconditioner.updateColoring(FormFunctionForColoring, this)); diff --git a/src/solver/impls/snes/snes.hxx b/src/solver/impls/snes/snes.hxx index 6cdaac082c..4c749f75e0 100644 --- a/src/solver/impls/snes/snes.hxx +++ b/src/solver/impls/snes/snes.hxx @@ -78,10 +78,16 @@ class SNESSolver : public Solver { public: explicit SNESSolver(Options* opts = nullptr); ~SNESSolver() override = default; + using Solver::addJacobianPattern; int init() override; int run() override; + bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) override { + return queueJacobianPattern(op, out_var, in_var); + } + /// Nonlinear function. This is called by PETSc SNES object /// via a static C-style function. For implicit /// time integration this function calculates: diff --git a/src/solver/solver.cxx b/src/solver/solver.cxx index 1f7b654427..b18f30c045 100644 --- a/src/solver/solver.cxx +++ b/src/solver/solver.cxx @@ -47,8 +47,15 @@ #include "bout/vector2d.hxx" #include "bout/vector3d.hxx" +#if BOUT_HAS_PETSC +#include "bout/petsc_jacobian.hxx" +#include "bout/petsc_operators.hxx" +#include "bout/petsclib.hxx" +#endif + #include +#include #include #include #include @@ -350,6 +357,111 @@ void Solver::add(Vector3D& v, const std::string& name, const std::string& descri v3d.emplace_back(std::move(d)); } +Solver::VarRef Solver::getVarRef(std::string_view name) const { + int index = 0; + for (const auto& field : f2d) { + if (field.name == name) { + return VarRef(index); + } + ++index; + } + for (const auto& field : f3d) { + if (field.name == name) { + return VarRef(index); + } + ++index; + } + return VarRef::Invalid(); +} + +#if BOUT_HAS_PETSC +bool Solver::addJacobianPattern(const PetscCellOperator& op) { + return addJacobianPattern(op, VarRef::All(), VarRef::All()); +} + +bool Solver::addJacobianPattern(const PetscCellOperator& UNUSED(op), + VarRef UNUSED(out_var), VarRef UNUSED(in_var)) { + return false; +} + +bool Solver::queueJacobianPattern(const PetscCellOperator& op, VarRef out_var, + VarRef in_var) { + if (initialised || !out_var.isValid() || !in_var.isValid()) { + return false; + } + + auto out_mapping = + std::dynamic_pointer_cast(op.getOutMapping()); + auto in_mapping = std::dynamic_pointer_cast(op.getInMapping()); + if (!out_mapping || !in_mapping || (out_mapping != in_mapping)) { + return false; + } + + DeferredJacobianPattern pattern; + *pattern.submatrix = out_mapping->extractEvolvingSubmatrix(op); + pattern.out_var = out_var; + pattern.in_var = in_var; + deferred_jacobian_patterns.emplace_back(std::move(pattern)); + return true; +} + +bool Solver::canApplyQueuedJacobianPatterns() const { + if (deferred_jacobian_patterns.empty()) { + return true; + } + + if (n2Dvars() != 0 || n3Dvars() == 0) { + return false; + } + + const auto evolve_boundary = [](const auto& field) { return field.evolve_bndry; }; + return std::none_of(f2d.begin(), f2d.end(), evolve_boundary) + && std::none_of(f3d.begin(), f3d.end(), evolve_boundary); +} + +void Solver::applyQueuedJacobianPatterns(Mat Jfd) const { + if (deferred_jacobian_patterns.empty()) { + return; + } + + if (!canApplyQueuedJacobianPatterns()) { + throw BoutException( + "Queued Jacobian patterns require a uniform evolving-cell layout with only " + "Field3D variables and no evolving boundary cells"); + } + + const int nvars = n2Dvars() + n3Dvars(); + ASSERT1(nvars > 0); + + BOUT_DO_PETSC(MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); + + for (const auto& pattern : deferred_jacobian_patterns) { + if (pattern.out_var.isAll() && pattern.in_var.isAll()) { + addOperatorSparsity(Jfd, *pattern.submatrix); + continue; + } + if (pattern.out_var.isAll()) { + for (int out_var = 0; out_var < nvars; ++out_var) { + addOperatorSparsity(Jfd, *pattern.submatrix, out_var, pattern.in_var.index()); + } + continue; + } + if (pattern.in_var.isAll()) { + for (int in_var = 0; in_var < nvars; ++in_var) { + addOperatorSparsity(Jfd, *pattern.submatrix, pattern.out_var.index(), in_var); + } + continue; + } + addOperatorSparsity(Jfd, *pattern.submatrix, pattern.out_var.index(), + pattern.in_var.index()); + } + + BOUT_DO_PETSC(MatAssemblyBegin(Jfd, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(Jfd, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE)); +} +#endif + /************************************************************************** * Constraints **************************************************************************/ diff --git a/tests/integrated/CMakeLists.txt b/tests/integrated/CMakeLists.txt index d87bd329f0..5272ea3378 100644 --- a/tests/integrated/CMakeLists.txt +++ b/tests/integrated/CMakeLists.txt @@ -34,6 +34,7 @@ add_subdirectory(test-naulin-laplace) add_subdirectory(test-options-netcdf) add_subdirectory(test-petsc_laplace) add_subdirectory(test-petsc_laplace_MAST-grid) +add_subdirectory(test-petsc-ordering) add_subdirectory(test-petsc-operators) add_subdirectory(test-restart-io) add_subdirectory(test-restarting) diff --git a/tests/integrated/test-petsc-ordering/CMakeLists.txt b/tests/integrated/test-petsc-ordering/CMakeLists.txt new file mode 100644 index 0000000000..4366c38657 --- /dev/null +++ b/tests/integrated/test-petsc-ordering/CMakeLists.txt @@ -0,0 +1,7 @@ +bout_add_integrated_test( + test-petsc-ordering + SOURCES test_petsc_ordering.cxx + REQUIRES BOUT_HAS_PETSC + USE_RUNTEST USE_DATA_BOUT_INP + PROCESSORS 4 +) diff --git a/tests/integrated/test-petsc-ordering/data/BOUT.inp b/tests/integrated/test-petsc-ordering/data/BOUT.inp new file mode 100644 index 0000000000..9b51b9f296 --- /dev/null +++ b/tests/integrated/test-petsc-ordering/data/BOUT.inp @@ -0,0 +1,14 @@ +# Input file for test_petsc_ordering integrated test. +# Uses a small grid large enough to exercise MPI decomposition on 4 ranks +# (nx=10 gives 2 interior x-points per rank when NXPE=4; ny and nz similarly). + +[mesh] +nx = 10 # 8 interior + 2 guard (mxg=1 each side) +ny = 8 +nz = 4 + +ixseps1 = nx +ixseps2 = nx + +[output] +enabled = false # Suppress data file output; we only need stdout diff --git a/tests/integrated/test-petsc-ordering/runtest b/tests/integrated/test-petsc-ordering/runtest new file mode 100755 index 0000000000..26540d5338 --- /dev/null +++ b/tests/integrated/test-petsc-ordering/runtest @@ -0,0 +1,154 @@ +#!/usr/bin/env python3 +# requires: petsc +# cores: 4 +""" +Integrated test: Ordering equivalence between PetscCellMapping +and the SNES solver's globalIndex traversal. + +Runs test_petsc_ordering with 1, 2, and 4 MPI ranks and checks that: + - The executable exits with status 0. + - The output contains "ordering_check=PASS". + - No "MISMATCH" or "SHIFT_NOT_CONSTANT" lines appear. +""" + +import argparse +import re +import subprocess +import sys + +from boututils.run_wrapper import build_and_log + + +def parse_args(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--executable", + default="./test_petsc_ordering", + help="Path to the test executable", + ) + parser.add_argument( + "--mpirun", default="mpirun", help="MPI launcher (mpirun, srun, ...)" + ) + parser.add_argument( + "--nprocs", + type=int, + nargs="+", + default=[1, 2, 4], + help="List of MPI rank counts to test", + ) + parser.add_argument( + "--timeout", type=int, default=120, help="Per-run timeout in seconds" + ) + return parser.parse_args() + + +def run_case(mpirun, executable, nproc, nxpe, timeout): + """Launch one run and return (stdout+stderr, returncode).""" + cmd = [mpirun, "-n", str(nproc), executable, f"nxpe={nxpe}"] + try: + result = subprocess.run( + cmd, + stdout=subprocess.PIPE, + stderr=subprocess.STDOUT, + timeout=timeout, + text=True, + ) + return result.stdout, result.returncode + except subprocess.TimeoutExpired: + return f"TIMED OUT after {timeout}s\n", -1 + except FileNotFoundError as e: + return f"Could not launch: {e}\n", -1 + + +def check_output(output, nproc): + """ + Inspect the combined output for pass/fail markers. + Returns a list of failure strings (empty means pass). + """ + failures = [] + + # Must contain the summary pass marker + if "ordering_check=PASS" not in output: + if "ordering_check=FAIL" in output: + failures.append("ordering_check=FAIL found in output") + else: + failures.append("ordering_check marker not found in output") + + # Must not contain any mismatch lines + mismatch_lines = [ + line for line in output.splitlines() if line.startswith("MISMATCH") + ] + if mismatch_lines: + failures.append( + f"{len(mismatch_lines)} MISMATCH line(s) found:\n" + + "\n".join(f" {line}" for line in mismatch_lines) + ) + + # Must not contain any shift-not-constant lines + shift_lines = [ + line for line in output.splitlines() if line.startswith("SHIFT_NOT_CONSTANT") + ] + if shift_lines: + failures.append( + f"{len(shift_lines)} SHIFT_NOT_CONSTANT line(s) found:\n" + + "\n".join(f" {line}" for line in shift_lines) + ) + + # Extract and report the numeric summary for informational purposes + for marker in ("total_mismatches", "total_shift_failures"): + m = re.search(rf"{marker}=(\d+)", output) + if m: + count = int(m.group(1)) + if count != 0: + failures.append(f"{marker}={count} (expected 0)") + + return failures + + +def main(): + args = parse_args() + + build_and_log("PETSc ordering test") + + overall_pass = True + + for nproc in args.nprocs: + for nxpe in [1, 2, 4]: + if nxpe > nproc: + break + print(f"\n{'=' * 60}") + print(f"Running with {nproc} MPI rank(s), nxpe={nxpe}") + print(f"{'=' * 60}") + + output, returncode = run_case( + args.mpirun, args.executable, nproc, nxpe, args.timeout + ) + + case_failures = [] + + if returncode != 0: + # Note: MPI task can exit non-zero for reasons not connected to test + print(f"Warning: Non-zero exit code: {returncode}") + + case_failures.extend(check_output(output, nproc)) + + if case_failures: + print(output) # Only print output on failure + overall_pass = False + print(f"FAIL (nproc={nproc}, nxpe={nxpe}):") + for f in case_failures: + print(f" - {f}") + else: + print(f"PASS (nproc={nproc}, nxpe={nxpe})") + + print(f"\n{'=' * 60}") + if overall_pass: + print("ALL CASES PASSED") + return 0 + else: + print("ONE OR MORE CASES FAILED") + return 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx new file mode 100644 index 0000000000..6a0fa25d0b --- /dev/null +++ b/tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx @@ -0,0 +1,221 @@ +/// Integrated test for ordering equivalence between PetscCellMapping +/// and the SNES solver's globalIndex traversal. +/// +/// For every locally owned evolving interior cell this program checks: +/// +/// 1. The local offset into mapOwnedInteriorCells equals the local SNES +/// index from globalIndex (ordering equivalence). +/// +/// 2. The shift (rowStart() - Istart) is the same constant for every cell +/// on this rank (scalar-shift property). +/// +/// Results are written to the BOUT++ output so the runtest script can parse +/// them. The program exits with a non-zero status if any check fails. + +#include + +#if !BOUT_HAS_PETSC +// Without PETSc there is nothing to test. Exit cleanly so CTest skips rather +// than fails. +#include +int main() { return EXIT_SUCCESS; } +#else + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include + +#include +#include +#include +#include + +int main(int argc, char** argv) { + int total_mismatches = 0; + int total_shift_fail = 0; + + // Initialise BOUT++ (reads BOUT.inp, sets up mesh, MPI, PETSc, etc.) + BoutInitialise(argc, argv); + { + Mesh* mesh = bout::globals::mesh; + + // ── Build cell_number field ──────────────────────────────────────────────── + // Assign stored indices to interior cells exactly as the Python weights + // module does: consecutive integers over the evolving region. + Field3D cell_number{-1.0, mesh}; + { + const int ngy = mesh->GlobalNyNoBoundaries; + const int ngz = mesh->GlobalNzNoBoundaries; + + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + const int global_x = + mesh->getGlobalXIndexNoBoundaries(x); // Exclude boundary cells + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + const int global_y = mesh->getGlobalYIndexNoBoundaries(y); + for (int z = 0; z < mesh->LocalNz; ++z) { + const int global_z = mesh->getGlobalZIndexNoBoundaries(z); + cell_number(x, y, z) = ((global_x * ngy) + global_y) * ngz + global_z; + } + } + } + // Communicate so guard cells are filled (needed by some mesh operations) + mesh->communicate(cell_number); + } + + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + // Count total evolving cells across all ranks + int n_local_evolving = 0; + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + if (cell_number[i] >= 0) { + ++n_local_evolving; + } + } + int n_global_evolving = 0; + MPI_Allreduce(&n_local_evolving, &n_global_evolving, 1, MPI_INT, MPI_SUM, + BoutComm::get()); + + // ── Build PetscCellMapping ───────────────────────────────────────────────── + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, n_global_evolving); + + // ── Build globalIndex field (replicates FDJinitialise logic) ────────────── + // globalIndex(0) gives 0-based local indices over RGN_NOBNDRY. + // After shifting by Istart these become the SNES global indices. + Field3D snes_index{-1.0, mesh}; + { + int local_idx = 0; + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + if (cell_number[i] >= 0) { + snes_index[i] = local_idx++; + } + } + } + + // Determine SNES global offset for this rank by constructing a dummy PETSc + // matrix of size n_global_evolving and reading its ownership range. + Mat dummy{nullptr}; + MatCreate(BoutComm::get(), &dummy); + MatSetSizes(dummy, n_local_evolving, n_local_evolving, PETSC_DETERMINE, + PETSC_DETERMINE); + MatSetType(dummy, MATMPIAIJ); + MatSetUp(dummy); + PetscInt Istart, Iend; + MatGetOwnershipRange(dummy, &Istart, &Iend); + MatDestroy(&dummy); + + // ── Compare traversal orderings ─────────────────────────────────────────── + // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that + // the local offset produced by each is identical for every cell. + + struct Mismatch { + int x, y, z; + PetscInt mapping_local; // local offset from mapOwnedInteriorCells + int snes_local; // local offset from globalIndex + }; + std::vector mismatches; + + // Collect (Ind3D -> mapping local offset) from mapOwnedInteriorCells + // so we can look up by field index. + // mapping_local[Ind3D linear] = local offset (0-based on this rank) + const PetscInt rstart = mapping->rowStart(); + std::vector mapping_local_by_ind( + static_cast(mesh->LocalNx * mesh->LocalNy * mesh->LocalNz), -1); + { + PetscInt local_offset = 0; + mapping->mapOwnedInteriorCells( + [&](PetscInt petsc_row, const Ind3D& i, int /*stored*/) { + mapping_local_by_ind[i.ind] = static_cast(petsc_row - rstart); + ++local_offset; + }); + } + + // Walk RGN_NOBNDRY and compare + bool shift_initialised = false; + PetscInt expected_shift = 0; + bool shift_constant = true; + + int snes_local_idx = 0; + for (const auto& i : mesh->getRegion3D("RGN_NOBNDRY")) { + if (cell_number[i] < 0) { + continue; + } + + const int map_local = mapping_local_by_ind[i.ind]; + if (map_local < 0) { + // Cell is in RGN_NOBNDRY but not visited by mapOwnedInteriorCells + mismatches.push_back({i.x(), i.y(), i.z(), -1, snes_local_idx}); + ++snes_local_idx; + continue; + } + + // Check ordering equivalence: local offsets must match + if (map_local != snes_local_idx) { + mismatches.push_back( + {i.x(), i.y(), i.z(), static_cast(map_local), snes_local_idx}); + } + + // Check scalar-shift property: petsc_global - snes_global must be constant + const PetscInt petsc_global = rstart + map_local; + const PetscInt snes_global = Istart + snes_local_idx; + const PetscInt shift = petsc_global - snes_global; + if (!shift_initialised) { + expected_shift = shift; + shift_initialised = true; + } else if (shift != expected_shift) { + shift_constant = false; + } + + ++snes_local_idx; + } + + // ── Gather results across ranks ─────────────────────────────────────────── + const int rank_mismatches = static_cast(mismatches.size()); + MPI_Allreduce(&rank_mismatches, &total_mismatches, 1, MPI_INT, MPI_SUM, + BoutComm::get()); + + int rank_shift_fail = shift_constant ? 0 : 1; + MPI_Allreduce(&rank_shift_fail, &total_shift_fail, 1, MPI_INT, MPI_SUM, + BoutComm::get()); + + // ── Report ──────────────────────────────────────────────────────────────── + const int my_rank = BoutComm::rank(); + + // Each rank writes its own diagnostics; rank 0 also writes the summary. + for (const auto& m : mismatches) { + output.write("MISMATCH rank={} cell=({},{},{}) mapping_local={} snes_local={}\n", + my_rank, m.x, m.y, m.z, m.mapping_local, m.snes_local); + } + if (!shift_constant) { + output.write("SHIFT_NOT_CONSTANT rank={} expected_shift={}\n", my_rank, + expected_shift); + } + + if (my_rank == 0) { + output.write("total_mismatches={}\n", total_mismatches); + output.write("total_shift_failures={}\n", total_shift_fail); + if (total_mismatches == 0 && total_shift_fail == 0) { + output.write("ordering_check=PASS\n"); + } else { + output.write("ordering_check=FAIL\n"); + } + } + } + BoutFinalise(); + + return (total_mismatches == 0 && total_shift_fail == 0) ? EXIT_SUCCESS : EXIT_FAILURE; +} + +#endif // BOUT_HAS_PETSC diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index e8c4da31ca..98d3f4c816 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -100,6 +100,9 @@ set(serial_tests_source ./mesh/test_mesh.cxx ./mesh/test_paralleltransform.cxx ./mesh/test_petsc_operators.cxx + ./mesh/test_petsc_operators_make_evolving_is.cxx + ./mesh/test_petsc_operators_extract_evolving_submatrix.cxx + ./mesh/test_petsc_jacobian.cxx ./solver/test_fakesolver.cxx ./solver/test_fakesolver.hxx ./solver/test_nvector.cxx diff --git a/tests/unit/mesh/test_petsc_jacobian.cxx b/tests/unit/mesh/test_petsc_jacobian.cxx new file mode 100644 index 0000000000..ac2e8e5d08 --- /dev/null +++ b/tests/unit/mesh/test_petsc_jacobian.cxx @@ -0,0 +1,638 @@ +// ============================================================================ +// addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) +// +// Concept: the nonzero pattern of `sub` is scattered into the variable block +// (out_var, in_var) of the pre-allocated Jacobian `Jfd`, with the correct +// index stride. +// +// `sub` is a square matrix of size n_evolving x n_evolving. +// `Jfd` is a square matrix of size (n_evolving * nvars) x (n_evolving * nvars). +// nvars is inferred internally as Jfd_global / sub_global. +// +// A nonzero at (r, c) in sub produces an entry at +// (r * nvars + out_var, c * nvars + in_var) +// in Jfd. +// +// ============================================================================ + +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +#include "bout/petsc_jacobian.hxx" + +#include + +#include +#include +#include +#include + +// ── Helpers ────────────────────────────────────────────────────────────────── + +namespace { + +// Build and assemble a square MPIAIJ matrix of global size n x n, +// with the given (global_row, global_col) nonzero positions (value = 1.0). +// nlocal is the number of locally owned rows on this rank. +// rstart is the first globally owned row. +Mat buildPatternMat(PetscInt n, PetscInt nlocal, PetscInt rstart, + const std::vector>& entries) { + Mat A; + MatCreate(BoutComm::get(), &A); + MatSetSizes(A, nlocal, nlocal, n, n); + MatSetType(A, MATMPIAIJ); + // Generous preallocation: allow up to n nonzeros per row + MatMPIAIJSetPreallocation(A, n, nullptr, n, nullptr); + MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + const PetscScalar one = 1.0; + for (const auto& [r, c] : entries) { + MatSetValues(A, 1, &r, 1, &c, &one, INSERT_VALUES); + } + MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); + return A; +} + +// Collect all (global_row, global_col) nonzero positions in a matrix, +// iterating only over locally owned rows. +std::vector> matNonzeroPositions(Mat A) { + PetscInt rstart, rend; + MatGetOwnershipRange(A, &rstart, &rend); + std::vector> result; + for (PetscInt row = rstart; row < rend; ++row) { + PetscInt ncols; + const PetscInt* cols; + const PetscScalar* vals; + MatGetRow(A, row, &ncols, &cols, &vals); + for (PetscInt k = 0; k < ncols; ++k) { + result.emplace_back(row, cols[k]); + } + MatRestoreRow(A, row, &ncols, &cols, &vals); + } + return result; +} + +// Build a Jfd of size (n_sub * nvars) x (n_sub * nvars), with nlocal_sub +// rows per variable per rank, and rstart_sub the sub-matrix row offset. +// Returns the matrix plus the Jfd local size and row start. +struct JfdInfo { + Mat jfd; + PetscInt nlocal; + PetscInt rstart; +}; + +JfdInfo buildEmptyJfd(PetscInt n_sub, PetscInt nlocal_sub, PetscInt rstart_sub, + int nvars) { + const PetscInt n = n_sub * nvars; + const PetscInt nlocal = nlocal_sub * nvars; + + Mat jfd; + MatCreate(BoutComm::get(), &jfd); + MatSetSizes(jfd, nlocal, nlocal, n, n); + MatSetType(jfd, MATMPIAIJ); + MatMPIAIJSetPreallocation(jfd, n, nullptr, n, nullptr); + MatSetOption(jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE); + + PetscInt jfd_rstart, jfd_rend; + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + MatGetOwnershipRange(jfd, &jfd_rstart, &jfd_rend); + return {jfd, nlocal, jfd_rstart}; +} + +// Call addOperatorSparsity and assemble Jfd. +void applyAndAssemble(Mat jfd, Mat sub, int out_var, int in_var) { + addOperatorSparsity(jfd, sub, out_var, in_var); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); +} + +} // namespace + +// ── Fixture ────────────────────────────────────────────────────────────────── + +using PetscAddSparsityTest = FakeMeshFixture; + +// ============================================================================ +// single_variable_diagonal_sub_produces_diagonal_jfd +// +// With nvars=1 a diagonal sub must produce a diagonal Jfd. +// The simplest case: no variable interleaving. +// ============================================================================ +TEST_F(PetscAddSparsityTest, single_variable_diagonal_sub_produces_diagonal_jfd) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; // single rank + const PetscInt rstart_sub = 0; + const int nvars = 1; + + // Diagonal sub: (0,0), (1,1), (2,2) + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, rstart_sub, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = + buildEmptyJfd(n_sub, nlocal_sub, rstart_sub, nvars); + + applyAndAssemble(jfd, sub, /*out_var=*/0, /*in_var=*/0); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(3U, nzs.size()); + for (const auto& [r, c] : nzs) { + EXPECT_EQ(r, c) << "Expected diagonal entry, got (" << r << "," << c << ")"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// correct_block_is_filled_diagonal_coupling +// +// With nvars=3, out_var=in_var=1 (self-coupling of variable 1), +// a diagonal sub must populate only the (1,1) variable block of Jfd. +// ============================================================================ +TEST_F(PetscAddSparsityTest, correct_block_is_filled_diagonal_coupling) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + const int out_var = 1; + const int in_var = 1; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(3U, nzs.size()); + + for (const auto& [r, c] : nzs) { + // Row must be in the out_var=1 slot: r % nvars == 1 + EXPECT_EQ(out_var, r % nvars) + << "Entry (" << r << "," << c << ") is in wrong row variable"; + // Col must be in the in_var=1 slot: c % nvars == 1 + EXPECT_EQ(in_var, c % nvars) + << "Entry (" << r << "," << c << ") is in wrong col variable"; + // Cell index must match: r / nvars == c / nvars for diagonal sub + EXPECT_EQ(r / nvars, c / nvars) + << "Entry (" << r << "," << c << ") cell indices don't match"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// correct_block_is_filled_off_diagonal_coupling +// +// With nvars=3, out_var=0, in_var=2 (cross-variable coupling), +// all Jfd entries must be in the (0,2) variable block. +// ============================================================================ +TEST_F(PetscAddSparsityTest, correct_block_is_filled_off_diagonal_coupling) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + const int out_var = 0; + const int in_var = 2; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(3U, nzs.size()); + + for (const auto& [r, c] : nzs) { + EXPECT_EQ(out_var, r % nvars) + << "Entry (" << r << "," << c << ") is in wrong row variable"; + EXPECT_EQ(in_var, c % nvars) + << "Entry (" << r << "," << c << ") is in wrong col variable"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// stride_is_correct_for_nvars_2 +// +// With nvars=2 and a sub containing entries (0,0) and (0,1), the Jfd +// entry positions must follow the stride formula: +// Jfd row = sub_row * nvars + out_var +// Jfd col = sub_col * nvars + in_var +// ============================================================================ +TEST_F(PetscAddSparsityTest, stride_is_correct_for_nvars_2) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + const int out_var = 0; + const int in_var = 1; + + // Sub entries: (0,0), (1,2), (2,1) — arbitrary non-diagonal pattern + const std::vector> sub_entries = {{0, 0}, {1, 2}, {2, 1}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + // Build expected Jfd positions using the stride formula + std::vector> expected; + for (const auto& [sr, sc] : sub_entries) { + expected.emplace_back((sr * nvars) + out_var, (sc * nvars) + in_var); + } + std::sort(expected.begin(), expected.end()); + + auto nzs = matNonzeroPositions(jfd); + std::sort(nzs.begin(), nzs.end()); + + ASSERT_EQ(expected.size(), nzs.size()); + for (std::size_t k = 0; k < expected.size(); ++k) { + EXPECT_EQ(expected[k].first, nzs[k].first) << "Row mismatch at entry " << k; + EXPECT_EQ(expected[k].second, nzs[k].second) << "Col mismatch at entry " << k; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// no_other_blocks_written +// +// With nvars=3 and only (out_var=1, in_var=1) populated, the other eight +// variable blocks must remain empty. +// ============================================================================ +TEST_F(PetscAddSparsityTest, no_other_blocks_written) { + const PetscInt n_sub = 4; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + const int out_var = 1; + const int in_var = 1; + + const std::vector> sub_entries = { + {0, 0}, {1, 1}, {2, 2}, {3, 3}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, out_var, in_var); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()) << "Wrong number of entries in Jfd"; + for (const auto& [r, c] : nzs) { + EXPECT_EQ(out_var, r % nvars) << "Unexpected entry in row-variable " << r % nvars + << " at (" << r << "," << c << ")"; + EXPECT_EQ(in_var, c % nvars) << "Unexpected entry in col-variable " << c % nvars + << " at (" << r << "," << c << ")"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// multiple_calls_union_patterns +// +// Two calls with different (out_var, in_var) pairs must union their patterns: +// the total nonzero count equals the sum of each operator's count, and the +// entries from each call are in the correct block. +// ============================================================================ +TEST_F(PetscAddSparsityTest, multiple_calls_union_patterns) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + // Sub A: diagonal + const std::vector> entries_a = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub_a = buildPatternMat(n_sub, nlocal_sub, 0, entries_a); + + // Sub B: one off-diagonal entry per row + const std::vector> entries_b = {{0, 1}, {1, 2}, {2, 0}}; + Mat sub_b = buildPatternMat(n_sub, nlocal_sub, 0, entries_b); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + // Insert both without assembling between calls + addOperatorSparsity(jfd, sub_a, /*out_var=*/0, /*in_var=*/0); + addOperatorSparsity(jfd, sub_b, /*out_var=*/1, /*in_var=*/1); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + + // Total count: 3 from sub_a + 3 from sub_b = 6 + EXPECT_EQ(6U, nzs.size()); + + // Check each entry is in the correct variable block + for (const auto& [r, c] : nzs) { + const int rv = static_cast(r % nvars); + const int cv = static_cast(c % nvars); + const bool in_a_block = (rv == 0 && cv == 0); + const bool in_b_block = (rv == 1 && cv == 1); + EXPECT_TRUE(in_a_block || in_b_block) + << "Entry (" << r << "," << c << ") is in unexpected variable block (" << rv + << "," << cv << ")"; + } + + MatDestroy(&sub_a); + MatDestroy(&sub_b); + MatDestroy(&jfd); +} + +// ============================================================================ +// empty_sub_adds_no_entries +// +// A zero-pattern sub must not modify Jfd. +// ============================================================================ +TEST_F(PetscAddSparsityTest, empty_sub_adds_no_entries) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, {}); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, 0, 0); + + const auto nzs = matNonzeroPositions(jfd); + EXPECT_EQ(0U, nzs.size()); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// nonzero_count_equals_sub_nonzero_count +// +// For a single call the total Jfd nonzero count must equal the sub nonzero +// count exactly — no extras, no missing. +// ============================================================================ +TEST_F(PetscAddSparsityTest, nonzero_count_equals_sub_nonzero_count) { + const PetscInt n_sub = 4; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + + // Arbitrary non-trivial pattern + const std::vector> sub_entries = { + {0, 0}, {0, 1}, {1, 1}, {1, 3}, {2, 2}, {3, 0}, {3, 3}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + applyAndAssemble(jfd, sub, 0, 0); + + MatInfo info; + MatGetInfo(jfd, MAT_GLOBAL_SUM, &info); + EXPECT_EQ(static_cast(sub_entries.size()), + static_cast(info.nz_used)); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// first_and_last_variable_blocks_correct +// +// Explicitly check out_var=0,in_var=0 and out_var=nvars-1,in_var=nvars-1 +// to confirm the stride formula is correct at the boundaries. +// ============================================================================ +TEST_F(PetscAddSparsityTest, first_and_last_variable_blocks_correct) { + const PetscInt n_sub = 2; + const PetscInt nlocal_sub = n_sub; + const int nvars = 4; + + const std::vector> sub_entries = { + {0, 0}, {0, 1}, {1, 0}, {1, 1}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + // Test first block (0,0) + { + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + applyAndAssemble(jfd, sub, 0, 0); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()) << "First block: wrong number of entries"; + for (const auto& [r, c] : nzs) { + EXPECT_EQ(0, r % nvars) << "First block row variable wrong"; + EXPECT_EQ(0, c % nvars) << "First block col variable wrong"; + } + MatDestroy(&jfd); + } + + // Test last block (nvars-1, nvars-1) + { + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + applyAndAssemble(jfd, sub, nvars - 1, nvars - 1); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()) << "Last block: wrong number of entries"; + for (const auto& [r, c] : nzs) { + EXPECT_EQ(nvars - 1, r % nvars) << "Last block row variable wrong"; + EXPECT_EQ(nvars - 1, c % nvars) << "Last block col variable wrong"; + } + MatDestroy(&jfd); + } + + MatDestroy(&sub); +} + +// ============================================================================ +// all_pairs_single_variable +// +// With nvars=1 the all-pairs overload is identical to addOperatorSparsity +// with out_var=0, in_var=0. +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_single_variable) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 1; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + ASSERT_EQ(sub_entries.size(), nzs.size()); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// all_pairs_fills_every_variable_block +// +// With nvars=3 every one of the 9 variable blocks must contain exactly the +// same pattern as sub. The total nonzero count must equal +// nvars * nvars * nnz(sub). +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_fills_every_variable_block) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + + const std::vector> sub_entries = {{0, 0}, {1, 1}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + const std::size_t expected_total = + static_cast(nvars * nvars) * sub_entries.size(); + ASSERT_EQ(expected_total, nzs.size()); + + // Every (out_var, in_var) pair must appear exactly once per sub entry. + // Count how many Jfd entries land in each variable block. + std::vector> block_counts(nvars, std::vector(nvars, 0)); + for (const auto& [r, c] : nzs) { + block_counts[r % nvars][c % nvars]++; + } + for (int ov = 0; ov < nvars; ++ov) { + for (int iv = 0; iv < nvars; ++iv) { + EXPECT_EQ(static_cast(sub_entries.size()), block_counts[ov][iv]) + << "Block (" << ov << "," << iv << ") has wrong entry count"; + } + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// all_pairs_entries_match_single_pair_calls +// +// The all-pairs overload must produce exactly the same Jfd nonzero positions +// as calling the single-pair overload for every (out_var, in_var). +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_entries_match_single_pair_calls) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + // Non-trivial pattern so the stride mapping is exercised + const std::vector> sub_entries = { + {0, 0}, {0, 2}, {1, 1}, {2, 0}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + // Build expected positions by calling the single-pair version for each block + std::vector> expected; + { + auto [jfd_ref, nlocal_ref, rstart_ref] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + for (int ov = 0; ov < nvars; ++ov) { + for (int iv = 0; iv < nvars; ++iv) { + addOperatorSparsity(jfd_ref, sub, ov, iv); + } + } + MatAssemblyBegin(jfd_ref, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd_ref, MAT_FINAL_ASSEMBLY); + expected = matNonzeroPositions(jfd_ref); + MatDestroy(&jfd_ref); + } + + // Build actual positions using the all-pairs overload + std::vector> actual; + { + auto [jfd_act, nlocal_act, rstart_act] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + addOperatorSparsity(jfd_act, sub); + MatAssemblyBegin(jfd_act, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd_act, MAT_FINAL_ASSEMBLY); + actual = matNonzeroPositions(jfd_act); + MatDestroy(&jfd_act); + } + + std::sort(expected.begin(), expected.end()); + std::sort(actual.begin(), actual.end()); + + ASSERT_EQ(expected.size(), actual.size()); + for (std::size_t k = 0; k < expected.size(); ++k) { + EXPECT_EQ(expected[k], actual[k]) << "Mismatch at position " << k; + } + + MatDestroy(&sub); +} + +// ============================================================================ +// all_pairs_empty_sub_fills_nothing +// +// An empty sub must leave Jfd empty regardless of nvars. +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_empty_sub_fills_nothing) { + const PetscInt n_sub = 4; + const PetscInt nlocal_sub = n_sub; + const int nvars = 3; + + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, {}); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + EXPECT_EQ(0U, nzs.size()); + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +// ============================================================================ +// all_pairs_no_entries_outside_expected_positions +// +// Every entry in Jfd must lie at a position consistent with the stride formula +// applied to some (out_var, in_var) pair. Specifically, for each entry +// (r, c) in Jfd there must exist a (sr, sc) in sub such that +// r == sr * nvars + (r % nvars) and c == sc * nvars + (c % nvars). +// ============================================================================ +TEST_F(PetscAddSparsityTest, all_pairs_no_entries_outside_expected_positions) { + const PetscInt n_sub = 3; + const PetscInt nlocal_sub = n_sub; + const int nvars = 2; + + std::vector> sub_entries = { + {0, 1}, {1, 0}, {1, 2}, {2, 2}}; + Mat sub = buildPatternMat(n_sub, nlocal_sub, 0, sub_entries); + + // Pre-compute sub cell pairs for fast lookup + const std::set> sub_set(sub_entries.begin(), + sub_entries.end()); + + auto [jfd, jfd_nlocal, jfd_rstart] = buildEmptyJfd(n_sub, nlocal_sub, 0, nvars); + + addOperatorSparsity(jfd, sub); + MatAssemblyBegin(jfd, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(jfd, MAT_FINAL_ASSEMBLY); + + const auto nzs = matNonzeroPositions(jfd); + for (const auto& [r, c] : nzs) { + const PetscInt sr = r / nvars; + const PetscInt sc = c / nvars; + EXPECT_TRUE(sub_set.count({sr, sc}) > 0) + << "Entry (" << r << "," << c << ") maps to sub cell (" << sr << "," << sc + << ") which is not in sub"; + } + + MatDestroy(&sub); + MatDestroy(&jfd); +} + +#endif // BOUT_HAS_PETSC diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx index 2ba63f4528..853578ba15 100644 --- a/tests/unit/mesh/test_petsc_operators.cxx +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -7,7 +7,6 @@ #include "bout/array.hxx" #include "bout/bout_types.hxx" #include "bout/output.hxx" -#include "bout/output_bout_types.hxx" #include "bout/petsc_operators.hxx" #include "bout/region.hxx" diff --git a/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx b/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx new file mode 100644 index 0000000000..4ef74318f2 --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators_extract_evolving_submatrix.cxx @@ -0,0 +1,489 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#include "bout/array.hxx" +#include "bout/bout_types.hxx" +#include "bout/globals.hxx" +#include "bout/petsc_operators.hxx" +#include "bout/region.hxx" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +#include +#include +#include +#include + +#include + +using bout::globals::mesh; + +// ============================================================================ +// PetscCellMapping::extractEvolvingSubmatrix() +// +// Concept: given a PetscCellOperator with a known nonzero pattern, the +// submatrix returned by extractEvolvingSubmatrix() contains exactly the rows +// and columns that correspond to evolving interior cells, with values +// preserved, and no entries from boundary cells. +// +// The fixture is FakeMeshFixture (LocalNx=5, LocalNy=4, LocalNz=2, mxg=1). +// Interior cells: x in [1,3], y in [1,2], z in [0,1] => 12 evolving cells. +// ============================================================================ + +// ── Helpers ───────────────────────────────────────────────────────────────── + +namespace { + +// Build a mapping where evolving cells have stored indices [0, n_evolving), +// with optional extra boundary cells appended after. +struct MappingWithCells { + std::shared_ptr mapping; + int n_evolving; + // Global PETSc index of each evolving cell, in mapOwnedInteriorCells order + std::vector evolving_petsc; + // Stored indices of boundary cells (xin, xout, yup, ydown) + std::vector boundary_stored; +}; + +MappingWithCells buildMappingWithBoundaries(Mesh* mesh, bool add_xin = false, + bool add_xout = false, bool add_yup = false, + bool add_ydown = false) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + const int n_evolving = next; + + std::vector boundary_stored; + if (add_xin) { + cell_number(0, mesh->ystart, 0) = next; + boundary_stored.push_back(next++); + } + if (add_xout) { + cell_number(mesh->LocalNx - 1, mesh->ystart, 0) = next; + boundary_stored.push_back(next++); + } + if (add_yup) { + forward_cell_number(mesh->xstart, mesh->yend, 0) = next; + boundary_stored.push_back(next++); + } + if (add_ydown) { + backward_cell_number(mesh->xstart, mesh->ystart, 0) = next; + boundary_stored.push_back(next++); + } + + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, next); + + std::vector evolving_petsc; + mapping->mapOwnedInteriorCells( + [&](PetscInt row, const Ind3D&, int) { evolving_petsc.push_back(row); }); + + return {std::move(mapping), n_evolving, std::move(evolving_petsc), + std::move(boundary_stored)}; +} + +// Build a PetscCellOperator from explicit (row, col, value) triples over a +// given mapping. Rows and cols are stored indices (not PETSc global). +// The operator is assembled from a hand-built CSR representation. +PetscCellOperator +buildOperatorFromTriples(std::shared_ptr mapping, + const std::vector>& triples) { + + const int n = mapping->globalSize(); + // Build dense CSR: rows array size n+1 + std::vector row_counts(n, 0); + for (const auto& [r, c, v] : triples) { + row_counts[r]++; + } + // rows (CSR row pointer) + Array rows(n + 1); + rows[0] = 0; + for (int i = 0; i < n; ++i) { + rows[i + 1] = rows[i] + row_counts[i]; + } + // cols and weights + Array cols(static_cast(triples.size())); + Array weights(static_cast(triples.size())); + std::vector fill(n, 0); + for (const auto& [r, c, v] : triples) { + const int pos = rows[r] + fill[r]++; + cols[pos] = c; + weights[pos] = v; + } + + return PetscCellOperator(mapping, mapping, rows, cols, weights); +} + +// Return a flat vector of (global_row, global_col, value) for every nonzero +// in a matrix, using MatGetRow. Rows are iterated over the local ownership +// range only. +struct Nonzero { + PetscInt row, col; + PetscScalar val; +}; + +std::vector matNonzeros(Mat A) { + PetscInt rstart, rend; + MatGetOwnershipRange(A, &rstart, &rend); + std::vector nzs; + for (PetscInt row = rstart; row < rend; ++row) { + PetscInt ncols; + const PetscInt* col_idx; + const PetscScalar* vals; + MatGetRow(A, row, &ncols, &col_idx, &vals); + for (PetscInt k = 0; k < ncols; ++k) { + nzs.push_back({row, col_idx[k], vals[k]}); + } + MatRestoreRow(A, row, &ncols, &col_idx, &vals); + } + return nzs; +} + +} // namespace + +// ── Fixture ───────────────────────────────────────────────────────────────── + +using PetscExtractEvolvingTest = FakeMeshFixture; + +// ============================================================================ +// submatrix_global_size_equals_n_evolving_squared +// +// MatGetSize on the result must report n_evolving × n_evolving. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_global_size_equals_n_evolving_squared) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/true); + + // Diagonal operator over the full cell space (safe: every cell maps to itself) + const auto& local = mapping->localStoredIndices(); + std::vector> triples; + for (int s : local) { + if (s >= 0) { + triples.emplace_back(s, s, 1.0); + } + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + + EXPECT_EQ(n_evolving, sub_rows); + EXPECT_EQ(n_evolving, sub_cols); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_local_size_equals_n_evolving_locally +// +// MatGetLocalSize must report the local evolving count on each rank. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_local_size_equals_n_evolving_locally) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, true, true, true, true); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt local_rows, local_cols; + MatGetLocalSize(sub, &local_rows, &local_cols); + + // On a single rank local == global + EXPECT_EQ(n_evolving, local_rows); + EXPECT_EQ(n_evolving, local_cols); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_values_preserved_for_evolving_block +// +// For entries (r, c) where both r and c are evolving cells, the value in the +// submatrix must equal the value in the original operator. +// +// Strategy: build an operator with known values on evolving-to-evolving +// entries, extract the submatrix, and check each value via MatGetValues. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_values_preserved_for_evolving_block) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/true, + /*yup=*/true, /*ydown=*/true); + + // Insert entries with distinct values between pairs of evolving cells. + // Use stored indices directly (0..n_evolving-1). + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + // Diagonal entry + triples.emplace_back(i, i, static_cast(10 + i)); + // One off-diagonal entry per row (wrap around) + const int j = (i + 1) % n_evolving; + triples.emplace_back(i, j, static_cast(100 + i)); + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + const auto nzs = matNonzeros(sub); + + // Build expected map: submatrix uses 0-based indices within the evolving + // block. The IS preserves mapOwnedInteriorCells order, which matches + // stored index order (0..n_evolving-1) on a single rank. + // Sub row/col i corresponds to evolving stored index i. + PetscInt sub_rstart, sub_rend; + MatGetOwnershipRange(sub, &sub_rstart, &sub_rend); + + for (const auto& nz : nzs) { + // nz.row and nz.col are global indices in the submatrix space. + // On a single rank these equal the stored indices of the evolving cells. + const int stored_row = static_cast(nz.row - sub_rstart); + const int stored_col = static_cast(nz.col); + + // Find the expected value from our triples + BoutReal expected = 0.0; + for (const auto& [r, c, v] : triples) { + if (r == stored_row && c == stored_col) { + expected = v; + break; + } + } + EXPECT_DOUBLE_EQ(expected, nz.val) + << "Wrong value at submatrix (" << nz.row << ", " << nz.col << ")"; + } + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_nonzero_count_matches_evolving_entries +// +// The total nonzero count in the submatrix must equal the number of entries +// in the original operator whose row AND column are both evolving cells. +// Entries coupling to boundary cells must be absent. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_nonzero_count_matches_evolving_entries) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/false, + /*yup=*/true, /*ydown=*/false); + + std::vector> triples; + // Entries entirely within the evolving block + int expected_in_sub = 0; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + ++expected_in_sub; + } + // Entries from an evolving row into a boundary column — must be excluded + for (int bnd : boundary_stored) { + triples.emplace_back(0, bnd, 99.0); // row 0 (evolving) -> boundary col + } + // Entries from a boundary row into an evolving column — must be excluded + for (int bnd : boundary_stored) { + triples.emplace_back(bnd, 0, 88.0); // boundary row -> col 0 (evolving) + } + + const auto op = buildOperatorFromTriples(mapping, triples); + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt nnz; + MatInfo info; + MatGetInfo(sub, MAT_GLOBAL_SUM, &info); + nnz = static_cast(info.nz_used); + + EXPECT_EQ(expected_in_sub, nnz); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_excludes_boundary_columns +// +// No column index in the submatrix should correspond to a boundary cell's +// PETSc index in the original cell space. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_excludes_boundary_columns) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/true, + /*yup=*/true, /*ydown=*/true); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + // Add evolving-row → boundary-column coupling + for (int bnd : boundary_stored) { + triples.emplace_back(0, bnd, 7.0); + } + + const auto op = buildOperatorFromTriples(mapping, triples); + + // Collect the PETSc indices of the boundary cells so we can check against them. + // The submatrix uses a re-indexed column space [0, n_evolving), so boundary + // columns from the full space simply will not appear. + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + + // The boundary-coupling entries were given the sentinel value 7.0. + // If any survive into the submatrix, the extraction has failed to exclude them. + const auto nzs = matNonzeros(sub); + for (const auto& nz : nzs) { + EXPECT_NE(7.0, nz.val) + << "Boundary-column entry (sentinel value 7.0) survived into the " + "submatrix at (" + << nz.row << ", " << nz.col << ")"; + } + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_excludes_boundary_rows +// +// The submatrix must have no rows for boundary cells. Verified by checking +// that MatGetLocalSize equals n_evolving (not n_evolving + n_boundary). +// Also checked by confirming no row index outside [sub_rstart, sub_rend) +// appears in the nonzeros. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_excludes_boundary_rows) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh, /*xin=*/true, /*xout=*/false, + /*yup=*/false, /*ydown=*/true); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + // boundary-row → evolving-column entries + for (int bnd : boundary_stored) { + triples.emplace_back(bnd, 0, 5.0); + } + + const auto op = buildOperatorFromTriples(mapping, triples); + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt local_rows, local_cols; + MatGetLocalSize(sub, &local_rows, &local_cols); + EXPECT_EQ(n_evolving, local_rows); + + PetscInt sub_rstart, sub_rend; + MatGetOwnershipRange(sub, &sub_rstart, &sub_rend); + const auto nzs = matNonzeros(sub); + for (const auto& nz : nzs) { + EXPECT_GE(nz.row, sub_rstart); + EXPECT_LT(nz.row, sub_rend); + } + + MatDestroy(&sub); +} + +// ============================================================================ +// zero_operator_gives_zero_submatrix +// +// A zero operator (no nonzeros) must produce an assembled zero submatrix of +// the correct size. No crash, no uninitialized memory. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, zero_operator_gives_zero_submatrix) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh); + + // Empty triples: zero operator + const auto op = buildOperatorFromTriples(mapping, {}); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + EXPECT_EQ(n_evolving, sub_rows); + EXPECT_EQ(n_evolving, sub_cols); + + MatInfo info; + MatGetInfo(sub, MAT_GLOBAL_SUM, &info); + EXPECT_EQ(0, static_cast(info.nz_used)); + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_of_identity_is_identity +// +// An identity operator over the evolving cells only (no boundary entries) +// should produce an identity submatrix: diagonal entries of 1.0, nothing else. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_of_identity_is_identity) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 1.0); + } + const auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + const auto nzs = matNonzeros(sub); + + PetscInt sub_rstart, sub_rend; + MatGetOwnershipRange(sub, &sub_rstart, &sub_rend); + + ASSERT_EQ(n_evolving, static_cast(nzs.size())); + for (const auto& nz : nzs) { + EXPECT_EQ(nz.row, nz.col) << "Off-diagonal entry in identity submatrix"; + EXPECT_DOUBLE_EQ(1.0, nz.val); + } + + MatDestroy(&sub); +} + +// ============================================================================ +// submatrix_is_independent_of_original_operator +// +// Modifying the original operator's underlying matrix after extraction must +// not affect the already-returned submatrix. MatCreateSubMatrix with +// MAT_INITIAL_MATRIX produces an independent copy. +// ============================================================================ +TEST_F(PetscExtractEvolvingTest, submatrix_is_independent_of_original_operator) { + auto [mapping, n_evolving, evolving_petsc, boundary_stored] = + buildMappingWithBoundaries(mesh); + + std::vector> triples; + for (int i = 0; i < n_evolving; ++i) { + triples.emplace_back(i, i, 2.0); + } + auto op = buildOperatorFromTriples(mapping, triples); + + Mat sub = mapping->extractEvolvingSubmatrix(op); + + // Scale the original operator's matrix in place + MatScale(op.raw(), 0.0); + + // The submatrix must still hold the original values + const auto nzs = matNonzeros(sub); + for (const auto& nz : nzs) { + EXPECT_DOUBLE_EQ(2.0, nz.val) + << "Submatrix was affected by modification of the source operator"; + } + + MatDestroy(&sub); +} + +#endif // BOUT_HAS_PETSC diff --git a/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx b/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx new file mode 100644 index 0000000000..9dbaed584e --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators_make_evolving_is.cxx @@ -0,0 +1,484 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#include "bout/output_bout_types.hxx" +#include "bout/petsc_operators.hxx" +#include "bout/region.hxx" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +#include +#include +#include + +#include + +using bout::globals::mesh; + +// ============================================================================ +// Helper functions shared across unit tests +// ============================================================================ + +namespace { + +/// Build a PetscCellMapping whose evolving region has exactly the cells listed +/// in cell_number with non-negative values, plus whatever boundary cells are in +/// forward_cell_number / backward_cell_number. total_cells must equal the +/// number of distinct non-negative entries across all three fields. +std::shared_ptr +makeMappingFromFields(const Field3D& cell_number, const Field3D& forward_cell_number, + const Field3D& backward_cell_number, int total_cells) { + return std::make_shared(cell_number, forward_cell_number, + backward_cell_number, total_cells); +} + +/// Collect the global PETSc indices that mapOwnedInteriorCells visits, in order. +std::vector collectOwnedInteriorIndices(const PetscCellMapping& mapping) { + std::vector indices; + mapping.mapOwnedInteriorCells( + [&](PetscInt row, const Ind3D& /*i*/, int /*stored*/) { indices.push_back(row); }); + return indices; +} + +/// Unwrap an IS into a sorted std::vector. The IS is not destroyed. +std::vector isToSortedVector(IS is) { + PetscInt n; + ISGetLocalSize(is, &n); + const PetscInt* idx; + ISGetIndices(is, &idx); + std::vector v(idx, idx + n); + ISRestoreIndices(is, &idx); + std::sort(v.begin(), v.end()); + return v; +} +} // Namespace + +// ============================================================================ +// Fixture +// ============================================================================ + +// FakeMeshFixture provides a mesh with: +// LocalNx = 5 (xstart=1, xend=3, so mxg=1, 3 interior x-points) +// LocalNy = 4 (ystart=1, yend=2, so mgy=1, 2 interior y-points) +// LocalNz = 2 +// Interior cells: x in [1,3], y in [1,2], z in [0,1] => 3*2*2 = 12 evolving cells. +using PetscEvolvingISTest = FakeMeshFixture; + +// ============================================================================ +// Helper: build the standard mapping used by most tests in this suite. +// +// Stored cell numbers are assigned to interior cells only (indices 0..11). +// No forward or backward boundary cells. +// ============================================================================ +namespace { +std::shared_ptr buildStandardMapping(Mesh* mesh) { + Field3D cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + return makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, + next); +} +} // Namespace + +// ============================================================================ +// IS_size_equals_n_evolving +// +// The IS produced by makeEvolvingIS() must contain exactly one index per +// evolving interior cell. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_size_equals_n_evolving) { + auto mapping = buildStandardMapping(mesh); + + // FakeMesh: 3 x-interior * 2 y-interior * 2 z = 12 evolving cells + const PetscInt expected_n_evolving = + (mesh->xend - mesh->xstart + 1) * (mesh->yend - mesh->ystart + 1) * mesh->LocalNz; + + IS is = mapping->makeEvolvingIS(); + + PetscInt local_size; + ISGetLocalSize(is, &local_size); + + EXPECT_EQ(expected_n_evolving, local_size); + + ISDestroy(&is); +} + +// ============================================================================ +// IS_indices_match_mapOwnedInteriorCells_order +// +// The IS indices must appear in exactly the same order as the global PETSc +// row indices produced by mapOwnedInteriorCells. Both traversals must visit +// the same cells in the same sequence. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_indices_match_mapOwnedInteriorCells_order) { + auto mapping = buildStandardMapping(mesh); + + const std::vector from_map = collectOwnedInteriorIndices(*mapping); + + IS is = mapping->makeEvolvingIS(); + PetscInt n; + ISGetLocalSize(is, &n); + const PetscInt* idx; + ISGetIndices(is, &idx); + const std::vector from_is(idx, idx + n); + ISRestoreIndices(is, &idx); + ISDestroy(&is); + + ASSERT_EQ(from_map.size(), from_is.size()); + for (std::size_t i = 0; i < from_map.size(); ++i) { + EXPECT_EQ(from_map[i], from_is[i]) << "Index mismatch at position " << i; + } +} + +// ============================================================================ +// IS_indices_are_in_global_petsc_range +// +// On a single MPI rank the global PETSc row range is [rowStart(), rowEnd()). +// Every index in the IS must lie in that range. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_indices_are_in_global_petsc_range) { + auto mapping = buildStandardMapping(mesh); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + const PetscInt row_start = mapping->rowStart(); + const PetscInt row_end = mapping->rowEnd(); + + for (PetscInt idx : v) { + EXPECT_GE(idx, row_start) << "Index " << idx << " is below rowStart()"; + EXPECT_LT(idx, row_end) << "Index " << idx << " is at or above rowEnd()"; + } +} + +// ============================================================================ +// IS_excludes_xin_boundary_cells +// +// Cells in the inner-X boundary region have stored indices assigned by the +// mapping but must NOT appear in the evolving IS. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_xin_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + // Assign evolving cells + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // Assign one inner-X boundary cell a stored index + const int xin_stored = next++; + cell_number(0, mesh->ystart, 0) = xin_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt xin_petsc = mapping->storedToPetsc(xin_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), xin_petsc)) + << "Inner-X boundary cell (PETSc index " << xin_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_excludes_xout_boundary_cells +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_xout_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // Outer-X boundary: x == LocalNx - 1 + const int xout_stored = next++; + cell_number(mesh->LocalNx - 1, mesh->ystart, 0) = xout_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt xout_petsc = mapping->storedToPetsc(xout_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), xout_petsc)) + << "Outer-X boundary cell (PETSc index " << xout_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_excludes_forward_boundary_cells +// +// Virtual forward-parallel boundary cells (yup) live in forward_cell_number, +// not cell_number. They must not appear in the evolving IS. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_forward_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // One virtual forward boundary cell + const int fwd_stored = next++; + forward_cell_number(mesh->xstart, mesh->yend, 0) = fwd_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt fwd_petsc = mapping->storedToPetsc(fwd_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), fwd_petsc)) + << "Forward boundary cell (PETSc index " << fwd_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_excludes_backward_boundary_cells +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_excludes_backward_boundary_cells) { + Field3D cell_number{-1.0, mesh}; + const Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + const int bwd_stored = next++; + backward_cell_number(mesh->xstart, mesh->ystart, 0) = bwd_stored; + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + const PetscInt bwd_petsc = mapping->storedToPetsc(bwd_stored); + + IS is = mapping->makeEvolvingIS(); + const auto v = isToSortedVector(is); + ISDestroy(&is); + + EXPECT_EQ(v.end(), std::find(v.begin(), v.end(), bwd_petsc)) + << "Backward boundary cell (PETSc index " << bwd_petsc + << ") should not appear in the evolving IS"; +} + +// ============================================================================ +// IS_is_empty_when_no_evolving_cells +// +// If cell_number is entirely -1, there are no evolving cells and the IS must +// have size zero. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_is_empty_when_no_evolving_cells) { + const Field3D cell_number{-1.0, mesh}; + const Field3D forward_cell_number{-1.0, mesh}; + const Field3D backward_cell_number{-1.0, mesh}; + + // total_cells = 0: no cells at all + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, 0); + + IS is = mapping->makeEvolvingIS(); + PetscInt local_size; + ISGetLocalSize(is, &local_size); + ISDestroy(&is); + + EXPECT_EQ(0, local_size); +} + +// ============================================================================ +// IS_covers_all_evolving_stored_indices +// +// The set of stored indices reachable from the IS (via storedToPetsc inverse) +// must exactly equal the set of stored indices of interior cells as visited by +// mapOwnedInteriorCells. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_covers_all_evolving_stored_indices) { + auto mapping = buildStandardMapping(mesh); + + // Collect stored indices from mapOwnedInteriorCells + std::vector stored_from_map; + mapping->mapOwnedInteriorCells([&](PetscInt /*row*/, const Ind3D& /*i*/, int stored) { + stored_from_map.push_back(stored); + }); + std::sort(stored_from_map.begin(), stored_from_map.end()); + + // Collect PETSc indices from IS, then map back to stored via localStoredIndices + IS is = mapping->makeEvolvingIS(); + PetscInt n; + ISGetLocalSize(is, &n); + const PetscInt* idx; + ISGetIndices(is, &idx); + + // localStoredIndices[i] is the stored index for PETSc row rowStart()+i + const auto& local_stored = mapping->localStoredIndices(); + const PetscInt rstart = mapping->rowStart(); + std::vector stored_from_is; + for (PetscInt k = 0; k < n; ++k) { + const PetscInt local_row = idx[k] - rstart; + ASSERT_GE(local_row, 0); + ASSERT_LT(local_row, static_cast(local_stored.size())); + stored_from_is.push_back(local_stored[local_row]); + } + ISRestoreIndices(is, &idx); + ISDestroy(&is); + + std::sort(stored_from_is.begin(), stored_from_is.end()); + + EXPECT_EQ(stored_from_map, stored_from_is); +} + +// ============================================================================ +// IS_indices_are_unique +// +// No PETSc index should appear more than once in the IS. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_indices_are_unique) { + auto mapping = buildStandardMapping(mesh); + + IS is = mapping->makeEvolvingIS(); + const auto sorted = isToSortedVector(is); + ISDestroy(&is); + + for (std::size_t i = 1; i < sorted.size(); ++i) { + EXPECT_LT(sorted[i - 1], sorted[i]) + << "Duplicate index " << sorted[i] << " found at positions " << i - 1 << " and " + << i; + } +} + +// ============================================================================ +// IS_with_mixed_boundary_cells_has_correct_size +// +// When forward, backward, and X-boundary cells are all present alongside +// evolving cells, the IS size equals only the evolving cell count. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_with_mixed_boundary_cells_has_correct_size) { + Field3D cell_number{-1.0, mesh}; + Field3D forward_cell_number{-1.0, mesh}; + Field3D backward_cell_number{-1.0, mesh}; + + int next = 0; + // Evolving interior cells + const int n_evolving = + (mesh->xend - mesh->xstart + 1) * (mesh->yend - mesh->ystart + 1) * mesh->LocalNz; + for (int x = mesh->xstart; x <= mesh->xend; ++x) { + for (int y = mesh->ystart; y <= mesh->yend; ++y) { + for (int z = 0; z < mesh->LocalNz; ++z) { + cell_number(x, y, z) = next++; + } + } + } + // One of each boundary type + cell_number(0, mesh->ystart, 0) = next++; // xin + cell_number(mesh->LocalNx - 1, mesh->ystart, 0) = next++; // xout + forward_cell_number(mesh->xstart, mesh->yend, 0) = next++; // yup + backward_cell_number(mesh->xstart, mesh->ystart, 0) = next++; // ydown + + auto mapping = + makeMappingFromFields(cell_number, forward_cell_number, backward_cell_number, next); + + IS is = mapping->makeEvolvingIS(); + PetscInt local_size; + ISGetLocalSize(is, &local_size); + ISDestroy(&is); + + EXPECT_EQ(n_evolving, local_size); +} + +// ============================================================================ +// IS_can_be_used_to_create_submatrix +// +// Regression guard: the IS must be acceptable to MatCreateSubMatrix without +// error. Uses a simple identity-pattern matrix over the full cell space. +// ============================================================================ +TEST_F(PetscEvolvingISTest, IS_can_be_used_to_create_submatrix) { + auto mapping = buildStandardMapping(mesh); + + // Build a diagonal matrix over the full cell space + const PetscInt n = mapping->globalSize(); + const PetscInt nlocal = mapping->localSize(); + Mat A; + MatCreate(BoutComm::get(), &A); + MatSetSizes(A, nlocal, nlocal, n, n); + MatSetType(A, MATMPIAIJ); + MatMPIAIJSetPreallocation(A, 1, nullptr, 0, nullptr); + const PetscInt rstart = mapping->rowStart(); + for (PetscInt i = 0; i < nlocal; ++i) { + const PetscInt global = rstart + i; + const PetscScalar val = 1.0; + MatSetValues(A, 1, &global, 1, &global, &val, INSERT_VALUES); + } + MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); + + IS is = mapping->makeEvolvingIS(); + + Mat sub = nullptr; + // This must not throw or return a PETSc error + const PetscErrorCode ierr = MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &sub); + + EXPECT_EQ(PETSC_SUCCESS, ierr); + + // The submatrix size should equal n_evolving × n_evolving + if (sub != nullptr) { + PetscInt sub_rows, sub_cols; + MatGetSize(sub, &sub_rows, &sub_cols); + const PetscInt n_evolving = + (mesh->xend - mesh->xstart + 1) * (mesh->yend - mesh->ystart + 1) * mesh->LocalNz; + EXPECT_EQ(n_evolving, sub_rows); + EXPECT_EQ(n_evolving, sub_cols); + MatDestroy(&sub); + } + + ISDestroy(&is); + MatDestroy(&A); +} + +#endif // BOUT_HAS_PETSC diff --git a/tests/unit/solver/test_solver.cxx b/tests/unit/solver/test_solver.cxx index 525d52767f..9f8367b684 100644 --- a/tests/unit/solver/test_solver.cxx +++ b/tests/unit/solver/test_solver.cxx @@ -381,6 +381,39 @@ TEST_F(SolverTest, AddVector3D) { EXPECT_EQ(solver.listVector3DNames(), expected_names); } +TEST_F(SolverTest, GetVarRef) { + Options options; + FakeSolver solver{&options}; + + Field2D field2d{}; + Field3D field3d{}; + Vector3D vector{}; + + solver.add(field2d, "phi"); + solver.add(field3d, "n"); + solver.add(vector, "v"); + + const auto phi = solver.getVarRef("phi"); + const auto n = solver.getVarRef("n"); + const auto vx = solver.getVarRef("v_x"); + const auto vy = solver.getVarRef("v_y"); + const auto vz = solver.getVarRef("v_z"); + const auto missing = solver.getVarRef("missing"); + + EXPECT_TRUE(phi.isConcrete()); + EXPECT_EQ(phi.index(), 0); + EXPECT_TRUE(n.isConcrete()); + EXPECT_EQ(n.index(), 1); + EXPECT_TRUE(vx.isConcrete()); + EXPECT_EQ(vx.index(), 2); + EXPECT_TRUE(vy.isConcrete()); + EXPECT_EQ(vy.index(), 3); + EXPECT_TRUE(vz.isConcrete()); + EXPECT_EQ(vz.index(), 4); + EXPECT_TRUE(missing.isInvalid()); + EXPECT_TRUE(Solver::VarRef::All().isAll()); +} + TEST_F(SolverTest, ConstraintField2D) { Options options; FakeSolver solver{&options};