Petsc matrix operators for FCI#3330
Conversation
Read operators from the mesh in CSR format. PETSc handles communication pattern for MPI decomposition.
Handles the mapping between mesh indices and PETSc global indices.
Simple wrapper functions around the Options functionality, so that Array types can be read from the mesh file.
Reads row, column and weights from Mesh, constructs Petsc Mat objects. Working on composition and operation on Field3D objects.
Map from Field3D, including yup/ydown boundaries, to PETSc Vec, perform matrix-vector product, and map back to a Field3D.
Docstrings to explain PetscMapping, PetscOperator and PetscOperators classes. Extend `BOUT_DO_PETSC` macro to print error message, file and line number.
Alternative to initializer_list constructor that avoids
confusion with e.g. `Array{2}` : Is it an array of size 2,
or an array of size 1 containing value 2?
`std::vector` uses `vector(2)` means size 2; `vector{2}`
means value 2.
With this commit, both `Array{2}` and `Array(2)` mean size 2;
`Array<T>::fromValues({2})` means size 1, value 2. It's more verbose
but less prone to typos.
Some care needed with yup/down and yp/ym indices.
`diagonal` creates diagonal matrices from Field3D. `transpose` creates a new matrix that is a transpose of a PetscOperator. This is useful for converting gradients into divergences.
Allows input operators to contain entries in boundary rows. This is important so that the gradient and divergence operators depend on boundary cell values.
![github-actions[bot]](https://avatars.github.com/in/15368?s=60&v=4){kind=small}
| /// The local index is the offset into the packed local PETSc vector. | ||
| template <typename Function> | ||
| void map_local_yup(Function func) const { | ||
| PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size(); |
There was a problem hiding this comment.
warning: narrowing conversion from 'unsigned int' to signed type 'PetscInt' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]
PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size();
^| /// The local index is the offset into the packed local PETSc vector. | ||
| template <typename Function> | ||
| void map_local_ydown(Function func) const { | ||
| PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size() |
There was a problem hiding this comment.
warning: narrowing conversion from 'unsigned int' to signed type 'PetscInt' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]
PetscInt row = evolving_region.size() + xin_region.size() + xout_region.size()
^| /// - applied to a `Field3D`; | ||
| /// - composed with other operators; and | ||
| /// - added or subtracted. | ||
| class PetscOperator { |
There was a problem hiding this comment.
warning: class 'PetscOperator' defines a destructor but does not define a copy constructor, a copy assignment operator, a move constructor or a move assignment operator [cppcoreguidelines-special-member-functions]
class PetscOperator {
^| int mesh_total_cells; | ||
| if (mesh->get(mesh_total_cells, "total_cells") == 0) { | ||
| // Check total number of cells | ||
| if (mesh_total_cells != mapping->size()) { |
There was a problem hiding this comment.
warning: comparison of integers of different signs: 'int' and 'unsigned int' [clang-diagnostic-sign-compare]
if (mesh_total_cells != mapping->size()) {
^| IMEXBDF2* s; | ||
| ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s)); | ||
| CHKERRQ(ierr); | ||
| PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s))); |
There was a problem hiding this comment.
warning: do not use reinterpret_cast [cppcoreguidelines-pro-type-reinterpret-cast]
PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));
^| IMEXBDF2* s; | ||
| ierr = PCShellGetContext(pc, reinterpret_cast<void**>(&s)); | ||
| CHKERRQ(ierr); | ||
| PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s))); |
There was a problem hiding this comment.
warning: multilevel pointer conversion from 'void **' to 'void *', please use explicit cast [bugprone-multi-level-implicit-pointer-conversion]
PetscCall(PCShellGetContext(pc, reinterpret_cast<void**>(&s)));
^Not available unless configured with PETSc.
| return Region<Ind3D>(xout_indices); | ||
| } | ||
|
|
||
| PetscMapping::PetscMapping(const Field3D& cell_number, const Field3D& forward_cell_number, |
There was a problem hiding this comment.
warning: constructor does not initialize these fields: row_start, row_end, mat_mesh_to_petsc, mat_petsc_to_mesh [cppcoreguidelines-pro-type-member-init]
include/bout/petsc_operators.hxx:236:
- PetscInt row_start, row_end;
+ PetscInt row_start{}, row_end{};include/bout/petsc_operators.hxx:239:
- Mat mat_mesh_to_petsc;
+ Mat mat_mesh_to_petsc{};include/bout/petsc_operators.hxx:242:
- Mat mat_petsc_to_mesh;
+ Mat mat_petsc_to_mesh{};| MatDestroy(&this->mat_petsc_to_mesh); | ||
| } | ||
|
|
||
| PetscOperator::PetscOperator(PetscMappingPtr mapping, Array<int> rows, Array<int> cols, |
There was a problem hiding this comment.
warning: constructor does not initialize these fields: mat_operator [cppcoreguidelines-pro-type-member-init]
include/bout/petsc_operators.hxx:335:
- Mat mat_operator;
+ Mat mat_operator{};| BOUT_DO_PETSC(MatDestroy(&mat)); | ||
| } | ||
|
|
||
| void PetscOperator::copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec) { |
There was a problem hiding this comment.
warning: the parameter 'mapping' is copied for each invocation but only used as a const reference; consider making it a const reference [performance-unnecessary-value-param]
| void PetscOperator::copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec) { | |
| void PetscOperator::copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec) { |
include/bout/petsc_operators.hxx:345:
- static void copyToVec(PetscMappingPtr mapping, const Field3D& f, Vec vec);
+ static void copyToVec(const PetscMappingPtr& mapping, const Field3D& f, Vec vec);| return PetscOperator(this->mapping, mat); | ||
| } | ||
|
|
||
| Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const { |
There was a problem hiding this comment.
warning: method 'meshGetField3D' can be made static [readability-convert-member-functions-to-static]
| Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) const { | |
| Field3D PetscOperators::meshGetField3D(Mesh* mesh, const std::string& name) { |
include/bout/petsc_operators.hxx:392:
- Field3D meshGetField3D(Mesh* mesh, const std::string& name) const;
+ static Field3D meshGetField3D(Mesh* mesh, const std::string& name) ;Check number of cells against global mapping size, not local size.
Handle non-standard CSR inputs: - Negative indices for missing rows - Missing final entry in rows array
Changes PetscMapping::size() to PetscMapping::localSize()
Defined in `bout/petsc_interface.hxx`, `bout::petsc::UniqueVec and `UniqueMat` are `std::unique_ptr` wrappers around PETSc Vec and Mat pointers, with custom deleters. PetscOperator copy constructor and assignment deleted, move assignment set to default. This enables some optimizations to operate in-place on temporaries.
|
This looks amazing, thanks @bendudson! I've not finished going through this yet, but how much is this a replacement for |
| struct VecDeleter { | ||
| void operator()(Vec* vec_ptr) const { | ||
| VecDestroy(vec_ptr); | ||
| delete vec_ptr; |
There was a problem hiding this comment.
warning: deleting a pointer through a type that is not marked 'gsl::owner<>'; consider using a smart pointer instead [cppcoreguidelines-owning-memory]
delete vec_ptr;
^Additional context
include/bout/petsc_interface.hxx:67: variable declared here
void operator()(Vec* vec_ptr) const {
^| struct MatDeleter { | ||
| void operator()(Mat* mat_ptr) const { | ||
| MatDestroy(mat_ptr); | ||
| delete mat_ptr; |
There was a problem hiding this comment.
warning: deleting a pointer through a type that is not marked 'gsl::owner<>'; consider using a smart pointer instead [cppcoreguidelines-owning-memory]
delete mat_ptr;
^Additional context
include/bout/petsc_interface.hxx:80: variable declared here
void operator()(Mat* mat_ptr) const {
^Shows how to use the PetscOperators, and compares against the yup/ydown operators. Includes a script to generate the mesh, which currently requires a branch of Zoidberg to run. The code to calculate Grad_par, Div_par and Div_par_Grad_par from forward and backward operators is now in `PetscOperators::getParallel()`. That calculates and returns a collection of `Parallel` operators. This example could be turned into a test: - `Grad_par` should be close between the two methods - Volume integrals of `Div_par` and `Div_par_Grad_par` results should be near zero (flux conserved).
|
Thanks for going through this @ZedThree ! Sorry there's quite a lot, though most is self-contained. The implementation is independent of I've added an example of how these operators are used (https://github.com/boutproject/BOUT-dev/pull/3330/changes#diff-d24773f587abce8257e893cfcfec54002aafaf875fb9cf93f54076b5c60ec7abR16). There is some amount of setup needed to calculate differential operators from the forward / backward maps, that I've put into a helper function that returns a struct of operators (https://github.com/boutproject/BOUT-dev/pull/3330/changes#diff-62e947100b5e999ca6d6d57a6f1281449adc7e2db9e8b7fc5e9273fea887ad88R302). Not sure this is the right API; suggestions welcome! |
I had thought that's essentially what |
The |
|
Ahhh, I see! Thanks! |
This is a function rather than a PetscOperator because it requires at least two matrix-vector products. This implementation requires six Mat-Vec products that could probably be optimized.
Trace multiple field lines per cell, to construct smoother maps. A complication is cells that partially intersect a boundary.
Mesh not available in repo, only generation script.
bendudson
changed the title
Since correctness depends on a combination of mesh generation and BOUT++ implementation, this integrated test installs Zoidberg in a virtual environment to generate a mesh. Uses pytest for forward compatibility
Test passes locally, fails in CI.
One error occurs when allocating:
```
Option total_cells = 4828 (grid.fci.nc)
Option backward_cell_number = Tensor<BoutReal> (grid.fci.nc)
Option backward_cell_number = Tensor<BoutReal> (grid.fci.nc)
Option forward_cell_number = Tensor<BoutReal> (grid.fci.nc)
Option forward_cell_number = Tensor<BoutReal> (grid.fci.nc)
Option cell_number = Tensor<BoutReal> (grid.fci.nc)
Option cell_number = Tensor<BoutReal> (grid.fci.nc)
STDERR:
terminate called after throwing an instance of 'std::bad_array_new_length'
what(): std::bad_array_new_length
[runnervm3jyl0:12838] *** Process received signal ***
[runnervm3jyl0:12838] Signal: Aborted (6)
[runnervm3jyl0:12838] Signal code: (-6)
malloc(): invalid size (unsorted)
```
A second error is due to missing `netCDF4` library
```
FAILED test_petsc_operators.py::test_petsc_operators - ImportError:
DataFile: No supported NetCDF modules available -- requires netCDF4
```
but netCDF4 seems to be installed:
```
Successfully installed boutdata-0.4.0 certifi-2026.5.20 cftime-1.6.5 contourpy-1.3.3 cycler-0.12.1 fonttools-4.63.0 fortranformat-2.0.3 freeqdsk-0.5.2 kiwisolver-1.5.0 matplotlib-3.10.9 mpmath-1.3.0 natsort-8.4.0 netCDF4-1.7.4 numpy-2.4.6 packaging-26.2 pillow-12.2.0 pyparsing-3.3.2 python-dateutil-2.9.0.post0 scipy-1.17.1 shapely-2.1.2 six-1.17.0 sympy-1.14.0 zoidberg-0.3.1.dev3+gfc5e428
```
Attempting to resolve netCDF4 issue
Trying to get stack trace on allocation failure.
Enclose PetscOperators in a scope that ends before BoutFinalise is called, so that the last PetscLib destructor is called before MPI_Finalize.
dschwoerer
force-pushed
the
petsc-operators
branch
from
45bdc3f to
b47b904
Compare
|
I noticed we run 56 test in the "static, optimized case", including the test-petsc-operator, while fedora runs 52 + 8 + 1 - and the The Fedora CI run has 3 additional tests disabled, according to the configure stage: So on the non-fedora CI should have 64 tests run - so 8 tests are missing. Any idea why some tests are not run @ZedThree ? |
| Ind3D _ind() const { return region->bndry_points[pos].index; } | ||
| signed char _boundary_width() const { | ||
| return region->localmesh->ystart - region->bndry_points[pos].abs_offset; | ||
| return region->localmesh->ystart - region->bndry_points[pos].abs_offset + 1; |
There was a problem hiding this comment.
warning: narrowing conversion from 'int' to signed type 'signed char' is implementation-defined [bugprone-narrowing-conversions]
{
^
dschwoerer
force-pushed
the
petsc-operators
branch
2 times, most recently
from
562e38c to
8019a94
Compare
4 participants
Reads matrix operators from the mesh file in CSR format, creating PETSc matrices that apply to Field3Ds. Mainly intended to implement FCI parallel operators as an alternative to yup/down communication.
This enables:
snessolver needs to know which cells are connected, so that PETSc can perform coloring and calculate the system Jacobian efficiently. Knowing the forward and backward matrices will enable this to be done.There are a couple of downsides: