diff --git a/CMakeLists.txt b/CMakeLists.txt index 5f054c1c1f..0b14e5c3d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -177,6 +177,7 @@ set(BOUT_SOURCES ./include/bout/parallel_boundary_region.hxx ./include/bout/paralleltransform.hxx ./include/bout/petsc_interface.hxx + ./include/bout/petsc_operators.hxx ./include/bout/petsclib.hxx ./include/bout/physicsmodel.hxx ./include/bout/rajalib.hxx @@ -308,6 +309,7 @@ set(BOUT_SOURCES ./src/mesh/parallel/shiftedmetricinterp.hxx ./src/mesh/parallel_boundary_op.cxx ./src/mesh/parallel_boundary_region.cxx + ./src/mesh/petsc_operators.cxx ./src/mesh/surfaceiter.cxx ./src/physics/gyro_average.cxx ./src/physics/physicsmodel.cxx diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index f2444e6bf6..5eabaf045d 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -18,6 +18,7 @@ add_subdirectory(dalf3) add_subdirectory(eigen-box) add_subdirectory(elm-pb) add_subdirectory(elm-pb-outerloop) +add_subdirectory(fci-operators) add_subdirectory(fci-wave) add_subdirectory(finite-volume/diffusion) add_subdirectory(finite-volume/fluid) diff --git a/examples/fci-operators/CMakeLists.txt b/examples/fci-operators/CMakeLists.txt new file mode 100644 index 0000000000..3a19ea9dcc --- /dev/null +++ b/examples/fci-operators/CMakeLists.txt @@ -0,0 +1,14 @@ +cmake_minimum_required(VERSION 3.13) + +project(fci-operators LANGUAGES CXX) + +if(NOT TARGET bout++::bout++) + find_package(bout++ REQUIRED) +endif() + +bout_add_example( + fci-operators + SOURCES fci_operators_example.cxx + DATA_DIRS data + EXTRA_FILES rotating-ellipse.py makeplots.py +) diff --git a/examples/fci-operators/data/BOUT.inp b/examples/fci-operators/data/BOUT.inp new file mode 100644 index 0000000000..182d49cb96 --- /dev/null +++ b/examples/fci-operators/data/BOUT.inp @@ -0,0 +1,5 @@ +[mesh] +file = "rotating-ellipse-20x8x20.fci.nc" + +[mesh:paralleltransform] +type = fci diff --git a/examples/fci-operators/fci_operators_example.cxx b/examples/fci-operators/fci_operators_example.cxx new file mode 100644 index 0000000000..3366622ff0 --- /dev/null +++ b/examples/fci-operators/fci_operators_example.cxx @@ -0,0 +1,61 @@ +#include "bout/bout.hxx" +#include "bout/difops.hxx" +#include "bout/field.hxx" +#include "bout/field3d.hxx" +#include "bout/field_factory.hxx" +#include "bout/globals.hxx" +#include "bout/options.hxx" +#include "bout/petsc_operators.hxx" + +int main(int argc, char** argv) { + BoutInitialise(argc, argv); + + const PetscOperators ops(bout::globals::mesh); + + // Parallel operators + const auto parallel = ops.getParallel(); + + // Test field + Field3D f = FieldFactory::get()->create3D("sin(y) + cos(z)"); + bout::globals::mesh->communicate(f); + f.applyParallelBoundary("parallel_dirichlet_o2"); + + Field3D f_neumann = FieldFactory::get()->create3D("sin(y) + cos(z)"); + bout::globals::mesh->communicate(f_neumann); + f_neumann.applyParallelBoundary("parallel_neumann_o1"); + + auto forward_op = parallel.Restrict_minus * parallel.Forward; + + Options dump; + dump["f"] = f; + dump["dV"] = parallel.dV; + + // Test1 : (parallel.Restrict_minus * parallel.Inject_plus) is the identity in the interior (not X boundaries) + + dump["grad_par_op"] = parallel.Grad_par(f); + Field3D forward{0.0}; + BOUT_FOR(i, forward.getRegion("RGN_NOBNDRY")) { forward[i] = f.yup()[i.yp()]; } + + dump["forward"] = forward; + + dump["forward_op"] = forward_op(f); + + dump["grad_par_yud"] = Grad_par(f); + + // Test2 : Sum of dV * Div_par over domain is zero + + dump["div_par_op"] = parallel.Div_par(f); + + auto* coords = bout::globals::mesh->getCoordinates(); + bout::globals::mesh->communicate(coords->Bxy); + coords->Bxy.applyParallelBoundary("parallel_neumann_o1"); + dump["div_par_yud"] = Div_par(f); + + dump["div_par_grad_par_op"] = parallel.Div_par_Grad_par(f_neumann); + dump["div_par_grad_par_yud"] = Div_par_K_Grad_par(1.0, f_neumann); + + bout::writeDefaultOutputFile(dump); + + BoutFinalise(); + return 0; +} diff --git a/examples/fci-operators/makeplots.py b/examples/fci-operators/makeplots.py new file mode 100644 index 0000000000..53c3d38f05 --- /dev/null +++ b/examples/fci-operators/makeplots.py @@ -0,0 +1,90 @@ +from boutdata import collect +import matplotlib.pyplot as plt +import numpy as np + +path = "data" + +forward_op = collect("forward_op", path=path) +forward = collect("forward", path=path) + +grad_par_op = collect("grad_par_op", path=path) +grad_par_yud = collect("grad_par_yud", path=path) + +div_par_op = collect("div_par_op", path=path) +div_par_yud = collect("div_par_yud", path=path) + +dV = collect("dV", path=path) # Cell volume + +# Integrate divergence +op_div_int = np.sum((dV * div_par_op)[2:-2, :, :]) +yud_div_int = np.sum((dV * div_par_yud)[2:-2, :, :]) + +print(f"Div_par volume integral Op: {op_div_int} Yup/down: {yud_div_int}") + +div_par_grad_par_op = collect("div_par_grad_par_op", path=path) +div_par_grad_par_yud = collect("div_par_grad_par_yud", path=path) + +# Integrate divergence +op_div2_int = np.sum((dV * div_par_grad_par_op)[2:-2, :, :]) +yud_div2_int = np.sum((dV * div_par_grad_par_yud)[2:-2, :, :]) + +print(f"Div_par_Grad_par volume integral Op: {op_div2_int} Yup/down: {yud_div2_int}") + + +def plot_comparison(a, alabel, b, blabel): + fig, axs = plt.subplots(1, 3, figsize=(10, 3)) + vmax = max([np.amax(a), np.amax(b)]) + vmin = min([np.amin(a), np.amin(b)]) + + cs = axs[0].contourf(a, vmin=vmin, vmax=vmax) + fig.colorbar(cs, ax=axs[0]) + axs[0].set_title(alabel) + + axs[1].contourf(b, vmin=vmin, vmax=vmax) + axs[1].set_title(blabel) + + cs = axs[2].contourf(a - b) + axs[2].set_title("Difference") + fig.colorbar(cs, ax=axs[2]) + + return fig + + +fig = plot_comparison( + forward_op[2:-2, 2, :], "Operator", forward[2:-2, 2, :], "Yup/down" +) +fig.suptitle("Forward") +fig.tight_layout() + +fig = plot_comparison( + grad_par_op[2:-2, 2, :], "Operator", grad_par_yud[2:-2, 2, :], "Yup/down" +) +fig.suptitle("Parallel Gradient") +fig.tight_layout() +fig.savefig("gradient.pdf") +fig.savefig("gradient.png") + +fig = plot_comparison( + div_par_op[2:-2, 2, :], + f"Operator [integral {op_div_int:.2e}]", + div_par_yud[2:-2, 2, :], + f"Yup/down [integral {yud_div_int:.2e}]", +) +fig.suptitle("Parallel Divergence") +fig.tight_layout() +fig.savefig("divergence.pdf") +fig.savefig("divergence.png") + +fig = plot_comparison( + div_par_grad_par_op[2:-2, 2, :], + f"Operator [integral {op_div2_int:.2e}]", + div_par_grad_par_yud[2:-2, 2, :], + f"Yup/down [integral {yud_div2_int:.2e}]", +) +fig.suptitle("Parallel diffusion") +fig.tight_layout() +fig.savefig("diffusion.pdf") +fig.savefig("diffusion.png") + + +plt.show() diff --git a/examples/fci-operators/rotating-ellipse.py b/examples/fci-operators/rotating-ellipse.py new file mode 100755 index 0000000000..3aa485c093 --- /dev/null +++ b/examples/fci-operators/rotating-ellipse.py @@ -0,0 +1,102 @@ +#!/usr/bin/env python + +# Create grids for a rotating ellipse stellarator, based on curvilinear grids + +import numpy as np + +import zoidberg + +############################################################################# +# Define the magnetic field + +# Length in y after which the coils return to their starting (R,Z) locations +yperiod = 2.0 * np.pi / 5 + +xcentre = 2.0 + +magnetic_field = zoidberg.field.RotatingEllipse( + xcentre=xcentre, # Major radius of axis [m] + radius=0.8, # Minor radius of the coils [m] + I_coil=0.4, # Coil current + yperiod=yperiod, + Btor=1.0, # Toroidal magnetic field +) + +############################################################################# +# Create the inner flux surface, starting at a point at phi=0 +# To do this we need to define the y locations of the poloidal points +# where we will construct grids + +start_r = xcentre + 0.4 +start_z = 0.0 + +nslices = 8 # Number of poloidal slices +ycoords = np.linspace(0, yperiod, nslices) +npoints = 20 # Points per poloidal slice + +rzcoord, ycoords = zoidberg.fieldtracer.trace_poincare( + magnetic_field, start_r, start_z, yperiod, y_slices=ycoords, revs=npoints, nover=1 +) + +inner_lines = [] +for i in range(nslices): + r = rzcoord[:, i, 0, 0] + z = rzcoord[:, i, 0, 1] + line = zoidberg.rzline.line_from_points(r, z) + # Re-map the points so they're approximately uniform in distance along the surface + # Note that this results in some motion of the line + line = line.equallySpaced() + inner_lines.append(line) + +# Now have a list of y coordinates (ycoords) and inner lines (inner_lines) + +############################################################################# +# Generate a fixed circle for the outer boundary + +outer_line = zoidberg.rzline.circle(R0=xcentre, r=0.6) + +############################################################################# +# Now have inner and outer boundaries for each poloidal grid +# Generate a grid on each poloidal slice using the elliptic grid generator + +nx = 20 +nz = 20 + +pol_grids = [ + zoidberg.poloidal_grid.grid_elliptic(inner_line, outer_line, nx, nz) + for inner_line in inner_lines +] + +############################################################################# +# Create a grid, then calculate forward and backward maps + +grid = zoidberg.grid.Grid(pol_grids, ycoords, yperiod, yperiodic=True) + +samples_per_dim = 5 # Sub-sampling in each cell + +maps = zoidberg.make_maps(grid, magnetic_field, samples_per_dim=samples_per_dim) + +############################################################################# +# Interpolation weights + +weight_vars = zoidberg.weights.calculate_parallel_map_operators(maps) +maps.update(weight_vars) +# Remove sub-sampled maps from output +for var in [ + "forward_xt_primes", + "forward_zt_primes", + "backward_xt_primes", + "backward_zt_primes", + "subcell_weights", +]: + maps.pop(var, None) + +############################################################################# +# Write grid file + +filename = f"rotating-ellipse-{nx}x{nslices}x{nz}-s{samples_per_dim}.fci.nc" + +print(f"Writing to grid file '{filename}'") +zoidberg.write_maps( + grid, magnetic_field, maps, gridfile=filename, new_names=False, metric2d=False +) diff --git a/include/bout/array.hxx b/include/bout/array.hxx index 82677fd38a..d8ce207dac 100644 --- a/include/bout/array.hxx +++ b/include/bout/array.hxx @@ -27,6 +27,7 @@ #define BOUT_ARRAY_H #include +#include #include #include #include @@ -187,9 +188,22 @@ public: Array() noexcept : ptr(nullptr) {} /*! - * Create an array of given length + * Create an array of given length. */ - Array(size_type len) { ptr = get(len); } + Array(size_type len) : ptr(get(len)) {} + + /*! + * Create an array with initializer list. + * This is explicit to avoid confusion with (size_type) constructor. + */ + static Array fromValues(std::initializer_list init) { + if (init.size() == 0) { + return Array(); + } + Array array(init.size()); + std::copy(init.begin(), init.end(), array.begin()); + return array; + } /*! * Destructor. Releases the underlying dataBlock @@ -199,7 +213,7 @@ public: /*! * Copy constructor */ - Array(const Array& other) noexcept { ptr = other.ptr; } + Array(const Array& other) noexcept : ptr(other.ptr) {} /*! * Assignment operator diff --git a/include/bout/field3d.hxx b/include/bout/field3d.hxx index 7ba6ae9575..5e97e3ebc9 100644 --- a/include/bout/field3d.hxx +++ b/include/bout/field3d.hxx @@ -1,7 +1,7 @@ /************************************************************************** - * 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++. * @@ -349,6 +349,14 @@ public: return std::end(getRegion("RGN_ALL")); }; + /// Convert (x, y, z) to an Ind3D + /// Requires mesh size. + Ind3D indexAt(int x, int y, int z) const { + const int ny = this->getNy(); + const int nz = this->getNz(); + return Ind3D{(((x * ny) + y) * nz) + z, ny, nz}; + } + BoutReal& operator[](const Ind3D& d) { return data[d.ind]; } const BoutReal& operator[](const Ind3D& d) const { return data[d.ind]; } diff --git a/include/bout/griddata.hxx b/include/bout/griddata.hxx index 29a32e5779..e9d9214846 100644 --- a/include/bout/griddata.hxx +++ b/include/bout/griddata.hxx @@ -2,9 +2,9 @@ * Functions to read input variables on a grid/mesh * ************************************************************************** - * 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++. * @@ -29,12 +29,16 @@ class GridDataSource; #define BOUT_GRIDDATA_H #include "mesh.hxx" +#include "bout/array.hxx" #include "bout/bout_types.hxx" #include "bout/options.hxx" #include #include +#include +#include + /// Interface class to serve grid data /*! * Provides a generic interface for sources of @@ -71,6 +75,15 @@ public: static constexpr Direction Y = Direction::Y; static constexpr Direction Z = Direction::Z; + virtual bool get([[maybe_unused]] Array& var, + [[maybe_unused]] const std::string& name) { + return false; + } + virtual bool get([[maybe_unused]] Array& var, + [[maybe_unused]] const std::string& name) { + return false; + } + virtual bool get(Mesh* m, std::vector& var, const std::string& name, int len, int offset = 0, Direction dir = GridDataSource::X) = 0; virtual bool get(Mesh* m, std::vector& var, const std::string& name, int len, @@ -115,6 +128,9 @@ public: return getField(m, var, name, def, location); } + bool get(Array& var, const std::string& name) override; + bool get(Array& var, const std::string& name) override; + bool get(Mesh* m, std::vector& var, const std::string& name, int len, int offset = 0, GridDataSource::Direction dir = GridDataSource::X) override; bool get(Mesh* m, std::vector& var, const std::string& name, int len, diff --git a/include/bout/mesh.hxx b/include/bout/mesh.hxx index 25e66fe5d2..c1750fa414 100644 --- a/include/bout/mesh.hxx +++ b/include/bout/mesh.hxx @@ -43,6 +43,7 @@ class Mesh; #ifndef BOUT_MESH_H #define BOUT_MESH_H +#include "bout/array.hxx" #include "bout/bout_enum_class.hxx" #include "bout/bout_types.hxx" #include "bout/coordinates.hxx" // Coordinates class @@ -169,6 +170,9 @@ public: /// @returns zero if successful, non-zero on failure int get(bool& bval, const std::string& name, bool def = false); + int get(Array& var, const std::string& name); + int get(Array& var, const std::string& name); + /// Get a Field2D from the input source /// including communicating guard cells /// diff --git a/include/bout/parallel_boundary_region.hxx b/include/bout/parallel_boundary_region.hxx index 4d5278d00f..0833c8d3bc 100644 --- a/include/bout/parallel_boundary_region.hxx +++ b/include/bout/parallel_boundary_region.hxx @@ -108,6 +108,8 @@ public: void next() final { ++bndry_position; } bool isDone() final { return (bndry_position == end(bndry_points)); } + std::size_t size() { return this->bndry_points.size(); } + // getter Ind3D ind() const { return bndry_position->index; } BoutReal s_x() const { return bndry_position->intersection.s_x; } diff --git a/include/bout/petsc_interface.hxx b/include/bout/petsc_interface.hxx index 2ce71d0549..cfe27597ff 100644 --- a/include/bout/petsc_interface.hxx +++ b/include/bout/petsc_interface.hxx @@ -62,6 +62,36 @@ #include +namespace bout { +namespace petsc { +struct VecDeleter { + void operator()(Vec* vec_ptr) const { + VecDestroy(vec_ptr); + delete vec_ptr; + } +}; + +/// Unique pointer wrapper around PETSc Vec +/// Assumes that the Vec has been allocated on the heap. +/// +/// e.g. UniqueVec(new Vec{}); +using UniqueVec = std::unique_ptr; + +struct MatDeleter { + void operator()(Mat* mat_ptr) const { + MatDestroy(mat_ptr); + delete mat_ptr; + } +}; + +/// Unique pointer wrapper around PETSc Mat +/// Assumes that the Mat has been allocated on the heap. +/// +/// e.g. UniqueMat(new Mat{}); +using UniqueMat = std::unique_ptr; +} // namespace petsc +} // namespace bout + /*! * A class which wraps PETSc vector objects, allowing them to be * indexed using the BOUT++ scheme. Note that boundaries are only @@ -70,7 +100,7 @@ template class PetscVector; template -void swap(PetscVector& first, PetscVector& second); +void swap(PetscVector& first, PetscVector& second) noexcept; template inline MPI_Comm getComm([[maybe_unused]] const T& field) { @@ -88,13 +118,6 @@ public: static_assert(bout::utils::is_Field_v, "PetscVector only works with Fields"); using ind_type = typename T::ind_type; - struct VectorDeleter { - void operator()(Vec* vec) const { - VecDestroy(vec); - delete vec; - } - }; - PetscVector() : vector(new Vec{}) {} ~PetscVector() = default; @@ -160,7 +183,7 @@ public: return *this; } - friend void swap(PetscVector& first, PetscVector& second); + friend void swap(PetscVector& first, PetscVector& second) noexcept; BoutReal& operator()(const ind_type& index) { #if CHECKLEVEL >= 1 @@ -239,23 +262,23 @@ public: const Vec* get() const { return vector.get(); } private: - PetscLib lib{}; - std::unique_ptr vector = nullptr; + PetscLib lib; + bout::petsc::UniqueVec vector = nullptr; IndexerPtr indexConverter{}; CELL_LOC location = CELL_LOC::deflt; bool initialised = false; - Array vector_values{}; + Array vector_values; }; /*! - * A class which wraps PETSc vector objects, allowing them to be + * A class which wraps PETSc matrix objects, allowing them to be * indexed using the BOUT++ scheme. It provides the option of setting * a y-offset that interpolates onto field lines. */ template class PetscMatrix; template -void swap(PetscMatrix& first, PetscMatrix& second); +void swap(PetscMatrix& first, PetscMatrix& second) noexcept; template class PetscMatrix { @@ -263,13 +286,6 @@ public: static_assert(bout::utils::is_Field_v, "PetscMatrix only works with Fields"); using ind_type = typename T::ind_type; - struct MatrixDeleter { - void operator()(Mat* mat) const { - MatDestroy(mat); - delete mat; - } - }; - /// Default constructor does nothing PetscMatrix() : matrix(new Mat()) {} ~PetscMatrix() = default; @@ -289,7 +305,7 @@ public: } // Construct a matrix capable of operating on the specified field, - // preallocating memory if requeted and possible. + // preallocating memory if requested and possible. PetscMatrix(IndexerPtr indConverter, bool preallocate = true) : matrix(new Mat()), indexConverter(indConverter), pt(&indConverter->getMesh()->getCoordinates()->getParallelTransform()) { @@ -327,7 +343,7 @@ public: rhs.initialised = false; return *this; } - friend void swap(PetscMatrix& first, PetscMatrix& second); + friend void swap(PetscMatrix& first, PetscMatrix& second) noexcept; /*! * A class which is used to assign to a particular element of a PETSc @@ -540,7 +556,7 @@ private: * vice versa. */ template -void swap(PetscVector& first, PetscVector& second) { +void swap(PetscVector& first, PetscVector& second) noexcept { std::swap(first.vector, second.vector); std::swap(first.indexConverter, second.indexConverter); std::swap(first.location, second.location); @@ -552,7 +568,7 @@ void swap(PetscVector& first, PetscVector& second) { * vice versa. */ template -void swap(PetscMatrix& first, PetscMatrix& second) { +void swap(PetscMatrix& first, PetscMatrix& second) noexcept { std::swap(first.matrix, second.matrix); std::swap(first.indexConverter, second.indexConverter); std::swap(first.pt, second.pt); diff --git a/include/bout/petsc_operators.hxx b/include/bout/petsc_operators.hxx new file mode 100644 index 0000000000..0d9c6c6912 --- /dev/null +++ b/include/bout/petsc_operators.hxx @@ -0,0 +1,1024 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "bout/array.hxx" +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/field3d.hxx" +#include "bout/mesh.hxx" +#include "bout/output.hxx" +#include "bout/output_bout_types.hxx" +#include "bout/petsc_interface.hxx" +#include "bout/petsclib.hxx" +#include "bout/region.hxx" +#include "bout/utils.hxx" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +/// @brief Tag type identifying the cell space C. +/// +/// Used as a compile-time template parameter to distinguish vectors and operators +/// that live in the cell space from those in the forward or backward leg spaces. +/// Mismatched space tags produce a compile error rather than a runtime failure. +struct CellSpaceTag {}; + +/// @brief Tag type identifying the forward leg space L+. +/// +/// Operators and vectors tagged with this type map to or from the forward (y+1) +/// leg space. See @ref BackwardLegSpaceTag and @ref CellSpaceTag. +struct ForwardLegSpaceTag {}; + +/// @brief Tag type identifying the backward leg space L-. +/// +/// Operators and vectors tagged with this type map to or from the backward (y-1) +/// leg space. See @ref ForwardLegSpaceTag and @ref CellSpaceTag. +struct BackwardLegSpaceTag {}; + +/// @brief Bidirectional index mapping between mesh-file stored numbering and PETSc +/// row ownership. +/// +/// Every vector space (cell space, forward leg space, backward leg space) is defined +/// by a global stored numbering written to the mesh file by the Python weights module. +/// PETSc distributes rows contiguously across MPI ranks, so the local PETSc row range +/// [rowStart(), rowEnd()) does not in general coincide with the stored indices. +/// +/// This class maintains: +/// - the list of stored indices owned by the local rank, in PETSc row order; +/// - a hash map for O(1) stored → PETSc lookup; +/// - a permutation matrix @c mat_stored_to_petsc (PETSc ordering → stored numbering) +/// and its transpose @c mat_petsc_to_stored (stored numbering → PETSc ordering), +/// used to translate operator column indices during matrix assembly. +/// +/// Subclasses specialise the constructor to populate these structures for a particular +/// space from mesh metadata. +class PetscIndexMapping { +public: + PetscIndexMapping() = default; + virtual ~PetscIndexMapping(); + + PetscIndexMapping(const PetscIndexMapping&) = delete; + PetscIndexMapping& operator=(const PetscIndexMapping&) = delete; + PetscIndexMapping(PetscIndexMapping&&) = delete; + PetscIndexMapping& operator=(PetscIndexMapping&&) = delete; + + /// @brief Number of rows owned locally by this MPI rank. + PetscInt localSize() const { + return static_cast(local_stored_indices.size()); + } + + /// @brief Total number of rows in this space across all MPI ranks. + PetscInt globalSize() const { return global_size; } + + /// @brief Global PETSc index of the first row owned by this MPI rank. + PetscInt rowStart() const { return row_start; } + + /// @brief Global PETSc index one past the last row owned by this MPI rank. + /// + /// The locally owned rows are the half-open interval [rowStart(), rowEnd()). + PetscInt rowEnd() const { return row_end; } + + /// @brief Stored mesh-file indices owned locally, listed in PETSc row order. + /// + /// Entry @c i corresponds to global PETSc row rowStart() + i. + const std::vector& localStoredIndices() const { return local_stored_indices; } + + /// @brief Convert a locally owned stored index to its global PETSc row index. + /// + /// @param stored_index A mesh-file stored index that must be owned by this rank. + /// @returns The corresponding global PETSc row index. + /// @throws BoutException if @p stored_index is not owned locally. + PetscInt storedToPetsc(int stored_index) const; + + /// @brief Return the permutation matrix that maps PETSc row ordering to stored + /// mesh-file numbering. + /// + /// This is the transpose of @c mat_stored_to_petsc and is used in + /// @ref PetscOperator::assembleFromCSR to translate stored column indices into + /// PETSc column indices via a post-multiplication. + Mat getPetscToStored() const { return mat_petsc_to_stored; } + +protected: + /// @brief BOUT++ PETSc library handle; ensures PETSc is initialised. + PetscLib lib; + + /// @brief Construct the permutation matrices and populate internal lookup structures. + /// + /// Assigns contiguous PETSc rows starting from the rank-local offset determined by + /// PETSc's ownership partitioning. Each stored index in @p stored_indices maps to + /// one PETSc row, in the order given. + /// + /// @param nlocal Number of rows owned by this rank. + /// @param nglobal Total number of rows across all ranks. + /// @param stored_indices Mesh-file stored indices for the locally owned rows, in + /// the order they should appear in the PETSc vector. + void buildPermutation(PetscInt nlocal, PetscInt nglobal, + const std::vector& stored_indices); + + PetscInt global_size{0}; ///< Total global row count. + PetscInt row_start{0}; ///< First global PETSc row owned locally. + PetscInt row_end{0}; ///< One-past-last global PETSc row owned locally. + std::vector local_stored_indices; ///< Stored indices in local PETSc row order. + std::unordered_map local_stored_to_petsc; ///< Stored → PETSc lookup. + Mat mat_stored_to_petsc{nullptr}; ///< Permutation: stored numbering → PETSc ordering. + Mat mat_petsc_to_stored{nullptr}; ///< Permutation: PETSc ordering → stored numbering. +}; + +/// @brief Cell-space index mapping that bridges stored cell numbers, PETSc row +/// ownership, and BOUT++ Field3D index arithmetic. +/// +/// The cell space C contains five disjoint subsets of cells, appended in order: +/// 1. Evolving interior cells (the primary unknowns). +/// 2. Inner radial X-boundary cells. +/// 3. Outer radial X-boundary cells. +/// 4. Forward parallel boundary cells (one virtual cell per interior source cell +/// whose forward map intersects the radial boundary). +/// 5. Backward parallel boundary cells (analogous for the backward map). +/// +/// The local PETSc vector stores these in the same order. The @c mapLocal* methods +/// iterate over subsets 1–3 (interior + X-boundary) while @c mapLocalYup and +/// @c mapLocalYdown cover subsets 4 and 5 respectively. +/// +/// @c mapOwnedInteriorCells iterates only the evolving interior cells (subset 1) and +/// provides global PETSc row indices, making it suitable for assembling matrices. +class PetscCellMapping : public PetscIndexMapping { +public: + /// @brief Construct the cell-space mapping from mesh metadata fields. + /// + /// Populates all five cell subsets and calls @ref buildPermutation with the + /// concatenated list of stored indices. + /// + /// @param cell_number Field3D of stored cell numbers for interior and + /// X-boundary cells. Negative entries are skipped. + /// @param forward_cell_number Field3D of stored cell numbers for forward parallel + /// boundary virtual cells. Negative entries are skipped. + /// @param backward_cell_number Field3D of stored cell numbers for backward parallel + /// boundary virtual cells. Negative entries are skipped. + /// @param total_cells Total number of cells in the global cell space, + /// as computed by the Python weights module. + PetscCellMapping(const Field3D& cell_number, const Field3D& forward_cell_number, + const Field3D& backward_cell_number, int total_cells); + + /// @brief Build the region of evolving interior cells from a cell-number field. + /// + /// Iterates RGN_NOBNDRY and retains indices where @p cell_number is non-negative. + /// + /// @param cell_number Field3D of stored cell numbers. + /// @returns Region containing the Ind3D indices of all valid interior cells. + static Region create_region(const Field3D& cell_number); + + /// @brief Build the region of inner radial X-boundary cells. + /// + /// Iterates RGN_INNER_X and retains indices where @p cell_number is non-negative. + /// + /// @param cell_number Field3D of stored cell numbers. + /// @returns Region containing the Ind3D indices of valid inner-X boundary cells. + static Region create_region_xin(const Field3D& cell_number); + + /// @brief Build the region of outer radial X-boundary cells. + /// + /// Iterates RGN_OUTER_X and retains indices where @p cell_number is non-negative. + /// + /// @param cell_number Field3D of stored cell numbers. + /// @returns Region containing the Ind3D indices of valid outer-X boundary cells. + static Region create_region_xout(const Field3D& cell_number); + + /// @brief Iterate over locally owned evolving interior cells, providing global + /// PETSc row indices. + /// + /// Calls @p func(petsc_row, field_index, stored_cell_number) for each cell in + /// the evolving region that is owned by this MPI rank. @p petsc_row is the global + /// PETSc index and is suitable for use as a row or column argument to MatSetValues. + /// + /// @tparam Function Callable with signature + /// (PetscInt row, Ind3D i, int stored) -> void. + /// @param func Callback invoked for each owned interior cell. + template + void mapOwnedInteriorCells(Function func) const { + PetscInt row = rowStart(); + BOUT_FOR_SERIAL(i, evolving_region) { + func(row, i, ROUND(cell_number[i])); + ++row; + } + } + + /// @brief Iterate over all locally stored field cells (interior + X-boundary), + /// providing local PETSc vector indices. + /// + /// Covers the evolving interior region, inner-X boundary, and outer-X boundary + /// in that order. The local index starts at zero on each rank and increases + /// contiguously. Use this to read from or write to a local PETSc Vec array + /// obtained via VecGetArray. + /// + /// @tparam Function Callable with signature + /// (PetscInt local_row, Ind3D i) -> void. + /// @param func Callback invoked for each locally stored field cell. + template + void mapLocalField(Function func) const { + PetscInt row = 0; + const std::vector>> regions = { + evolving_region, xin_region, xout_region}; + for (const auto& region : regions) { + BOUT_FOR_SERIAL(i, region.get()) { + func(row, i); + ++row; + } + } + } + + /// @brief Iterate over locally stored forward parallel boundary cells, providing + /// local PETSc vector indices. + /// + /// The local index starts immediately after the field cells covered by + /// @ref mapLocalField. The Field3D index passed to @p func is the y+1 neighbour + /// of the source cell (i.yp()), matching how field.yup() is indexed. + /// + /// @tparam Function Callable with signature + /// (PetscInt local_row, Ind3D i_yp) -> void. + /// @param func Callback invoked for each locally stored forward boundary cell. + template + void mapLocalYup(Function func) const { + PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); + BOUT_FOR_SERIAL(i, yup_region) { + func(row, i.yp()); + ++row; + } + } + + /// @brief Iterate over locally stored backward parallel boundary cells, providing + /// local PETSc vector indices. + /// + /// The local index starts immediately after the forward boundary cells covered by + /// @ref mapLocalYup. The Field3D index passed to @p func is the y-1 neighbour + /// of the source cell (i.ym()), matching how field.ydown() is indexed. + /// + /// @tparam Function Callable with signature + /// (PetscInt local_row, Ind3D i_ym) -> void. + /// @param func Callback invoked for each locally stored backward boundary cell. + template + void mapLocalYdown(Function func) const { + PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() + + yup_region.size(); + BOUT_FOR_SERIAL(i, ydown_region) { + func(row, i.ym()); + ++row; + } + } + +private: + Field3D cell_number; ///< Stored cell numbers for interior/X-boundary cells. + Field3D + forward_cell_number; ///< Stored cell numbers for forward boundary virtual cells. + Field3D + backward_cell_number; ///< Stored cell numbers for backward boundary virtual cells. + Region evolving_region; ///< Interior evolving cells. + Region xin_region; ///< Inner radial X-boundary cells. + Region xout_region; ///< Outer radial X-boundary cells. + Region yup_region; ///< Forward parallel boundary virtual cells. + Region ydown_region; ///< Backward parallel boundary virtual cells. +}; + +/// @brief Index mapping for forward (L+) or backward (L-) leg spaces. +/// +/// Each interior source cell contributes one or two leg rows depending on whether +/// its mapped sub-samples are entirely interior, entirely boundary, or mixed. +/// The leg indices are assigned by the Python weights module and stored in the mesh +/// file as Field3D arrays (e.g. forward_leg_interior_number). +/// +/// On construction, duplicate indices are removed and the remaining indices are +/// sorted before being passed to @ref buildPermutation. +class PetscLegMapping : public PetscIndexMapping { +public: + /// @brief Construct a leg-space mapping. + /// + /// @param total_legs Global number of leg rows in this space. + /// @param local_leg_indices Stored leg indices owned by this MPI rank. Duplicates + /// are removed internally before building the permutation. + PetscLegMapping(int total_legs, std::vector local_leg_indices); +}; + +/// @brief Shared-pointer alias for a const PetscCellMapping. +using PetscCellMappingPtr = std::shared_ptr; + +/// @brief Shared-pointer alias for a const PetscLegMapping. +using PetscLegMappingPtr = std::shared_ptr; + +/// @brief Type-safe wrapper around a PETSc Vec belonging to a particular vector space. +/// +/// The compile-time @p SpaceTag (one of @ref CellSpaceTag, @ref ForwardLegSpaceTag, +/// @ref BackwardLegSpaceTag) prevents accidental arithmetic between vectors from +/// incompatible spaces. All binary operations assert that both operands share the +/// same mapping pointer, catching dimension mismatches at runtime. +/// +/// Vectors are non-copyable but movable, matching the ownership semantics of the +/// underlying PETSc Vec. +/// +/// @tparam SpaceTag One of CellSpaceTag, ForwardLegSpaceTag, BackwardLegSpaceTag. +template +class PetscSpaceVector { + using UniqueVec = bout::petsc::UniqueVec; + +public: + /// @brief Default constructor; produces an empty (unusable) vector. + PetscSpaceVector() = default; + + /// @brief Construct a zero-initialised vector for the given space. + /// + /// Allocates a PETSc MPI vector sized according to @p mapping->localSize() and + /// @p mapping->globalSize(). + /// + /// @param mapping Shared ownership of the space mapping. + explicit PetscSpaceVector(std::shared_ptr mapping) + : mapping(std::move(mapping)), + vec(createVec(this->mapping->localSize(), this->mapping->globalSize())) {} + + /// @brief Construct from an existing mapping and a pre-built PETSc Vec. + /// + /// Takes ownership of @p vec. + /// + /// @param mapping Shared ownership of the space mapping. + /// @param vec An already-assembled PETSc Vec. Ownership is transferred. + PetscSpaceVector(std::shared_ptr mapping, UniqueVec vec) + : mapping(std::move(mapping)), vec(std::move(vec)) {} + + PetscSpaceVector(const PetscSpaceVector&) = delete; + PetscSpaceVector& operator=(const PetscSpaceVector&) = delete; + PetscSpaceVector(PetscSpaceVector&&) noexcept = default; + PetscSpaceVector& operator=(PetscSpaceVector&&) noexcept = default; + + /// @brief Return the underlying PETSc Vec handle (non-owning). + Vec raw() const { return *this->vec; } + + /// @brief Return the shared space mapping. + const std::shared_ptr& getMapping() const { return mapping; } + + /// @brief Return a deep copy of this vector. + /// + /// Allocates a new PETSc Vec with the same size and mapping, then copies values. + PetscSpaceVector copy() const { + UniqueVec out{new Vec{nullptr}}; + BOUT_DO_PETSC(VecDuplicate(*this->vec, out.get())); + BOUT_DO_PETSC(VecCopy(*this->vec, *out.get())); + return PetscSpaceVector(this->mapping, std::move(out)); + } + + /// @brief Return a new vector equal to @c (*this) * @p scalar. + /// + /// @param scalar Scalar multiplier. + PetscSpaceVector operator*(BoutReal scalar) const& { + auto out = this->copy(); + BOUT_DO_PETSC(VecScale(out.raw(), scalar)); + return out; + } + + /// Multiply by a scalar, modifying a temporary in-place + PetscSpaceVector operator*(BoutReal scalar) && { + BOUT_DO_PETSC(VecScale(this->raw(), scalar)); + return std::move(*this); + } + + /// @brief Return a new vector equal to @c (*this) + @p rhs. + /// + /// @param rhs Must share the same mapping as @c *this. + PetscSpaceVector operator+(const PetscSpaceVector& rhs) const& { + ASSERT0(mapping == rhs.mapping); + auto out = this->copy(); + BOUT_DO_PETSC(VecAXPY(out.raw(), 1.0, rhs.raw())); + return out; + } + + PetscSpaceVector operator+(const PetscSpaceVector& rhs) && { + ASSERT0(mapping == rhs.mapping); + BOUT_DO_PETSC(VecAXPY(this->raw(), 1.0, rhs.raw())); + return std::move(*this); + } + + /// @brief Return a new vector equal to @c (*this) - @p rhs. + /// + /// @param rhs Must share the same mapping as @c *this. + PetscSpaceVector operator-(const PetscSpaceVector& rhs) const& { + ASSERT0(mapping == rhs.mapping); + auto out = this->copy(); + BOUT_DO_PETSC(VecAXPY(out.raw(), -1.0, rhs.raw())); + return out; + } + + PetscSpaceVector operator-(const PetscSpaceVector& rhs) && { + ASSERT0(mapping == rhs.mapping); + BOUT_DO_PETSC(VecAXPY(this->raw(), -1.0, rhs.raw())); + return std::move(*this); + } + + /// @brief Return a new vector whose entries are the element-wise product + /// @c lhs[i] * rhs[i]. + /// + /// @param lhs Left operand; must share a mapping with @p rhs. + /// @param rhs Right operand. + static PetscSpaceVector pointwiseMultiply(const PetscSpaceVector& lhs, + const PetscSpaceVector& rhs) { + ASSERT0(lhs.mapping == rhs.mapping); + auto out = lhs.copy(); + BOUT_DO_PETSC(VecPointwiseMult(out.raw(), lhs.raw(), rhs.raw())); + return out; + } + + /// @brief Return a new vector whose entries are @c 1 / in[i]. + /// + /// @param in Input vector. No entry may be zero. + static PetscSpaceVector reciprocal(const PetscSpaceVector& in) { + auto out = in.copy(); + BOUT_DO_PETSC(VecReciprocal(out.raw())); + return out; + } + +private: + /// @brief Allocate a new zeroed PETSc MPI vector of the given size. + /// + /// @param local_size Number of locally owned entries. + /// @param global_size Total number of entries across all ranks. + static UniqueVec createVec(PetscInt local_size, PetscInt global_size) { + UniqueVec out{new Vec{nullptr}}; + BOUT_DO_PETSC(VecCreate(BoutComm::get(), out.get())); + BOUT_DO_PETSC(VecSetType(*out, VECMPI)); + BOUT_DO_PETSC(VecSetSizes(*out, local_size, global_size)); + return out; + } + + std::shared_ptr mapping; ///< Space mapping (shared ownership). + UniqueVec vec; ///< Underlying PETSc Vec (owned). +}; + +// Tagged vector aliases — these provide compile-time checks of compatible operations. +/// @brief PETSc vector in the cell space C. +using PetscCellVector = PetscSpaceVector; +/// @brief PETSc vector in the forward leg space L+. +using PetscForwardLegVector = PetscSpaceVector; +/// @brief PETSc vector in the backward leg space L-. +using PetscBackwardLegVector = PetscSpaceVector; + +/// @brief Populate a PetscCellVector from a Field3D and its parallel slices. +/// +/// The cell vector is filled in the order defined by @ref PetscCellMapping: +/// evolving interior cells (from @p field), inner/outer X-boundary cells +/// (from @p field), forward boundary virtual cells (from @p field.yup()), and +/// backward boundary virtual cells (from @p field.ydown()). +/// +/// @param mapping Shared cell-space mapping. +/// @param field Source field; must have parallel slices split and allocated. +/// @returns A fully assembled PetscCellVector. +PetscCellVector makePetscCellVector(const PetscCellMappingPtr& mapping, + const Field3D& field); + +/// @brief Extract a Field3D (with parallel slices) from a PetscCellVector. +/// +/// Reverses @ref makePetscCellVector: interior and X-boundary values are written to +/// @c result, forward boundary values to @c result.yup(), and backward boundary +/// values to @c result.ydown(). +/// +/// @param vec Source vector; must use a PetscCellMapping. +/// @param prototype Field3D whose mesh, coordinates, and metadata are copied to +/// initialise the result. +/// @returns A newly allocated Field3D with split parallel slices. +Field3D toField3D(const PetscCellVector& vec, const Field3D& prototype); + +/// @brief Type-safe sparse linear operator mapping from one space to another. +/// +/// Wraps a PETSc Mat and enforces compatible space pairings at compile time via +/// the @p OutSpaceTag and @p InSpaceTag template parameters. The underlying matrix +/// stores values in PETSc's internal ordering; translation from mesh-file stored +/// column indices to PETSc column indices is performed during construction via a +/// post-multiplication by @ref PetscIndexMapping::getPetscToStored. +/// +/// Operators are non-copyable but movable. Binary operations (addition, subtraction, +/// scalar multiplication, composition) return new operators without modifying +/// their operands, unless the rvalue-qualified @c operator*(BoutReal) overload is +/// selected, which scales in place. +/// +/// @tparam OutSpaceTag Tag for the output (row) space. +/// @tparam InSpaceTag Tag for the input (column) space. +template +class PetscOperator { + using UniqueMat = bout::petsc::UniqueMat; + +public: + /// @brief Default constructor; produces an empty (unusable) operator. + PetscOperator() = default; + ~PetscOperator() = default; + + /// @brief Construct an operator from a CSR-format sparse matrix stored in the mesh. + /// + /// Assembles the PETSc matrix from stored (non-PETSc) row and column indices by + /// calling @ref assembleFromCSR. + /// + /// @param out_mapping Shared mapping for the output (row) space. + /// @param in_mapping Shared mapping for the input (column) space. + /// @param rows CSR row-pointer array (size = number_of_rows + 1). + /// @param cols CSR column-index array, using stored (mesh-file) indices. + /// @param weights CSR nonzero-value array. + PetscOperator(std::shared_ptr out_mapping, + std::shared_ptr in_mapping, Array rows, + Array cols, Array weights) + : out_mapping(std::move(out_mapping)), in_mapping(std::move(in_mapping)) { + assembleFromCSR(rows, cols, weights); + } + + /// @brief Construct an operator from an already-assembled PETSc Mat. + /// + /// Takes ownership of @p mat. + /// + /// @param out_mapping Shared mapping for the output (row) space. + /// @param in_mapping Shared mapping for the input (column) space. + /// @param mat Pre-assembled PETSc Mat. Ownership is transferred. + PetscOperator(std::shared_ptr out_mapping, + std::shared_ptr in_mapping, UniqueMat mat) + : out_mapping(std::move(out_mapping)), in_mapping(std::move(in_mapping)), + mat_operator(std::move(mat)) {} + + // Operators are non-copyable (PETSc matrices are expensive to copy) but movable. + PetscOperator(const PetscOperator&) = delete; + PetscOperator& operator=(const PetscOperator&) = delete; + PetscOperator(PetscOperator&&) noexcept = default; + PetscOperator& operator=(PetscOperator&&) noexcept = default; + + /// @brief Apply the operator to a typed vector: @c result = A * rhs. + /// + /// @param rhs Input vector; must share the same mapping as this operator's + /// input space. + /// @returns A new output-space vector holding the result. + PetscSpaceVector + operator()(const PetscSpaceVector& rhs) const { + ASSERT0(in_mapping == rhs.getMapping()); + PetscSpaceVector result(out_mapping); + BOUT_DO_PETSC(MatMult(*this->mat_operator, rhs.raw(), result.raw())); + return result; + } + + /// @brief Apply the operator to a Field3D: converts to a cell vector then applies. + /// + /// Only available when @p InSpaceTag is @ref CellSpaceTag and @p OutSpaceTag is + /// not @ref CellSpaceTag. The Field3D is packed into a PetscCellVector via + /// @ref makePetscCellVector before applying the matrix. + /// + /// @param rhs Source field; must have parallel slices split and allocated. + /// @returns A new output-space vector holding the result. + template + && !std::is_same_v>> + PetscSpaceVector operator()(const Field3D& rhs) const { + auto in = makePetscCellVector( + std::static_pointer_cast(in_mapping), rhs); + return (*this)(in); + } + + /// Apply operator to a Field3D, returning a Field3D. + /// Enable if both input and output mapppings are in Cell Space + template + && std::is_same_v>> + Field3D operator()(const Field3D& rhs) const { + auto in = makePetscCellVector( + std::static_pointer_cast(in_mapping), rhs); + return toField3D((*this)(in), rhs); + } + + /// @brief Return the transpose operator mapping from OutSpaceTag to InSpaceTag. + /// + /// Allocates a new PETSc Mat via MatTranspose (MAT_INITIAL_MATRIX). The mappings + /// are swapped: the transposed operator's input mapping is this operator's output + /// mapping, and vice versa. + PetscOperator transpose() const { + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatTranspose(*mat_operator, MAT_INITIAL_MATRIX, out.get())); + return PetscOperator(in_mapping, out_mapping, + std::move(out)); + } + + /// @brief Return a new operator equal to @c (*this) + @p rhs. + /// + /// Both operands must share the same input and output mappings. Uses + /// UNKNOWN_NONZERO_PATTERN, so the sparsity structure of the result may differ + /// from either operand. + /// + /// @param rhs Operator to add; must have the same in/out mappings. + PetscOperator operator+(const PetscOperator& rhs) const { + ASSERT0(out_mapping == rhs.out_mapping); + ASSERT0(in_mapping == rhs.in_mapping); + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, out.get())); + BOUT_DO_PETSC(MatAXPY(*out, 1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); + return PetscOperator(out_mapping, in_mapping, std::move(out)); + } + + /// @brief Return a new operator equal to @c (*this) - @p rhs. + /// + /// Both operands must share the same input and output mappings. Uses + /// UNKNOWN_NONZERO_PATTERN; see @ref operator+. + /// + /// @param rhs Operator to subtract; must have the same in/out mappings. + PetscOperator operator-(const PetscOperator& rhs) const { + ASSERT0(out_mapping == rhs.out_mapping); + ASSERT0(in_mapping == rhs.in_mapping); + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, out.get())); + BOUT_DO_PETSC(MatAXPY(*out, -1.0, *rhs.mat_operator, UNKNOWN_NONZERO_PATTERN)); + return PetscOperator(out_mapping, in_mapping, std::move(out)); + } + + /// @brief Return a new operator equal to @c (*this) * @p scalar (lvalue overload). + /// + /// Duplicates the matrix before scaling, leaving @c *this unchanged. + /// + /// @param scalar Scalar multiplier. + PetscOperator operator*(BoutReal scalar) const& { + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatDuplicate(*this->mat_operator, MAT_COPY_VALUES, out.get())); + BOUT_DO_PETSC(MatScale(*out, scalar)); + return PetscOperator(out_mapping, in_mapping, std::move(out)); + } + + /// @brief Scale this operator in place and return it (rvalue overload). + /// + /// More efficient than the lvalue overload when @c *this is a temporary, + /// as no matrix duplication is needed. + /// + /// @param scalar Scalar multiplier. + PetscOperator operator*(BoutReal scalar) && { + BOUT_DO_PETSC(MatScale(*this->mat_operator, scalar)); + return std::move(*this); + } + + /// @brief Return the underlying PETSc Mat handle (non-owning). + /// + /// Intended for use by the free @c operator* (matrix composition) and other + /// low-level callers that need direct PETSc access. + Mat raw() const { return *this->mat_operator; } + + /// @brief Construct a diagonal operator from a vector of diagonal entries. + /// + /// Only available when @p OutSpaceTag == @p InSpaceTag (square, same-space + /// operator). Builds an MPIAIJ matrix with one nonzero per row via + /// MatDiagonalSet. + /// + /// @param mapping Shared mapping for both the row and column space. + /// @param diag Vector of diagonal values; must share @p mapping. + /// @returns A new diagonal operator. + static PetscOperator diagonal(std::shared_ptr mapping, + const PetscSpaceVector& diag) { + static_assert(std::is_same_v, + "Diagonal operators require the same input and output space"); + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatCreate(BoutComm::get(), out.get())); + BOUT_DO_PETSC(MatSetType(*out, MATMPIAIJ)); + BOUT_DO_PETSC(MatSetSizes(*out, mapping->localSize(), mapping->localSize(), + mapping->globalSize(), mapping->globalSize())); + // Note: Off-diagonal blocks are empty because this is a square matrix + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(*out, 1, nullptr, 0, nullptr)); + BOUT_DO_PETSC(MatDiagonalSet(*out, diag.raw(), INSERT_VALUES)); + BOUT_DO_PETSC(MatAssemblyBegin(*out, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(*out, MAT_FINAL_ASSEMBLY)); + return PetscOperator(mapping, mapping, std::move(out)); + } + +private: + /// @brief Assemble the operator matrix from a mesh-file CSR representation. + /// + /// Iterates over the locally owned output rows (identified via + /// @c out_mapping->localStoredIndices()), looks up each row's nonzeros in the + /// CSR arrays, and inserts them using global PETSc row indices. Column indices + /// are given in stored (mesh-file) numbering; the final matrix is obtained by + /// post-multiplying the assembled temporary by + /// @c in_mapping->getPetscToStored() to convert to PETSc column ordering. + /// + /// @param rows CSR row-pointer array (size = total_stored_rows + 1). + /// @param cols CSR column-index array in stored (mesh-file) numbering. + /// @param weights CSR nonzero-value array. + void assembleFromCSR(const Array& rows, const Array& cols, + const Array& weights) { + UniqueMat temp{new Mat{nullptr}}; + BOUT_DO_PETSC(MatCreate(BoutComm::get(), temp.get())); + BOUT_DO_PETSC(MatSetType(*temp, MATMPIAIJ)); + BOUT_DO_PETSC(MatSetSizes(*temp, out_mapping->localSize(), PETSC_DECIDE, + out_mapping->globalSize(), in_mapping->globalSize())); + + const auto& local_rows = out_mapping->localStoredIndices(); + for (PetscInt local_row = 0; local_row < static_cast(local_rows.size()); + ++local_row) { + const int stored_row = local_rows[local_row]; + if (stored_row < 0 || stored_row + 1 >= static_cast(rows.size())) { + continue; + } + const int start = rows[stored_row]; + const int end = rows[stored_row + 1]; + if (end <= start) { + continue; + } + const PetscInt petsc_row = out_mapping->rowStart() + local_row; + BOUT_DO_PETSC(MatSetValues(*temp, 1, &petsc_row, end - start, &cols[start], + &weights[start], INSERT_VALUES)); + } + BOUT_DO_PETSC(MatAssemblyBegin(*temp, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(*temp, MAT_FINAL_ASSEMBLY)); + + // Post-multiply by the stored→PETSc permutation to convert column numbering. + UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC(MatMatMult(*temp, in_mapping->getPetscToStored(), MAT_INITIAL_MATRIX, + PETSC_DETERMINE, out.get())); + mat_operator = std::move(out); + } + + std::shared_ptr out_mapping; ///< Row-space mapping. + std::shared_ptr in_mapping; ///< Column-space mapping. + UniqueMat mat_operator{new Mat{nullptr}}; ///< Underlying PETSc Mat. + + template + friend PetscOperator operator*(const PetscOperator& lhs, + const PetscOperator& rhs); +}; + +/// @brief Compose two compatible operators: @c result = lhs * rhs. +/// +/// The output space of @p rhs must match the input space of @p lhs, enforced at +/// compile time by the @p MidSpaceTag parameter. The resulting operator maps from +/// @p rhs's input space to @p lhs's output space. +/// +/// @tparam OutSpaceTag Output space of @p lhs and of the result. +/// @tparam MidSpaceTag Input space of @p lhs and output space of @p rhs. +/// @tparam InSpaceTag Input space of @p rhs and of the result. +/// @param lhs Left (outer) operator. +/// @param rhs Right (inner) operator. +/// @returns A new operator representing the composition. +template +PetscOperator +operator*(const PetscOperator& lhs, + const PetscOperator& rhs) { + ASSERT0(lhs.in_mapping == rhs.out_mapping); + typename PetscOperator::UniqueMat out{new Mat{nullptr}}; + BOUT_DO_PETSC( + MatMatMult(lhs.raw(), rhs.raw(), MAT_INITIAL_MATRIX, PETSC_DETERMINE, out.get())); + return PetscOperator(lhs.out_mapping, rhs.in_mapping, + std::move(out)); +} + +/// @brief Operator mapping from cell space C to forward leg space L+. +using PetscForwardOperator = PetscOperator; +/// @brief Operator mapping from cell space C to backward leg space L-. +using PetscBackwardOperator = PetscOperator; +/// @brief Operator mapping from forward leg space L+ to cell space C. +using PetscForwardToCellOperator = PetscOperator; +/// @brief Operator mapping from backward leg space L- to cell space C. +using PetscBackwardToCellOperator = PetscOperator; +/// @brief Operator mapping within the cell space C. +using PetscCellOperator = PetscOperator; +/// @brief Operator mapping within the forward leg space L+. +using PetscForwardLegOperator = PetscOperator; +/// @brief Operator mapping within the backward leg space L-. +using PetscBackwardLegOperator = PetscOperator; + +/// @brief Factory that constructs cell-space and leg-space operators from mesh metadata. +/// +/// On construction, reads all required fields and integer scalars from the mesh and +/// sets up the @ref PetscCellMapping, forward @ref PetscLegMapping, and backward +/// @ref PetscLegMapping. Individual primitive operators (interpolation, injection, +/// leg weights) are constructed on demand. The high-level set of parallel operators +/// needed for Support-Operator-Method diffusion is assembled by @ref getParallel. +class PetscOperators { +public: + /// @brief Construct the operator factory from a BOUT++ mesh. + /// + /// Reads the following fields from @p mesh: + /// - @c cell_number, @c forward_cell_number, @c backward_cell_number (Field3D) + /// - @c total_cells, @c n_forward_legs, @c n_backward_legs (int) + /// - @c forward_leg_interior_number, @c forward_leg_boundary_number (Field3D) + /// - @c backward_leg_interior_number, @c backward_leg_boundary_number (Field3D) + /// + /// @param mesh Non-owning pointer to the BOUT++ mesh. Must outlive this object. + /// @throws BoutException if any required field or scalar is missing from the mesh. + explicit PetscOperators(Mesh* mesh); + + /// @brief Build the forward interpolation operator F: L+ ← C. + /// + /// Maps cell-space values to forward leg space using the CSR arrays stored in the + /// mesh (@c forward_rows, @c forward_columns, @c forward_weights). + PetscForwardOperator forward() const; + + /// @brief Build the backward interpolation operator B: L- ← C. + /// + /// Maps cell-space values to backward leg space using the CSR arrays stored in the + /// mesh (@c backward_rows, @c backward_columns, @c backward_weights). + PetscBackwardOperator backward() const; + + /// @brief Build the forward injection operator I+: L+ ← C. + /// + /// Places the value of each interior source cell into its corresponding forward leg + /// row(s), with weight 1. Boundary leg rows point to the virtual boundary cell in C. + PetscForwardOperator injectForward() const; + + /// @brief Build the backward injection operator I-: L- ← C. + /// + /// Analogous to @ref injectForward for the backward leg space. + PetscBackwardOperator injectBackward() const; + + /// @brief Load the forward leg-weight vector w+: L+. + /// + /// Each entry holds the quadrature weight of the corresponding forward leg row, + /// equal to the interior or boundary fraction of that leg's sub-sample integral. + PetscForwardLegVector forwardLegWeights() const; + + /// @brief Load the backward leg-weight vector w-: L-. + /// + /// Analogous to @ref forwardLegWeights for the backward leg space. + PetscBackwardLegVector backwardLegWeights() const; + + /// @brief Build a diagonal cell-space operator from a Field3D of values. + /// + /// @param f Field3D of diagonal entries; must have parallel slices allocated. + /// @returns A diagonal PetscCellOperator. + PetscCellOperator diagonal(const Field3D& f) const; + + /// @brief Build a diagonal forward-leg-space operator from a leg vector. + /// + /// @param v Forward leg vector of diagonal entries. + /// @returns A diagonal PetscForwardLegOperator. + PetscForwardLegOperator diagonal(const PetscForwardLegVector& v) const; + + /// @brief Build a diagonal backward-leg-space operator from a leg vector. + /// + /// @param v Backward leg vector of diagonal entries. + /// @returns A diagonal PetscBackwardLegOperator. + PetscBackwardLegOperator diagonal(const PetscBackwardLegVector& v) const; + + /// @brief Collection of parallel differential operators assembled for the + /// Support Operator Method (SOM). + /// + /// All operators are constructed together in @ref getParallel so that intermediate + /// quantities (leg weights, interpolated dl, dV) are shared correctly. + /// + /// Naming convention: the subscript on Grad/Div/Restrict refers to which side of + /// the cell face the operator acts on, not the leg space it lives in. + /// - @c Grad_plus / @c Div_plus: forward-side half-step operators (in L-). + /// - @c Grad_minus / @c Div_minus: backward-side half-step operators (in L+). + /// - @c Restrict_minus = I+^T * W+: weighted restriction from L+ back to C, + /// paired with the minus-side gradient. + /// - @c Restrict_plus = I-^T * W-: weighted restriction from L- back to C, + /// paired with the plus-side gradient. + struct Parallel { + /// Central parallel gradient: C ← C. + PetscCellOperator Grad_par; + /// Support-operator parallel divergence paired with Grad_par: C ← C. + PetscCellOperator Div_par; + /// Combined SOM parallel Laplacian (Div_par ∘ Grad_par): C ← C. + PetscCellOperator Div_par_Grad_par; + /// Cell volume J * dx * dy * dz. + Field3D dV; + + /// Half-step gradient on the backward (y-1) side: L- ← C. + PetscBackwardOperator Grad_minus; + /// Half-step gradient on the forward (y+1) side: L+ ← C. + PetscForwardOperator Grad_plus; + /// SOM divergence paired with Grad_plus: C ← L+. + PetscForwardToCellOperator Div_minus; + /// SOM divergence paired with Grad_minus: C ← L-. + PetscBackwardToCellOperator Div_plus; + + /// Raw backward interpolation operator B: L- ← C. + PetscBackwardOperator Backward; + /// Raw forward interpolation operator F: L+ ← C. + PetscForwardOperator Forward; + + /// Injection operator into backward leg space I-: L- ← C. + PetscBackwardOperator Inject_minus; + /// Injection operator into forward leg space I+: L+ ← C. + PetscForwardOperator Inject_plus; + + /// Average interpolation to backward leg positions: L- ← C. + PetscBackwardOperator Interp_minus; + /// Average interpolation to forward leg positions: L+ ← C. + PetscForwardOperator Interp_plus; + + /// Weighted restriction from L+ to C: I+^T * W+. + PetscForwardToCellOperator Restrict_minus; + /// Weighted restriction from L- to C: I-^T * W-. + PetscBackwardToCellOperator Restrict_plus; + + /// @brief Compute the SOM parallel diffusion Div_par(K * Grad_par(f)). + /// + /// Evaluates the half-step fluxes separately on each side, multiplies by the + /// interpolated diffusion coefficient K, and averages the two divergence + /// contributions. This is the primary user-facing SOM operator. + /// + /// @param K Diffusion coefficient field; must have parallel slices allocated. + /// @param f Field to differentiate; must have parallel slices allocated. + /// @returns Field3D containing Div_par(K * Grad_par(f)) at interior cells. + Field3D Div_par_K_Grad_par(const Field3D& K, const Field3D& f) const; + }; + + /// @brief Assemble and return the full set of parallel SOM operators. + /// + /// Constructs all primitive operators (Forward, Backward, Inject_plus, + /// Inject_minus, leg weights), derives the interpolation, gradient, restriction, + /// and divergence operators, and returns them packaged in a @ref Parallel struct. + /// Intermediate Field3D quantities (dl, dV) are computed from the mesh + /// coordinates and boundary conditions are applied before use. + /// + /// @returns A @ref Parallel struct containing all assembled operators. + Parallel getParallel() const; + +private: + /// @brief Read a Field3D from the mesh, throwing if the field is absent. + /// + /// @param mesh Non-owning mesh pointer. + /// @param name Field name as stored in the mesh file. + /// @returns The requested Field3D. + /// @throws BoutException if the field is not found. + static Field3D meshGetField3D(Mesh* mesh, const std::string& name); + + /// @brief Read an Array from the mesh, throwing if absent. + /// + /// @tparam T Element type of the array. + /// @param name Array name as stored in the mesh file. + /// @returns The requested Array. + /// @throws BoutException if the array is not found. + template + Array meshGetArray(const std::string& name) const { + Array result; + if (mesh->get(result, name) != 0) { + throw BoutException("PetscOperators requires array '{}'", name); + } + return result; + } + + /// @brief Read a scalar integer from the mesh, throwing if absent. + /// + /// @param name Scalar name as stored in the mesh file. + /// @returns The requested integer value. + /// @throws BoutException if the scalar is not found. + int meshGetInt(const std::string& name) const; + + /// @brief Build a leg-space injection operator for the given direction. + /// + /// Iterates owned interior cells and inserts a unit entry for each interior or + /// boundary leg row, mapping to the corresponding cell-space column. Column + /// indices come from @p interior_leg_number and @p boundary_leg_number. + /// + /// @tparam LegTag @ref ForwardLegSpaceTag or @ref BackwardLegSpaceTag. + /// @param interior_leg_number Field3D of stored interior leg row numbers. + /// @param boundary_leg_number Field3D of stored boundary leg row numbers. + /// @param leg_mapping Shared mapping for the leg space. + /// @returns The assembled injection operator LegTag ← CellSpaceTag. + template + PetscOperator + buildInjection(const Field3D& interior_leg_number, const Field3D& boundary_leg_number, + std::shared_ptr leg_mapping) const; + + /// @brief Load a leg-weight vector from a named mesh array. + /// + /// Reads the flat array @p name from the mesh, then scatters values into the + /// local PETSc vector using the mapping's stored index list. Stored indices that + /// are out of range are assigned weight 0. + /// + /// @tparam LegVector @ref PetscForwardLegVector or @ref PetscBackwardLegVector. + /// @param name Array name in the mesh file (e.g. "forward_leg_weights"). + /// @param mapping Shared mapping for the leg space. + /// @returns A fully populated leg-weight vector. + template + LegVector loadLegWeights(const std::string& name, + std::shared_ptr mapping) const { + auto values = meshGetArray(name); + LegVector out(mapping); + PetscScalar* x{nullptr}; + BOUT_DO_PETSC(VecGetArray(out.raw(), &x)); + const auto& local = mapping->localStoredIndices(); + for (PetscInt i = 0; i < static_cast(local.size()); ++i) { + const int stored = local[i]; + x[i] = (stored >= 0 && stored < static_cast(values.size())) ? values[stored] + : 0.0; + } + BOUT_DO_PETSC(VecRestoreArray(out.raw(), &x)); + return out; + } + + Mesh* mesh{nullptr}; ///< Non-owning pointer to the BOUT++ mesh. + PetscCellMappingPtr cell_mapping; ///< Cell-space index mapping. + PetscLegMappingPtr forward_leg_mapping; ///< Forward leg-space index mapping. + PetscLegMappingPtr backward_leg_mapping; ///< Backward leg-space index mapping. + Field3D forward_leg_interior_number; ///< Stored interior leg numbers for L+. + Field3D forward_leg_boundary_number; ///< Stored boundary leg numbers for L+. + Field3D backward_leg_interior_number; ///< Stored interior leg numbers for L-. + Field3D backward_leg_boundary_number; ///< Stored boundary leg numbers for L-. +}; + +#else + +#warning PETSc not available. No PetscOperators. + +#endif diff --git a/include/bout/petsclib.hxx b/include/bout/petsclib.hxx index 41d7618fd2..637943bd12 100644 --- a/include/bout/petsclib.hxx +++ b/include/bout/petsclib.hxx @@ -1,31 +1,31 @@ /*!************************************************************************ * Provides access to the PETSc library, handling initialisation and - * finalisation. + * finalisation. * * Usage * ----- * * #include - * + * * class MyClass { * public: - * + * * private: * PetscLib lib; * }; - * + * * * This will then automatically initialise Petsc the first time an object * is created, and finalise it when the last object is destroyed. * * This header tries to workaround some annoying PETSc features, and * so it *must* be included before *any* PETSc header. - * + * ************************************************************************** - * Copyright 2012 - 2025 BOUT++ contributors + * Copyright 2012 - 2026 BOUT++ contributors * * Contact: Ben Dudson, dudson2@llnl.gov - * + * * This file is part of BOUT++. * * BOUT++ is free software: you can redistribute it and/or modify @@ -48,6 +48,8 @@ #include "bout/build_defines.hxx" +#include + class Options; #if BOUT_HAS_PETSC @@ -66,7 +68,7 @@ class Options; #include "bout/boutexception.hxx" // NOLINTNEXTLINE(cppcoreguidelines-macro-usage) -#define BOUT_DO_PETSC(cmd) PetscLib::assertIerr((cmd), #cmd) +#define BOUT_DO_PETSC(cmd) PetscLib::assertIerr((cmd), #cmd, __FILE__, __LINE__) /*! * Handles initialisation and finalisation of PETSc library. @@ -124,10 +126,17 @@ public: static void cleanup(); static inline void assertIerr(PetscErrorCode ierr, - const std::string& petsc_op = "PETSc operation") { - if (ierr != 0) { - throw BoutException("{:s} failed with {:d}", petsc_op, ierr); + const std::string& petsc_op = "PETSc operation", + const char* file = nullptr, int linenr = 0) { + if (PetscLikely(ierr == PETSC_SUCCESS)) { + return; } + + const char* err_msg = nullptr; + PetscErrorMessage(ierr, &err_msg, nullptr); + throw BoutException("{:s} failed at {}:{} with {:d} [{}]", petsc_op, + (file != nullptr) ? file : "unknown_file", linenr, ierr, + (err_msg != nullptr) ? err_msg : ""); } static BoutException SNESFailure(SNES& snes); diff --git a/include/bout/region.hxx b/include/bout/region.hxx index e00ad6d41d..d423c3ef4c 100644 --- a/include/bout/region.hxx +++ b/include/bout/region.hxx @@ -1,9 +1,9 @@ /************************************************************************** - * Copyright 2010 B.D.Dudson, S.Farley, M.V.Umansky, X.Q.Xu + * Copyright 2010 - 2026 BOUT++ contributors * * Region class by D. Dickinson 2018 * - * Contact: Ben Dudson, bd512@york.ac.uk + * Contact: Ben Dudson, dudson2@llnl.gov * * This file is part of BOUT++. * @@ -170,8 +170,8 @@ struct SpecificInd { int ny = -1, nz = -1; ///< Sizes of y and z dimensions SpecificInd() = default; - SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz){}; - explicit SpecificInd(int i) : ind(i){}; + SpecificInd(int i, int ny, int nz) : ind(i), ny(ny), nz(nz) {}; + explicit SpecificInd(int i) : ind(i) {}; /// Allow explicit conversion to an int explicit operator int() const { return ind; } @@ -491,10 +491,9 @@ template class Region { // Following prevents a Region being created with anything other // than Ind2D, Ind3D or IndPerp as template type - static_assert( - std::is_base_of_v< - Ind2D, T> || std::is_base_of_v || std::is_base_of_v, - "Region must be templated with one of IndPerp, Ind2D or Ind3D"); + static_assert(std::is_base_of_v || std::is_base_of_v + || std::is_base_of_v, + "Region must be templated with one of IndPerp, Ind2D or Ind3D"); public: using data_type = T; @@ -570,7 +569,7 @@ public: }; Region(RegionIndices& indices, int maxregionblocksize = MAXREGIONBLOCKSIZE) - : indices(indices), blocks(getContiguousBlocks(maxregionblocksize)){}; + : indices(indices), blocks(getContiguousBlocks(maxregionblocksize)) {}; // We need to first set the blocks, and only after that call getRegionIndices. // Do not put in the member initialisation diff --git a/src/invert/fft_fftw.cxx b/src/invert/fft_fftw.cxx index 2c633da6d5..2ada5bc98a 100644 --- a/src/invert/fft_fftw.cxx +++ b/src/invert/fft_fftw.cxx @@ -4,9 +4,9 @@ * \brief FFT routines using external libraries * ************************************************************************** - * 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++. * @@ -25,9 +25,12 @@ * **************************************************************************/ -#include "bout/build_defines.hxx" - +#include +#include +#include +#include #include +#include #include #include #include @@ -518,7 +521,7 @@ Array rfft(const Array& in) { ASSERT1(!in.empty()); int size{in.size()}; - Array out{(size / 2) + 1}; + Array out((size / 2) + 1); bout::fft::rfft(in.begin(), size, out.begin()); return out; @@ -528,7 +531,7 @@ Array irfft(const Array& in, int length) { ASSERT1(!in.empty()); ASSERT1(in.size() == (length / 2) + 1); - Array out{length}; + Array out(length); bout::fft::irfft(in.begin(), length, out.begin()); return out; diff --git a/src/mesh/data/gridfromfile.cxx b/src/mesh/data/gridfromfile.cxx index 10f6f6926e..e2e9850612 100644 --- a/src/mesh/data/gridfromfile.cxx +++ b/src/mesh/data/gridfromfile.cxx @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -14,7 +15,9 @@ #include #include +#include #include +#include GridFile::GridFile(std::string gridfilename) : GridDataSource(true), data(bout::OptionsIO::create(gridfilename)->read()), @@ -120,7 +123,7 @@ bool GridFile::get(Mesh* UNUSED(m), BoutReal& rval, const std::string& name, /*! * Reads a 2D, 3D or FieldPerp field variable from a file - * + * * Successfully reads Field2D or FieldPerp if the variable in the file is 0-D or 2-D. * Successfully reads Field3D if the variable in the file is 0-D, 2-D or 3-D. */ @@ -378,7 +381,7 @@ void GridFile::readField(Mesh* m, const std::string& name, int ys, int yd, int n for (int x = xs; x < xs + nx_to_read; ++x) { for (int y = ys; y < ys + ny_to_read; ++y) { - BoutReal const value = full_var(x, y); + const BoutReal value = full_var(x, y); for (int z = 0; z < var.getNz(); z++) { var(x - xs + xd, y - ys + yd, z) = value; } @@ -466,6 +469,22 @@ void GridFile::readField(Mesh* m, const std::string& name, int UNUSED(ys), int U } } +bool GridFile::get(Array& var, const std::string& name) { + if (not data.isSet(name)) { + return false; + } + var = data[name].as>(); + return true; +} + +bool GridFile::get(Array& var, const std::string& name) { + if (not data.isSet(name)) { + return false; + } + var = data[name].as>(); + return true; +} + bool GridFile::get([[maybe_unused]] Mesh* m, [[maybe_unused]] std::vector& var, [[maybe_unused]] const std::string& name, [[maybe_unused]] int len, [[maybe_unused]] int offset, diff --git a/src/mesh/mesh.cxx b/src/mesh/mesh.cxx index 9a98b5848a..03fa3934c1 100644 --- a/src/mesh/mesh.cxx +++ b/src/mesh/mesh.cxx @@ -179,13 +179,21 @@ int Mesh::get(bool& bval, const std::string& name, bool def) { if (source == nullptr) { warn_default_used(def, name); bval = def; - return true; + return 1; } int bval_as_int = 0; - bool success = source->get(this, bval_as_int, name, def); + const bool success = source->get(this, bval_as_int, name, int(def)); bval = bool(bval_as_int); - return !success; + return success ? 0 : 1; +} + +int Mesh::get(Array& var, const std::string& name) { + return source->get(var, name) ? 0 : 1; +} + +int Mesh::get(Array& var, const std::string& name) { + return source->get(var, name) ? 0 : 1; } int Mesh::get(Field2D& var, const std::string& name, BoutReal def, bool communicate, diff --git a/src/mesh/petsc_operators.cxx b/src/mesh/petsc_operators.cxx new file mode 100644 index 0000000000..0fe564bcd3 --- /dev/null +++ b/src/mesh/petsc_operators.cxx @@ -0,0 +1,423 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "bout/assert.hxx" +#include "bout/bout_types.hxx" +#include "bout/boutexception.hxx" +#include "bout/field3d.hxx" +#include "bout/output.hxx" +#include "bout/output_bout_types.hxx" +#include "bout/petsc_operators.hxx" +#include "bout/region.hxx" +#include "bout/utils.hxx" + +#include +#include +#include +#include +#include + +PetscIndexMapping::~PetscIndexMapping() { + if (mat_stored_to_petsc != nullptr) { + MatDestroy(&mat_stored_to_petsc); + } + if (mat_petsc_to_stored != nullptr) { + MatDestroy(&mat_petsc_to_stored); + } +} + +PetscInt PetscIndexMapping::storedToPetsc(int stored_index) const { + auto it = local_stored_to_petsc.find(stored_index); + if (it == local_stored_to_petsc.end()) { + throw BoutException("Stored index {} is not owned locally", stored_index); + } + return it->second; +} + +void PetscIndexMapping::buildPermutation(PetscInt nlocal, PetscInt nglobal, + const std::vector& stored_indices) { + global_size = nglobal; + local_stored_indices = stored_indices; + + // Renumbering matrix. + // Maps PETSc row (or column) indices to the global indices stored in + // the mesh. This is needed because the PETSc indices depend on the + // number of processors. + BOUT_DO_PETSC(MatCreate(BoutComm::get(), &mat_stored_to_petsc)); + BOUT_DO_PETSC(MatSetSizes(mat_stored_to_petsc, nlocal, nlocal, nglobal, nglobal)); + BOUT_DO_PETSC(MatSetType(mat_stored_to_petsc, MATMPIAIJ)); + + // Each row will have one non-zero entry, which could be in + // either the "diagonal" or "off-diagonal" block. + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(mat_stored_to_petsc, 1, nullptr, 1, nullptr)); + + // Get the range of global rows owned by this processor + BOUT_DO_PETSC(MatGetOwnershipRange(mat_stored_to_petsc, &row_start, &row_end)); + ASSERT1(row_end - row_start == nlocal); + + const PetscScalar one = 1.0; + for (PetscInt i = 0; i < nlocal; ++i) { + const PetscInt petsc_row = row_start + i; + const PetscInt stored_col = stored_indices[i]; + ASSERT1(stored_col >= 0); + local_stored_to_petsc[static_cast(stored_col)] = petsc_row; + BOUT_DO_PETSC(MatSetValues(mat_stored_to_petsc, 1, &petsc_row, 1, &stored_col, &one, + INSERT_VALUES)); + } + BOUT_DO_PETSC(MatAssemblyBegin(mat_stored_to_petsc, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(mat_stored_to_petsc, MAT_FINAL_ASSEMBLY)); + // Transpose to map petsc indices to stored mesh indices + BOUT_DO_PETSC( + MatTranspose(mat_stored_to_petsc, MAT_INITIAL_MATRIX, &mat_petsc_to_stored)); +} + +Region PetscCellMapping::create_region(const Field3D& cell_number) { + Region::RegionIndices indices; + BOUT_FOR(i, cell_number.getRegion("RGN_NOBNDRY")) { + if (cell_number[i] >= 0) { + indices.push_back(i); + } + } + return Region(indices); +} + +Region PetscCellMapping::create_region_xin(const Field3D& cell_number) { + Region::RegionIndices indices; + const auto& region = cell_number.getRegion("RGN_INNER_X"); + BOUT_FOR(i, region) { + if (cell_number[i] >= 0) { + indices.push_back(i); + } + } + return Region(indices); +} + +Region PetscCellMapping::create_region_xout(const Field3D& cell_number) { + Region::RegionIndices indices; + const auto& region = cell_number.getRegion("RGN_OUTER_X"); + BOUT_FOR(i, region) { + if (cell_number[i] >= 0) { + indices.push_back(i); + } + } + return Region(indices); +} + +PetscCellMapping::PetscCellMapping(const Field3D& cell_number, + const Field3D& forward_cell_number, + const Field3D& backward_cell_number, int total_cells) + : cell_number(cell_number), forward_cell_number(forward_cell_number), + backward_cell_number(backward_cell_number), + evolving_region(create_region(cell_number)), + xin_region(create_region_xin(cell_number)), + xout_region(create_region_xout(cell_number)), + yup_region(create_region(forward_cell_number)), + ydown_region(create_region(backward_cell_number)) { + std::vector local_indices; + const auto push_region = [&](const Region& region, const Field3D& values) { + BOUT_FOR_SERIAL(i, region) { local_indices.push_back(ROUND(values[i])); } + }; + push_region(evolving_region, this->cell_number); + push_region(xin_region, this->cell_number); + push_region(xout_region, this->cell_number); + push_region(yup_region, this->forward_cell_number); + push_region(ydown_region, this->backward_cell_number); + buildPermutation(static_cast(local_indices.size()), total_cells, + local_indices); +} + +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()), + local_leg_indices.end()); + buildPermutation(static_cast(local_leg_indices.size()), total_legs, + local_leg_indices); +} + +PetscCellVector makePetscCellVector(const PetscCellMappingPtr& mapping, + const Field3D& field) { + ASSERT1(field.hasParallelSlices()); + ASSERT1(field.yup().isAllocated()); + ASSERT1(field.ydown().isAllocated()); + PetscCellVector out(mapping); + PetscScalar* x{nullptr}; + BOUT_DO_PETSC(VecGetArray(out.raw(), &x)); + mapping->mapLocalField([&](PetscInt row, const Ind3D& i) { x[row] = field[i]; }); + const auto& yup = field.yup(); + mapping->mapLocalYup([&](PetscInt row, const Ind3D& i) { x[row] = yup[i]; }); + const auto& ydown = field.ydown(); + mapping->mapLocalYdown([&](PetscInt row, const Ind3D& i) { x[row] = ydown[i]; }); + BOUT_DO_PETSC(VecRestoreArray(out.raw(), &x)); + return out; +} + +Field3D toField3D(const PetscCellVector& vec, const Field3D& prototype) { + auto mapping = std::static_pointer_cast(vec.getMapping()); + Field3D result{emptyFrom(prototype)}; + result.splitParallelSlices(); + result.yup().allocate(); + result.ydown().allocate(); + + const PetscScalar* x{nullptr}; + BOUT_DO_PETSC(VecGetArrayRead(vec.raw(), &x)); + mapping->mapLocalField([&](PetscInt row, const Ind3D& i) { result[i] = x[row]; }); + mapping->mapLocalYup([&](PetscInt row, const Ind3D& i) { result.yup()[i] = x[row]; }); + mapping->mapLocalYdown( + [&](PetscInt row, const Ind3D& i) { result.ydown()[i] = x[row]; }); + BOUT_DO_PETSC(VecRestoreArrayRead(vec.raw(), &x)); + + return result; +} + +Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { + Field3D result; + if (mesh->get(result, name) != 0) { + throw BoutException("PetscOperators requires field '{}'", name); + } + return result; +} + +int PetscOperators::meshGetInt(const std::string& name) const { + int result{0}; + if (mesh->get(result, name) != 0) { + throw BoutException("PetscOperators requires int '{}'", name); + } + return result; +} + +PetscOperators::PetscOperators(Mesh* mesh) + : mesh(mesh), + cell_mapping(std::make_shared( + meshGetField3D(mesh, "cell_number"), + meshGetField3D(mesh, "forward_cell_number"), + meshGetField3D(mesh, "backward_cell_number"), meshGetInt("total_cells"))), + forward_leg_interior_number(meshGetField3D(mesh, "forward_leg_interior_number")), + forward_leg_boundary_number(meshGetField3D(mesh, "forward_leg_boundary_number")), + backward_leg_interior_number(meshGetField3D(mesh, "backward_leg_interior_number")), + backward_leg_boundary_number(meshGetField3D(mesh, "backward_leg_boundary_number")) { + std::vector local_forward_legs; + std::vector local_backward_legs; + cell_mapping->mapOwnedInteriorCells([&](PetscInt, const Ind3D& i, int) { + const int f_int = ROUND(forward_leg_interior_number[i]); + const int f_bnd = ROUND(forward_leg_boundary_number[i]); + const int b_int = ROUND(backward_leg_interior_number[i]); + const int b_bnd = ROUND(backward_leg_boundary_number[i]); + if (f_int >= 0) { + local_forward_legs.push_back(f_int); + } + if (f_bnd >= 0) { + local_forward_legs.push_back(f_bnd); + } + if (b_int >= 0) { + local_backward_legs.push_back(b_int); + } + if (b_bnd >= 0) { + local_backward_legs.push_back(b_bnd); + } + }); + + forward_leg_mapping = std::make_shared(meshGetInt("n_forward_legs"), + std::move(local_forward_legs)); + backward_leg_mapping = std::make_shared( + meshGetInt("n_backward_legs"), std::move(local_backward_legs)); +} + +PetscForwardOperator PetscOperators::forward() const { + return PetscForwardOperator( + forward_leg_mapping, cell_mapping, meshGetArray("forward_rows"), + meshGetArray("forward_columns"), meshGetArray("forward_weights")); +} + +PetscBackwardOperator PetscOperators::backward() const { + return PetscBackwardOperator( + backward_leg_mapping, cell_mapping, meshGetArray("backward_rows"), + meshGetArray("backward_columns"), meshGetArray("backward_weights")); +} + +template +PetscOperator +PetscOperators::buildInjection(const Field3D& interior_leg_number, + const Field3D& boundary_leg_number, + std::shared_ptr leg_mapping) const { + using Op = PetscOperator; + bout::petsc::UniqueMat mat{new Mat{nullptr}}; + BOUT_DO_PETSC(MatCreate(BoutComm::get(), mat.get())); + BOUT_DO_PETSC(MatSetType(*mat, MATMPIAIJ)); + BOUT_DO_PETSC(MatSetSizes(*mat, leg_mapping->localSize(), cell_mapping->localSize(), + leg_mapping->globalSize(), cell_mapping->globalSize())); + BOUT_DO_PETSC(MatMPIAIJSetPreallocation(*mat, 1, nullptr, 1, nullptr)); + const PetscScalar one = 1.0; + cell_mapping->mapOwnedInteriorCells([&](PetscInt cell_row, const Ind3D& i, int) { + const int interior = ROUND(interior_leg_number[i]); + const int boundary = ROUND(boundary_leg_number[i]); + if (interior >= 0) { + const PetscInt leg_row = leg_mapping->storedToPetsc(interior); + BOUT_DO_PETSC(MatSetValues(*mat, 1, &leg_row, 1, &cell_row, &one, INSERT_VALUES)); + } + if (boundary >= 0) { + const PetscInt leg_row = leg_mapping->storedToPetsc(boundary); + BOUT_DO_PETSC(MatSetValues(*mat, 1, &leg_row, 1, &cell_row, &one, INSERT_VALUES)); + } + }); + BOUT_DO_PETSC(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY)); + BOUT_DO_PETSC(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY)); + return Op(leg_mapping, cell_mapping, std::move(mat)); +} + +template PetscForwardOperator PetscOperators::buildInjection( + const Field3D&, const Field3D&, std::shared_ptr) const; +template PetscBackwardOperator PetscOperators::buildInjection( + const Field3D&, const Field3D&, std::shared_ptr) const; + +PetscForwardOperator PetscOperators::injectForward() const { + return buildInjection( + forward_leg_interior_number, forward_leg_boundary_number, forward_leg_mapping); +} + +PetscBackwardOperator PetscOperators::injectBackward() const { + return buildInjection( + backward_leg_interior_number, backward_leg_boundary_number, backward_leg_mapping); +} + +PetscForwardLegVector PetscOperators::forwardLegWeights() const { + return loadLegWeights("forward_leg_weights", + forward_leg_mapping); +} + +PetscBackwardLegVector PetscOperators::backwardLegWeights() const { + return loadLegWeights("backward_leg_weights", + backward_leg_mapping); +} + +PetscCellOperator PetscOperators::diagonal(const Field3D& f) const { + return PetscCellOperator::diagonal(cell_mapping, makePetscCellVector(cell_mapping, f)); +} + +PetscForwardLegOperator PetscOperators::diagonal(const PetscForwardLegVector& v) const { + return PetscForwardLegOperator::diagonal(forward_leg_mapping, v); +} + +PetscBackwardLegOperator PetscOperators::diagonal(const PetscBackwardLegVector& v) const { + return PetscBackwardLegOperator::diagonal(backward_leg_mapping, v); +} + +PetscOperators::Parallel PetscOperators::getParallel() const { + // Primitive operators from the mesh / metadata + auto Forward = this->forward(); // L+ <- C + auto Backward = this->backward(); // L- <- C + auto Inject_plus = this->injectForward(); // L+ <- C + auto Inject_minus = this->injectBackward(); // L- <- C + + // Leg weights + const auto w_plus_vec = forwardLegWeights(); + const auto w_minus_vec = backwardLegWeights(); + const auto W_plus = diagonal(w_plus_vec); // L+ <- L+ + const auto W_minus = diagonal(w_minus_vec); // L- <- L- + + // Weighted adjoints: C <- L+ and C <- L- + auto Restrict_minus = Inject_plus.transpose() * W_plus; + auto Restrict_plus = Inject_minus.transpose() * W_minus; + + auto* coords = mesh->getCoordinates(); + + // Parallel spacing in cell space + Field3D dl = coords->dy * sqrt(coords->g_22); + dl.splitParallelSlices(); + dl.yup() = 0.0; + dl.ydown() = 0.0; + dl.applyParallelBoundary("parallel_neumann_o1"); + + // Cell volume + Field3D dV = coords->J * coords->dx * coords->dy * coords->dz; + dV.splitParallelSlices(); + dV.yup() = 0.0; + dV.ydown() = 0.0; + dV.applyParallelBoundary("parallel_neumann_o1"); + + Field3D neg_inv_dV = -1.0 / dV; + neg_inv_dV.splitParallelSlices(); + neg_inv_dV.yup() = 0.0; + neg_inv_dV.ydown() = 0.0; + neg_inv_dV.applyParallelBoundary("parallel_neumann_o1"); + + const auto DV = diagonal(dV); + const auto Neg_inv_dV = diagonal(neg_inv_dV); + + // Unweighted interpolation to +1/2 and -1/2 legs + auto Interp_plus = (Forward + Inject_plus) * 0.5; // L+ <- C + auto Interp_minus = (Backward + Inject_minus) * 0.5; // L- <- C + + // Leg-centered dl + const auto dl_plus = Interp_plus(dl); + const auto dl_minus = Interp_minus(dl); + + const auto inv_dl_plus = PetscForwardLegVector::reciprocal(dl_plus); + const auto inv_dl_minus = PetscBackwardLegVector::reciprocal(dl_minus); + + const auto Inv_dl_plus = diagonal(inv_dl_plus); // L+ <- L+ + const auto Inv_dl_minus = diagonal(inv_dl_minus); // L- <- L- + + // Half-step gradients + auto Grad_plus = Inv_dl_plus * (Forward - Inject_plus); // L+ <- C + auto Grad_minus = Inv_dl_minus * (Inject_minus - Backward); // L- <- C + + // Cell-centered central gradient + auto Grad_par = + ((Restrict_minus * Grad_plus) + (Restrict_plus * Grad_minus)) * 0.5; // C <- C + + // Leg-centered volumes from unweighted interpolation + const auto dV_plus = Interp_plus(dV); // L+ + const auto dV_minus = Interp_minus(dV); // L- + + // Physical leg volumes include the leg weights + const auto dV_leg_plus = PetscForwardLegVector::pointwiseMultiply(dV_plus, w_plus_vec); + const auto dV_leg_minus = + PetscBackwardLegVector::pointwiseMultiply(dV_minus, w_minus_vec); + + const auto DV_leg_plus = diagonal(dV_leg_plus); // L+ <- L+ + const auto DV_leg_minus = diagonal(dV_leg_minus); // L- <- L- + + // Support-operator divergence from half-step gradients + auto Div_minus = Neg_inv_dV * Grad_plus.transpose() * DV_leg_plus; // C <- L+ + auto Div_plus = Neg_inv_dV * Grad_minus.transpose() * DV_leg_minus; // C <- L- + + // Cell-centered divergence paired with the cell-centered gradient + auto Div_par = Neg_inv_dV * Grad_par.transpose() * DV; // C <- C + + // Diffusion operator assembled from the half-step pieces + auto Div_par_Grad_par = + ((Div_minus * Grad_plus) + (Div_plus * Grad_minus)) * 0.5; // C <- C + + return Parallel{// Cell <- Cell operators + std::move(Grad_par), std::move(Div_par), std::move(Div_par_Grad_par), + std::move(dV), + + // Leg <-> Cell operators + std::move(Grad_minus), std::move(Grad_plus), std::move(Div_minus), + std::move(Div_plus), + + std::move(Backward), std::move(Forward), std::move(Inject_minus), + std::move(Inject_plus), + + std::move(Interp_minus), std::move(Interp_plus), + + std::move(Restrict_minus), std::move(Restrict_plus)}; +} + +Field3D PetscOperators::Parallel::Div_par_K_Grad_par(const Field3D& K, + const Field3D& f) const { + const auto grad_plus = Grad_plus(f); + const auto grad_minus = Grad_minus(f); + const auto K_plus = Interp_plus(K); + const auto K_minus = Interp_minus(K); + + const auto flux_plus = PetscForwardLegVector::pointwiseMultiply(K_plus, grad_plus); + const auto flux_minus = PetscBackwardLegVector::pointwiseMultiply(K_minus, grad_minus); + + const auto div_minus = Div_minus(flux_plus); + const auto div_plus = Div_plus(flux_minus); + return (toField3D(div_minus, K) + toField3D(div_plus, K)) * 0.5; +} + +#endif diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx index 4ccf27981b..56d21c91d0 100644 --- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx +++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx @@ -4,11 +4,15 @@ #include "imex-bdf2.hxx" +#include #include +#include #include #include #include +#include #include +#include #include #include @@ -119,12 +123,9 @@ static PetscErrorCode FormFunctionForColoring(void* UNUSED(snes), Vec x, Vec f, } static PetscErrorCode imexbdf2PCapply(PC pc, Vec x, Vec y) { - int ierr; - // Get the context IMEXBDF2* s; - ierr = PCShellGetContext(pc, reinterpret_cast(&s)); - CHKERRQ(ierr); + PetscCall(PCShellGetContext(pc, reinterpret_cast(&s))); PetscFunctionReturn(s->precon(x, y)); } @@ -197,8 +198,8 @@ int IMEXBDF2::init() { // Allocate memory and initialise structures u.reallocate(nlocal); for (int i = 0; i < maxOrder; i++) { - uV.emplace_back(Array{nlocal}); - fV.emplace_back(Array{nlocal}); + uV.emplace_back(Array(nlocal)); + fV.emplace_back(Array(nlocal)); timesteps.push_back(timestep); uFac.push_back(0.0); fFac.push_back(0.0); @@ -220,16 +221,12 @@ int IMEXBDF2::init() { } // Initialise PETSc components - int ierr; // Vectors - ierr = VecCreate(BoutComm::get(), &snes_x); - CHKERRQ(ierr); - ierr = VecSetSizes(snes_x, nlocal, PETSC_DECIDE); - CHKERRQ(ierr); - ierr = VecSetFromOptions(snes_x); - CHKERRQ(ierr); - VecDuplicate(snes_x, &snes_f); + PetscCall(VecCreate(BoutComm::get(), &snes_x)); + PetscCall(VecSetSizes(snes_x, nlocal, PETSC_DECIDE)); + PetscCall(VecSetFromOptions(snes_x)); + PetscCall(VecDuplicate(snes_x, &snes_f)); // The SNES solver object(s) constructSNES(&snes); diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index b347de7354..d3d690417b 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -97,6 +97,7 @@ set(serial_tests_source ./mesh/test_invert3x3.cxx ./mesh/test_mesh.cxx ./mesh/test_paralleltransform.cxx + ./mesh/test_petsc_operators.cxx ./solver/test_fakesolver.cxx ./solver/test_fakesolver.hxx ./solver/test_solver.cxx diff --git a/tests/unit/include/bout/test_array.cxx b/tests/unit/include/bout/test_array.cxx index 21be6fe8a5..4efd968f57 100644 --- a/tests/unit/include/bout/test_array.cxx +++ b/tests/unit/include/bout/test_array.cxx @@ -1,9 +1,7 @@ #include "gtest/gtest.h" #include "bout/array.hxx" -#include "bout/boutexception.hxx" -#include #include // In order to keep these tests independent, they need to use @@ -19,7 +17,7 @@ class ArrayTest : public ::testing::Test { }; TEST_F(ArrayTest, ArraySize) { - Array a(5); + const Array a(5); ASSERT_FALSE(a.empty()); EXPECT_EQ(a.size(), 5); @@ -36,7 +34,7 @@ TEST_F(ArrayTest, ArrayValues) { } TEST_F(ArrayTest, CopyArrayConstructor) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -53,7 +51,7 @@ TEST_F(ArrayTest, CopyArrayConstructor) { } TEST_F(ArrayTest, CopyArrayOperator) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -71,7 +69,7 @@ TEST_F(ArrayTest, CopyArrayOperator) { } TEST_F(ArrayTest, CopyArrayNonMemberFunction) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -86,11 +84,11 @@ TEST_F(ArrayTest, CopyArrayNonMemberFunction) { } TEST_F(ArrayTest, SwapArray) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); - Array b{10}; + Array b(10); std::iota(b.begin(), b.end(), 5); @@ -108,7 +106,7 @@ TEST_F(ArrayTest, SwapArray) { } TEST_F(ArrayTest, MoveArrayConstructor) { - Array a{15}; + Array a(15); std::iota(a.begin(), a.end(), 0); @@ -123,7 +121,7 @@ TEST_F(ArrayTest, MoveArrayConstructor) { } TEST_F(ArrayTest, Reallocate) { - Array a{}; + Array a; ASSERT_TRUE(a.empty()); @@ -235,7 +233,7 @@ TEST_F(ArrayTest, RetrieveData) { TEST_F(ArrayTest, Assignment) { Array a(35); - Array b(35); + const Array b(35); // Assign a = b; @@ -251,3 +249,12 @@ TEST_F(ArrayTest, OutOfBoundsThrow) { EXPECT_THROW(a[-1] = 1.0, BoutException); } #endif + +TEST_F(ArrayTest, FromValues) { + const auto a = Array::fromValues({3, 1, 4, 2}); + ASSERT_EQ(4, a.size()); + EXPECT_EQ(3, a[0]); + EXPECT_EQ(1, a[1]); + EXPECT_EQ(4, a[2]); + EXPECT_EQ(2, a[3]); +} diff --git a/tests/unit/include/test_cyclic_reduction.cxx b/tests/unit/include/test_cyclic_reduction.cxx index e057c41bef..5eb9acff0a 100644 --- a/tests/unit/include/test_cyclic_reduction.cxx +++ b/tests/unit/include/test_cyclic_reduction.cxx @@ -2,6 +2,7 @@ #include "test_extras.hxx" #include "bout/array.hxx" +#include "bout/bout_types.hxx" #include "bout/boutcomm.hxx" #include "bout/cyclic_reduction.hxx" @@ -19,7 +20,7 @@ Array makeArrayFromVector(const std::vector& values) { using std::begin; using std::end; - Array array{static_cast::size_type>(values.size())}; + Array array(static_cast::size_type>(values.size())); std::copy(begin(values), end(values), begin(array)); return array; } @@ -28,9 +29,9 @@ Matrix makeMatrixFromVector(const std::vector>& using std::begin; using std::end; - Matrix matrix{static_cast::size_type>(values.size()), - static_cast::size_type>(values[0].size())}; - auto start = begin(matrix); + Matrix matrix(static_cast::size_type>(values.size()), + static_cast::size_type>(values[0].size())); + auto* start = begin(matrix); for (const auto& sub_values : values) { start = std::copy(begin(sub_values), end(sub_values), start); } @@ -48,7 +49,7 @@ TEST(CyclicReduction, SerialSolveSingleArray) { reduce.setCoefs(a, b, c); auto rhs = makeArrayFromVector({0., 1., 2., 2., 3.}); - Array x{reduction_size}; + Array x(reduction_size); reduce.solve(rhs, x); @@ -70,7 +71,7 @@ TEST(CyclicReduction, SerialSolveSingleMatrix) { reduce.setCoefs(a, b, c); auto rhs = makeMatrixFromVector({{0., 1., 2., 2., 3.}}); - Matrix x{1, reduction_size}; + Matrix x(1, reduction_size); reduce.solve(rhs, x); diff --git a/tests/unit/invert/test_fft.cxx b/tests/unit/invert/test_fft.cxx index 1fb30c5eb7..36b344dbb1 100644 --- a/tests/unit/invert/test_fft.cxx +++ b/tests/unit/invert/test_fft.cxx @@ -4,12 +4,12 @@ #include "test_extras.hxx" #include "bout/array.hxx" +#include "bout/bout_types.hxx" #include "bout/constants.hxx" #include "bout/dcomplex.hxx" #include "bout/fft.hxx" #include -#include #include #if BOUT_HAS_FFTW @@ -19,7 +19,7 @@ class FFTTest : public ::testing::TestWithParam { : size(GetParam()), nmodes((size / 2) + 1), real_signal(size), fft_signal(nmodes) { // Make grid indices from [0, size - 1] - Array indices{size}; + Array indices(size); std::iota(indices.begin(), indices.end(), 0.0); // Calculate sin(x) + cos(2x) on [0, 2pi] @@ -48,7 +48,7 @@ INSTANTIATE_TEST_SUITE_P(FFTEvenAndOddSamples, FFTTest, ::testing::Values(8, 9)) TEST_P(FFTTest, rfft) { - Array output{nmodes}; + Array output(nmodes); // Compute forward real FFT rfft(real_signal.begin(), size, output.begin()); @@ -63,7 +63,7 @@ TEST_P(FFTTest, rfft) { TEST_P(FFTTest, irfft) { - Array output{size}; + Array output(size); // Compute inverse real FFT irfft(fft_signal.begin(), size, output.begin()); diff --git a/tests/unit/mesh/test_petsc_operators.cxx b/tests/unit/mesh/test_petsc_operators.cxx new file mode 100644 index 0000000000..2ba63f4528 --- /dev/null +++ b/tests/unit/mesh/test_petsc_operators.cxx @@ -0,0 +1,155 @@ +#include "bout/build_defines.hxx" + +#if BOUT_HAS_PETSC + +#include "gtest/gtest.h" + +#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" + +#include "fake_mesh_fixture.hxx" +#include "test_extras.hxx" + +#include + +// Reuse the "standard" fixture for FakeMesh +using PetscMappingTest = FakeMeshFixture; + +TEST_F(PetscMappingTest, create_region_empty) { + const Field3D cell_number{-1.0}; // No cells >= 0 + + auto rgn = PetscCellMapping::create_region(cell_number); + + ASSERT_EQ(0, rgn.size()); +} + +TEST_F(PetscMappingTest, create_region) { + Field3D cell_number{-1.0}; // No cells >= 0 + + // Note: 1 boundary cell in X and Y + cell_number(0, 0, 0) = 0; // Corner + cell_number(1, 1, 0) = 1; // In domain + cell_number(0, 1, 1) = 2; // Xin boundary + + auto rgn = PetscCellMapping::create_region(cell_number); + ASSERT_EQ(1, rgn.size()); + + const Ind3D expected_ind = cell_number.indexAt(1, 1, 0); + ASSERT_TRUE(expected_ind == *rgn.begin()); +} + +TEST_F(PetscMappingTest, create_region_xin) { + Field3D cell_number{-1.0}; // No cells >= 0 + + // Note: 1 boundary cell in X and Y + cell_number(0, 0, 0) = 0; // Corner + cell_number(1, 1, 0) = 1; // In domain + cell_number(0, 1, 1) = 2; // Xin boundary + + auto rgn = PetscCellMapping::create_region_xin(cell_number); + ASSERT_EQ(1, rgn.size()); + + const Ind3D expected_ind = cell_number.indexAt(0, 1, 1); + output.write("Expecting {} == {}\n", expected_ind, *rgn.begin()); + ASSERT_TRUE(expected_ind == *rgn.begin()); +} + +TEST_F(PetscMappingTest, mapping) { + Field3D cell_number{-1.0}; // No cells >= 0 + + // Note: 1 boundary cell in X and Y + cell_number(0, 0, 0) = 2; // Corner + cell_number(1, 1, 0) = 0; // In domain + cell_number(0, 1, 1) = 1; // Xin boundary + + const Field3D forward_cell_number{-1.0}; + const Field3D backward_cell_number{-1.0}; + + const PetscCellMapping mapping(cell_number, forward_cell_number, backward_cell_number, + 2); + + // Two cells: one evolving and one in xin + ASSERT_EQ(2, mapping.localSize()); + ASSERT_EQ(2, mapping.globalSize()); +} + +using PetscOperatorTest = FakeMeshFixture; + +TEST_F(PetscOperatorTest, identity) { + Field3D cell_number{-1.0}; + const Field3D forward_cell_number{-1.0}; + const Field3D backward_cell_number{-1.0}; + + // Three cells active + cell_number(1, 1, 0) = 0; + cell_number(1, 2, 0) = 1; + cell_number(1, 2, 1) = 2; + + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, 3); + ASSERT_EQ(3, mapping->localSize()); + + const auto rows = Array::fromValues({0, 1, 2, 3}); + const auto cols = Array::fromValues({0, 1, 2}); + const auto weights = Array::fromValues({1.0, 1.0, 1.0}); + + const PetscCellOperator identity(mapping, mapping, rows, cols, weights); + + Field3D value{0.0}; + value(1, 1, 0) = 10; + value(1, 2, 0) = 21; + value(1, 2, 1) = 32; + value.splitParallelSlices(); + value.yup() = -1.0; + value.ydown() = -1.0; + + const Field3D result = identity(value); + + EXPECT_EQ(10, result(1, 1, 0)); + EXPECT_EQ(21, result(1, 2, 0)); + EXPECT_EQ(32, result(1, 2, 1)); +} + +TEST_F(PetscOperatorTest, identity_yupdown) { + Field3D cell_number{-1.0}; + Field3D forward_cell_number{-1.0}; + Field3D backward_cell_number{-1.0}; + + // Three cells active + cell_number(1, 1, 0) = 0; + cell_number(1, 2, 0) = 1; + cell_number(1, 2, 1) = 2; + + forward_cell_number(1, 2, 0) = 3; + backward_cell_number(1, 1, 2) = 4; + + auto mapping = std::make_shared( + cell_number, forward_cell_number, backward_cell_number, 5); + ASSERT_EQ(5, mapping->localSize()); + + const auto rows = Array::fromValues({0, 1, 2, 3}); + const auto cols = Array::fromValues({0, 1, 2}); + const auto weights = Array::fromValues({1.0, 1.0, 1.0}); + + const PetscCellOperator identity(mapping, mapping, rows, cols, weights); + + Field3D value{0.0}; + value(1, 1, 0) = 10; + value(1, 2, 0) = 21; + value(1, 2, 1) = 32; + value.splitParallelSlices(); + value.yup() = -1.0; + value.ydown() = -1.0; + + const Field3D result = identity(value); + + EXPECT_EQ(10, result(1, 1, 0)); + EXPECT_EQ(21, result(1, 2, 0)); + EXPECT_EQ(32, result(1, 2, 1)); +} + +#endif // BOUT_HAS_PETSC