6#ifndef PRECICE_NO_PETSC
17#include "petscdrawtypes.h"
20#include "petscsystypes.h"
22#include "petscviewertypes.h"
32#ifndef PRECICE_NO_PETSC
41#ifndef PRECICE_NO_PETSC
42 PetscBool petscIsInitialized;
43 PetscInitialized(&petscIsInitialized);
44 if (not petscIsInitialized) {
45 PETSC_COMM_WORLD = comm;
47 PetscOptionsSetValue(
nullptr,
"-no_signal_handler",
nullptr);
50 char **argv =
nullptr;
51 ierr = PetscInitialize(&argc, &argv,
"",
nullptr);
60#ifndef PRECICE_NO_PETSC
64 PetscBool petscIsInitialized;
65 PetscInitialized(&petscIsInitialized);
66 if (petscIsInitialized) {
67 PetscOptionsSetValue(
nullptr,
"-options_left",
"no");
74#ifndef PRECICE_NO_PETSC
79#include "petscviewer.h"
86 PetscErrorCode ierr = 0;
87 if (format ==
ASCII) {
88 ierr = PetscViewerASCIIOpen(comm, filename.
c_str(), &
viewer);
90 }
else if (format ==
BINARY) {
91 ierr = PetscViewerBinaryOpen(comm, filename.
c_str(), FILE_MODE_WRITE, &
viewer);
99 PetscErrorCode ierr = 0;
100 ierr = PetscViewerCreate(comm, &
viewer);
102 ierr = PetscViewerSetType(
viewer, PETSCVIEWERASCII);
109 PetscViewerPopFormat(
viewer);
111 PetscViewerDestroy(&
viewer);
116 auto ierr = PetscViewerPushFormat(
viewer, format);
129 PetscObjectGetComm(
reinterpret_cast<PetscObject
>(obj), &comm);
136 PetscErrorCode ierr = 0;
137 ierr = PetscObjectSetName(
reinterpret_cast<PetscObject
>(obj),
name.
c_str());
145 PetscObjectGetName(
reinterpret_cast<PetscObject
>(obj), &cstr);
153 PetscErrorCode ierr = 0;
170 vector = other.vector;
171 other.vector =
nullptr;
184 PetscErrorCode ierr = 0;
198 PetscErrorCode ierr = 0;
199 PetscBool petscIsInitialized;
200 PetscInitialized(&petscIsInitialized);
201 if (petscIsInitialized &&
vector)
202 ierr = VecDestroy(&
vector);
219 PetscErrorCode ierr = 0;
220 ierr = VecDuplicate(other, &newvector);
221 [&] { CHKERRV(ierr); }();
234 PetscErrorCode ierr = 0;
236 ierr = MatCreateVecs(m,
nullptr, &newvector);
238 ierr = MatCreateVecs(m, &newvector,
nullptr);
240 [&] { CHKERRV(ierr); }();
247 swap(vector, other.vector);
250Vector::operator Vec &()
250Vector::operator Vec &() {
…}
257 PetscErrorCode ierr = 0;
258 ierr = VecSetSizes(
vector, PETSC_DECIDE, rows);
260 ierr = VecSetFromOptions(
vector);
266 PetscErrorCode ierr = 0;
268 ierr = VecGetSize(
vector, &size);
275 PetscErrorCode ierr = 0;
277 ierr = VecGetLocalSize(
vector, &size);
284 PetscErrorCode ierr = 0;
285 ierr = VecSetValue(
vector, row, value, INSERT_VALUES);
291 PetscErrorCode ierr = 0;
293 PetscInt range_start, range_end, size;
294 VecGetSize(
vector, &size);
295 VecGetOwnershipRange(
vector, &range_start, &range_end);
296 double step_size = (stop - start) / size;
297 ierr = VecGetArray(
vector, &a);
299 for (PetscInt i = range_start; i < range_end; i++) {
300 a[i - range_start] = (i + start) * step_size;
302 VecRestoreArray(
vector, &a);
307 PetscErrorCode ierr = 0;
314 PetscRandomSetType(rctx, PETSCRAND48);
315 PetscRandomSetSeed(rctx, dist(rd));
316 PetscRandomSeed(rctx);
317 ierr = VecSetRandom(
vector, rctx);
319 PetscRandomDestroy(&rctx);
330 PetscErrorCode ierr = 0;
331 ierr = VecGetArray(
vector, &data);
334 ierr = VecRestoreArray(
vector, &data);
348 PetscErrorCode ierr = 0;
349 ierr = VecGetArray(
vector, &data);
351 auto dataEnd =
std::next(data, localSize);
353 ierr = VecRestoreArray(
vector, &data);
360 PetscErrorCode ierr = 0;
363 ierr = VecGetArray(
vector, &a);
365 ierr = VecGetSize(
vector, &size);
367 ierr = PetscSortReal(size, a);
369 ierr = VecRestoreArray(
vector, &a);
375 PetscErrorCode ierr = 0;
376 ierr = VecAssemblyBegin(
vector);
378 ierr = VecAssemblyEnd(
vector);
384 PetscInt range_start, range_end;
385 VecGetOwnershipRange(
vector, &range_start, &range_end);
392 VecView(
vector, viewer.viewer);
398 VecNorm(
vector, NORM_2, &val);
405 VecLoad(
vector, viewer.viewer);
411 ierr = VecView(
vector, PETSC_VIEWER_STDOUT_WORLD);
424 PetscErrorCode ierr = 0;
432 PetscErrorCode ierr = 0;
433 PetscBool petscIsInitialized;
434 PetscInitialized(&petscIsInitialized);
435 if (petscIsInitialized &&
matrix)
436 ierr = MatDestroy(&
matrix);
440Matrix::operator Mat &()
440Matrix::operator Mat &() {
…}
447 PetscErrorCode ierr = 0;
448 ierr = MatAssemblyBegin(
matrix, type);
450 ierr = MatAssemblyEnd(
matrix, type);
454void Matrix::init(PetscInt localRows, PetscInt localCols, PetscInt globalRows, PetscInt globalCols,
455 MatType type,
bool doSetup)
457 PetscErrorCode ierr = 0;
458 if (type !=
nullptr) {
459 ierr = MatSetType(
matrix, type);
462 ierr = MatSetSizes(
matrix, localRows, localCols, globalRows, globalCols);
464 ierr = MatSetFromOptions(
matrix);
454void Matrix::init(PetscInt localRows, PetscInt localCols, PetscInt globalRows, PetscInt globalCols, {
…}
473 PetscErrorCode ierr = 0;
475 ierr = MatDestroy(&
matrix);
485 MatGetInfo(
matrix, flag, &info);
491 PetscErrorCode ierr = 0;
492 ierr = MatSetValue(
matrix, row, col, value, INSERT_VALUES);
498 PetscErrorCode ierr = 0;
505 PetscRandomSetType(rctx, PETSCRAND48);
506 PetscRandomSetSeed(rctx, dist(rd));
507 PetscRandomSeed(rctx);
508 ierr = MatSetRandom(
matrix, rctx);
510 PetscRandomDestroy(&rctx);
515 PetscErrorCode ierr = 0;
516 const PetscScalar *vec;
517 PetscInt range_start, range_end;
518 VecGetOwnershipRange(v.
vector, &range_start, &range_end);
522 ierr = VecGetArrayRead(v.
vector, &vec);
524 ierr = MatSetValues(
matrix, range_end - range_start, irow.
data(), 1, &col, vec, INSERT_VALUES);
526 ierr = VecRestoreArrayRead(v.
vector, &vec);
528 ierr = MatAssemblyBegin(
matrix, MAT_FINAL_ASSEMBLY);
530 ierr = MatAssemblyEnd(
matrix, MAT_FINAL_ASSEMBLY);
537 MatGetSize(
matrix, &m, &n);
544 MatGetLocalSize(
matrix, &m, &n);
550 PetscInt range_start, range_end;
551 MatGetOwnershipRange(
matrix, &range_start, &range_end);
557 PetscInt range_start, range_end;
558 MatGetOwnershipRangeColumn(
matrix, &range_start, &range_end);
564 PetscErrorCode ierr = 0;
566 ierr = MatGetBlockSize(
matrix, &bs);
573 PetscErrorCode ierr = 0;
575 ierr = MatView(
matrix, viewer.viewer);
581 PetscErrorCode ierr = 0;
583 ierr = MatLoad(
matrix, viewer.viewer);
590 viewer.pushFormat(PETSC_VIEWER_ASCII_DENSE);
591 PetscErrorCode ierr = MatView(
matrix, viewer.viewer);
598 viewer.pushFormat(PETSC_VIEWER_ASCII_DENSE);
599 PetscErrorCode ierr = MatView(
matrix, viewer.viewer);
601 ierr = MatView(
matrix, viewer.viewer);
605 ierr = PetscViewerDrawGetDraw(viewer.viewer, 0, &draw);
607 ierr = PetscDrawSetPause(draw, -1);
615 PetscErrorCode ierr = 0;
623 PetscErrorCode ierr = 0;
624 PetscBool petscIsInitialized;
625 PetscInitialized(&petscIsInitialized);
626 if (petscIsInitialized &&
ksp)
627 ierr = KSPDestroy(&
ksp);
631KSPSolver::operator KSP &()
631KSPSolver::operator KSP &() {
…}
638 PetscErrorCode ierr = 0;
639 ierr = KSPReset(
ksp);
645 KSPConvergedReason convReason;
646 PetscErrorCode ierr = 0;
647 ierr = KSPGetConvergedReason(
ksp, &convReason);
651 if (convReason > 0) {
654 if (convReason == KSP_DIVERGED_ITS) {
669 KSPSolveTranspose(
ksp, b, x);
677 KSPConvergedReason convReason;
678 KSPGetConvergedReason(
ksp, &convReason);
680 PetscReal rtol, atol, dtol;
682 KSPGetTolerances(
ksp, &rtol, &atol, &dtol, &miter);
686 bool converged = (convReason >= 0);
687 bool stopped = (convReason == KSP_DIVERGED_ITS);
688 oss <<
"Solver " << (converged ?
"converged" : (stopped ?
"stopped" :
"diverged"));
692 switch (convReason) {
693 case (KSP_CONVERGED_RTOL):
694 case (KSP_CONVERGED_RTOL_NORMAL):
695 oss <<
" sufficient relative convergence";
697 case (KSP_CONVERGED_ATOL):
698 case (KSP_CONVERGED_ATOL_NORMAL):
699 oss <<
" sufficient absolute convergence";
701 case (KSP_DIVERGED_ITS):
702 oss <<
" reaching the maximum iterations";
704 case (KSP_DIVERGED_DTOL):
705 oss <<
" sufficient divergence";
707 case (KSP_DIVERGED_NANORINF):
708 oss <<
" the residual norm becoming nan or inf";
710 case (KSP_DIVERGED_BREAKDOWN):
711 oss <<
" a generic breakdown of the method";
714 oss <<
" the PETSc reason " << KSPConvergedReasons[convReason];
718 double bnorm = b.
l2norm();
719 double dlim = bnorm * dtol;
720 double rlim = bnorm * rtol;
722 oss <<
". Last residual norm: " <<
getResidualNorm() <<
", limits: relative " << rlim <<
" (rtol " << rtol <<
"), absolute " << atol <<
", divergence " << dlim <<
"(dtol " << dtol <<
')';
729 PetscErrorCode ierr = 0;
731 ierr = KSPGetIterationNumber(
ksp, &its);
738 PetscErrorCode ierr = 0;
740 ierr = KSPGetResidualNorm(
ksp, &val);
747 PetscErrorCode ierr = 0;
749 ierr = KSPGetTolerances(
ksp, &rtol,
nullptr,
nullptr,
nullptr);
758 PetscErrorCode ierr = 0;
759 PetscBool petscIsInitialized;
760 PetscInitialized(&petscIsInitialized);
762 if (IS and petscIsInitialized) {
763 ierr = ISLocalToGlobalMappingDestroy(IS);
770 PetscErrorCode ierr = 0;
771 PetscBool petscIsInitialized;
772 PetscInitialized(&petscIsInitialized);
774 if (ao and petscIsInitialized) {
775 ierr = AODestroy(ao);
#define PRECICE_TRACE(...)
int MPI_Comm_size(MPI_Comm comm, int *size)
#define PRECICE_ASSERT(...)
This class provides a lightweight logger.
A C++ 11 implementation of the non-owning C++20 std::span type.
constexpr iterator begin() const noexcept
constexpr iterator end() const noexcept
constexpr size_type size() const noexcept
static CommStatePtr current()
Returns an owning pointer to the current CommState.
static bool weInitialized
Whether we have initialized Petsc or if it was initialized by an application calling us.
static void finalize()
Finalizes Petsc environment.
static void initialize(utils::Parallel::Communicator comm)
Initializes the Petsc environment.
static logging::Logger _log
SolverResult getSolverResult()
Returns the current convergence reason as a SolverRestult.
PetscReal getRealtiveTolerance()
Returns the relavtive tolerance of the KSP.
std::string summaryFor(Vector &b)
Returns a summary the KSP solving for b.
SolverResult solveTranspose(Vector &b, Vector &x)
Solves the transposed linear system, returns false it not converged.
void reset()
Destroys and recreates the ksp on the same communicator.
KSPSolver(const KSPSolver &)=delete
Delete copy and assignment constructor.
PetscReal getResidualNorm()
Returns the last residual norm of the KSP.
PetscInt getIterationNumber()
Returns the iteration number of solver, either during or after the solve call.
SolverResult solve(Vector &b, Vector &x)
Solves the linear system, returns false it not converged.
SolverResult
The state of the KSP after returning from solve()
@ Diverged
The solver diverged.
@ Stopped
The solver reached the maximum iterations.
@ Converged
The solver converged.
void setColumn(Vector &v, PetscInt col)
void read(const std::string &filename)
Reads the matrix from file, stored in PETSc binary format.
void write(const std::string &filename, VIEWERFORMAT format=ASCII) const
Writes the matrix to file.
void assemble(MatAssemblyType type=MAT_FINAL_ASSEMBLY)
void view() const
Prints the matrix.
std::pair< PetscInt, PetscInt > getSize() const
Returns (rows, cols) global size.
void setValue(PetscInt row, PetscInt col, PetscScalar value)
void viewDraw() const
Graphically draws the matrix structure.
void init(PetscInt localRows, PetscInt localCols, PetscInt globalRows, PetscInt globalCols, MatType type=nullptr, bool doSetup=true)
Initializes matrix of given size and type.
std::pair< PetscInt, PetscInt > ownerRange() const
Returns a pair that mark the beginning and end of the matrix' ownership range.
PetscInt blockSize() const
Returns the block size of the matrix.
void reset()
Destroys and recreates the matrix on the same communicator.
Matrix(const Matrix &)=delete
Delete copy and assignment constructor.
std::pair< PetscInt, PetscInt > getLocalSize() const
Returns (rows, cols) local size.
MatInfo getInfo(MatInfoType flag) const
Get the MatInfo struct for the matrix.
std::pair< PetscInt, PetscInt > ownerRangeColumn() const
Returns a pair that mark the beginning and end of the matrix' column ownership range.
Vector & copyTo(precice::span< double > destination)
double l2norm() const
returns the l2-norm of the vector
PetscInt getLocalSize() const
void arange(double start, double stop)
void read(const std::string &filename, VIEWERFORMAT format=ASCII)
Reads the vector from file.
void setValue(PetscInt row, PetscScalar value)
void sort()
Sorts the LOCAL partition of the vector.
static Vector allocate(const std::string &name="")
Allocates a new vector on the given MPI communicator.
Vector & copyFrom(precice::span< const double > source)
void view() const
Prints the vector.
std::pair< PetscInt, PetscInt > ownerRange() const
Returns a pair that mark the beginning and end of the vectors ownership range. Use first and second t...
void write(const std::string &filename, VIEWERFORMAT format=ASCII) const
Writes the vector to file.
static logging::Logger _log
void swap(Vector &other) noexcept
Swaps the ownership of two vectors.
Vector & operator=(const Vector &other)
Vector(const std::string &name="")
Creates a new vector on the given MPI communicator.
void init(PetscInt rows)
Sets the size and calls VecSetFromOptions.
void setName(T obj, const std::string &name)
std::string getName(T obj)
MPI_Comm getCommunicator(T obj)
void destroy(ISLocalToGlobalMapping *IS)
Destroys an ISLocalToGlobalMapping, if IS is not null and PetscIsInitialized.
void swap(Vector &lhs, Vector &rhs) noexcept
contains precice-related utilities.
Viewer(const std::string &filename, VIEWERFORMAT format, MPI_Comm comm)
void pushFormat(PetscViewerFormat format)
Viewer(PetscViewerType type, MPI_Comm comm)