Calculate Jacobian coloring using PetscOperator matrices#3350
Calculate Jacobian coloring using PetscOperator matrices#3350bendudson wants to merge 20 commits into
Conversation
This function should create an index set of the evolving cells. It will be used to extract the part of the PetscOperator that represents coupling between evolving cells. Implementation to follow.
Creates an index set of global PETSc rows that correspond to evolving cells. Passes unit tests.
Method that will extract the part of a PetscCellOperator that couples evolving cells to other evolving cells. Dummy implementation with failing tests.
Passes unit tests
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
| // Collect global PETSc indices in mapOwnedInteriorCells order. | ||
| // Reserve the known count up front to avoid reallocation. | ||
| std::vector<PetscInt> indices; | ||
| indices.reserve(static_cast<std::size_t>(evolving_region.size())); |
There was a problem hiding this comment.
warning: no header providing "std::size_t" is directly included [misc-include-cleaner]
src/mesh/petsc_operators.cxx:1:
+ #include <cstddef>| mapOwnedInteriorCells( | ||
| [&](PetscInt row, const Ind3D& /*i*/, int /*stored*/) { indices.push_back(row); }); | ||
|
|
||
| IS is; |
There was a problem hiding this comment.
warning: variable 'is' is not initialized [cppcoreguidelines-init-variables]
| IS is; | |
| IS is = nullptr; |
| const PetscOperator<CellSpaceTag, CellSpaceTag>& op) const { | ||
| IS is = makeEvolvingIS(); | ||
|
|
||
| Mat sub; |
There was a problem hiding this comment.
warning: variable 'sub' is not initialized [cppcoreguidelines-init-variables]
| Mat sub; | |
| Mat sub = nullptr; |
Checks the global ordering of cells for consistency as number of processors and nxpe is varied.
Include guards, header includes.
This function will add non-zeros to a Jacobian matrix, based on a submatrix from a PetscCellOperator. Dummy implementation so unit tests fail.
Sets a sparsity pattern in a given Jacobian matrix, based on an operator matrix and the variable indices. Passes unit tests.
Should apply the connectivity pattern to all pairs of variables in the system Jacobian.
Applies the same connectivity in `mat` to all pairs of variables in `J`. Passes unit tests.
|
|
||
| void addOperatorSparsity(Mat Jfd, Mat sub, int out_var, int in_var) { | ||
| // Infer nvars from global sizes | ||
| PetscInt jfd_global{0}, sub_global{0}; |
There was a problem hiding this comment.
warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]
| PetscInt jfd_global{0}, sub_global{0}; | |
| PetscInt jfd_global{0}; | |
| PetscInt sub_global{0}; |
| ASSERT1(in_var >= 0 && in_var < static_cast<int>(nvars)); | ||
|
|
||
| // Iterate over locally owned rows of sub and insert into Jfd | ||
| PetscInt rstart{0}, rend{0}; |
There was a problem hiding this comment.
warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]
| PetscInt rstart{0}, rend{0}; | |
| PetscInt rstart{0}; | |
| PetscInt rend{0}; |
| } | ||
|
|
||
| void addOperatorSparsity(Mat Jfd, Mat sub) { | ||
| PetscInt jfd_global{0}, sub_global{0}; |
There was a problem hiding this comment.
warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]
| PetscInt jfd_global{0}, sub_global{0}; | |
| PetscInt jfd_global{0}; | |
| PetscInt sub_global{0}; |
| #include <bout/petsc_operators.hxx> | ||
| #include <bout/region.hxx> | ||
| #include <bout/utils.hxx> | ||
|
|
There was a problem hiding this comment.
warning: included header utils.hxx is not used directly [misc-include-cleaner]
| MatDestroy(&dummy); | ||
|
|
||
| // ── Compare traversal orderings ─────────────────────────────────────────── | ||
| // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that |
There was a problem hiding this comment.
warning: variable 'Istart' is not initialized [cppcoreguidelines-init-variables]
RMINE);
^this fix will not be applied because it overlaps with another fix
|
|
||
| // ── Compare traversal orderings ─────────────────────────────────────────── | ||
| // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that | ||
| // the local offset produced by each is identical for every cell. |
There was a problem hiding this comment.
warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]
| // the local offset produced by each is identical for every cell. | |
| RMINE); | |
| t, Iend;PetscInt Istart; | |
| PetscInt Iend; |
|
|
||
| // ── Compare traversal orderings ─────────────────────────────────────────── | ||
| // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that | ||
| // the local offset produced by each is identical for every cell. |
There was a problem hiding this comment.
warning: variable 'Iend' is not initialized [cppcoreguidelines-init-variables]
t, Iend;
^this fix will not be applied because it overlaps with another fix
| } | ||
|
|
||
| if (my_rank == 0) { | ||
| output.write("total_mismatches={}\n", total_mismatches); |
There was a problem hiding this comment.
warning: variable 'my_rank' is not initialized [cppcoreguidelines-init-variables]
tests/integrated/test-petsc-ordering/test_petsc_ordering.cxx:201:
- k;
+ k = 0;Test passes locally but fails in CI
Use global numbering
| } | ||
| } | ||
| } | ||
| // Communicate so guard cells are filled (needed by some mesh operations) |
There was a problem hiding this comment.
warning: '*' has higher precedence than '+'; add parentheses to explicitly specify the order of operations [readability-math-missing-parentheses]
| // Communicate so guard cells are filled (needed by some mesh operations) | |
| ell_number(x, y, z) = (((global_x * ngy) + global_y) * ngz) + global_z; |
| MatDestroy(&dummy); | ||
|
|
||
| // ── Compare traversal orderings ─────────────────────────────────────────── | ||
| // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that |
There was a problem hiding this comment.
warning: variable 'Iend' is not initialized [cppcoreguidelines-init-variables]
MatSetType(dummy, MATMPIAIJ);
^this fix will not be applied because it overlaps with another fix
| MatDestroy(&dummy); | ||
|
|
||
| // ── Compare traversal orderings ─────────────────────────────────────────── | ||
| // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that |
There was a problem hiding this comment.
warning: variable 'Istart' is not initialized [cppcoreguidelines-init-variables]
MatSetType(dummy, MATMPIAIJ);
^this fix will not be applied because it overlaps with another fix
|
|
||
| // ── Compare traversal orderings ─────────────────────────────────────────── | ||
| // Walk mapOwnedInteriorCells and RGN_NOBNDRY in parallel, checking that | ||
| // the local offset produced by each is identical for every cell. |
There was a problem hiding this comment.
warning: multiple declarations in a single statement reduces readability [readability-isolate-declaration]
| // the local offset produced by each is identical for every cell. | |
| MatSetType(dummy, MATMPIAIJ);PetscInt Istart; | |
| PetscInt Iend; |
dschwoerer
force-pushed
the
petsc-operators
branch
from
562e38c to
8019a94
Compare
| // Collect global PETSc indices in mapOwnedInteriorCells order. | ||
| // Reserve the known count up front to avoid reallocation. | ||
| std::vector<PetscInt> indices; | ||
| indices.reserve(static_cast<std::size_t>(evolving_region.size())); |
There was a problem hiding this comment.
warning: no header providing "std::size_t" is directly included [misc-include-cleaner]
src/mesh/petsc_operators.cxx:5:
+ #include <cstddef>Field2D expressions need to be evaluated before assignment to Field3D.
Handle the case where `source` is `nullptr`.
Assumes no evolving Field2D or evolving boundary cells.
Allows the sparsity pattern in PetscPreconditioner for SNES, CVODE and PETSc interfaces to be set using a PetscOperator. This should enable preconditioning of 3D FCI simulations.
Some reorganisation, removing duplicate and obsolete sections.
| #if BOUT_HAS_PETSC | ||
| struct DeferredJacobianPattern { | ||
| bout::petsc::UniqueMat submatrix{new Mat{nullptr}}; | ||
| VarRef out_var{}; |
There was a problem hiding this comment.
warning: initializer for member 'out_var' is redundant [readability-redundant-member-init]
| VarRef out_var{}; | |
| VarRef out_var; |
| struct DeferredJacobianPattern { | ||
| bout::petsc::UniqueMat submatrix{new Mat{nullptr}}; | ||
| VarRef out_var{}; | ||
| VarRef in_var{}; |
There was a problem hiding this comment.
warning: initializer for member 'in_var' is redundant [readability-redundant-member-init]
| VarRef in_var{}; | |
| VarRef in_var; |
| int run() override; | ||
| BoutReal run(BoutReal tout); | ||
| #if BOUT_HAS_PETSC | ||
| bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, |
There was a problem hiding this comment.
warning: no header providing "PetscCellOperator" is directly included [misc-include-cleaner]
src/solver/impls/cvode/cvode.hxx:32:
- #include "bout/solver.hxx"
+ #include "bout/petsc_operators.hxx"
+ #include "bout/solver.hxx"|
|
||
| int init() override; | ||
| int run() override; | ||
| bool addJacobianPattern(const PetscCellOperator& op, VarRef out_var, |
There was a problem hiding this comment.
warning: no header providing "PetscCellOperator" is directly included [misc-include-cleaner]
src/solver/impls/petsc/petsc.hxx:29:
- #include "bout/solver.hxx"
+ #include "bout/petsc_operators.hxx"
+ #include "bout/solver.hxx"|
|
||
| int init() override; | ||
| int run() override; | ||
| #if BOUT_HAS_PETSC |
There was a problem hiding this comment.
warning: nested redundant #if; consider removing it [readability-redundant-preprocessor]
#if BOUT_HAS_PETSC
^Additional context
src/solver/impls/snes/snes.hxx:33: previous #if was here
#if BOUT_HAS_PETSC
^| } | ||
|
|
||
| #if BOUT_HAS_PETSC | ||
| bool Solver::addJacobianPattern(const PetscCellOperator& op) { |
There was a problem hiding this comment.
warning: no header providing "PetscCellOperator" is directly included [misc-include-cleaner]
src/solver/solver.cxx:39:
- #include "bout/region.hxx"
+ #include "bout/petsc_operators.hxx"
+ #include "bout/region.hxx"| } | ||
|
|
||
| auto out_mapping = | ||
| std::dynamic_pointer_cast<const PetscCellMapping>(op.getOutMapping()); |
There was a problem hiding this comment.
warning: no header providing "PetscCellMapping" is directly included [misc-include-cleaner]
std::dynamic_pointer_cast<const PetscCellMapping>(op.getOutMapping());
^| } | ||
|
|
||
| const auto evolve_boundary = [](const auto& field) { return field.evolve_bndry; }; | ||
| return std::none_of(f2d.begin(), f2d.end(), evolve_boundary) |
There was a problem hiding this comment.
warning: no header providing "std::none_of" is directly included [misc-include-cleaner]
src/solver/solver.cxx:48:
+ #include <algorithm>|
|
||
| const auto evolve_boundary = [](const auto& field) { return field.evolve_bndry; }; | ||
| return std::none_of(f2d.begin(), f2d.end(), evolve_boundary) | ||
| && std::none_of(f3d.begin(), f3d.end(), evolve_boundary); |
There was a problem hiding this comment.
warning: no header providing "std::none_of" is directly included [misc-include-cleaner]
&& std::none_of(f3d.begin(), f3d.end(), evolve_boundary);
^| const int nvars = n2Dvars() + n3Dvars(); | ||
| ASSERT1(nvars > 0); | ||
|
|
||
| BOUT_DO_PETSC(MatSetOption(Jfd, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); |
There was a problem hiding this comment.
warning: no header providing "BOUT_DO_PETSC" is directly included [misc-include-cleaner]
src/solver/solver.cxx:39:
- #include "bout/region.hxx"
+ #include "bout/petsclib.hxx"
+ #include "bout/region.hxx"
bendudson
changed the title
1 participant
Builds on #3330 , using
PetscOperatormatrices to infer connectivity between cells.The aim is to accelerate transport simulations with FCI by inferring the sparsity pattern of the Jacobian in
SNESSolverandPetscSolver. This will enable coloring to efficiently calculate finite difference Jacobian approximations for LU-like preconditioners. This is the method used for structured grid (e.g. Tokamak) simulations, that so far has not been possible with FCI.