Petsc.hpp

Go to the documentation of this file.

#pragma once
 
#include "logging/Logger.hpp"
#include "precice/span.hpp"
#include "utils/Parallel.hpp"
 
namespace precice {
namespace logging {
class Logger;
} // namespace logging
 
namespace utils {
 
class Petsc {
public:
  static void initialize(utils::Parallel::Communicator comm);
 
  static void finalize();
 
private:
  static bool weInitialized;
 
  static logging::Logger _log;
};
} // namespace utils
} // namespace precice
 
#ifndef PRECICE_NO_PETSC
 
#include <string>
#include <utility>
#include "petscao.h"
#include "petscis.h"
#include "petscksp.h"
#include "petscmat.h"
#include "petscvec.h"
 
namespace precice {
namespace utils {
namespace petsc {
 
enum VIEWERFORMAT { ASCII,
                    BINARY };
 
class Matrix;
 
class Vector {
public:
  Vec vector = nullptr;
 
  enum LEFTRIGHT { LEFT,
                   RIGHT };
 
  explicit Vector(const std::string &name = "");
 
  Vector(const Vector &other);
 
  Vector &operator=(const Vector &other);
 
  Vector(Vector &&other) noexcept;
 
  Vector &operator=(Vector &&other) noexcept;
 
  Vector(Vec &other, const std::string &name = "");
 
  ~Vector();
 
 
  static Vector allocate(const std::string &name = "");
 
  static Vector allocate(Vector &other, const std::string &name = "");
 
  static Vector allocate(Vec &other, const std::string &name = "");
 
  static Vector allocate(Matrix &m, const std::string &name = "", LEFTRIGHT type = LEFT);
 
  static Vector allocate(Mat &m, const std::string &name = "", LEFTRIGHT type = LEFT);
 
 
  void swap(Vector &other) noexcept;
 
  operator Vec &();
 
  void init(PetscInt rows);
 
  PetscInt getSize() const;
 
  PetscInt getLocalSize() const;
 
  void setValue(PetscInt row, PetscScalar value);
 
  void arange(double start, double stop);
 
  void fillWithRandoms();
 
  Vector &copyFrom(precice::span<const double> source);
 
  Vector &copyTo(precice::span<double> destination);
 
  void sort();
 
  void assemble();
 
  std::pair<PetscInt, PetscInt> ownerRange() const;
 
  void write(const std::string &filename, VIEWERFORMAT format = ASCII) const;
 
  void read(const std::string &filename, VIEWERFORMAT format = ASCII);
 
  void view() const;
 
  double l2norm() const;
 
private:
  static logging::Logger _log;
};
 
void swap(Vector &lhs, Vector &rhs) noexcept;
 
class Matrix {
public:
  Mat matrix = nullptr;
 
 
  Matrix(const Matrix &) = delete;
  Matrix &operator=(const Matrix &) = delete;
 
  explicit Matrix(std::string name = "");
 
  Matrix(Matrix &&) = default;
  Matrix &operator=(Matrix &&) = default;
 
  ~Matrix();
 
  operator Mat &();
 
  void assemble(MatAssemblyType type = MAT_FINAL_ASSEMBLY);
 
 
  void init(PetscInt localRows, PetscInt localCols, PetscInt globalRows, PetscInt globalCols,
            MatType type = nullptr, bool doSetup = true);
 
  void reset();
 
 
  MatInfo getInfo(MatInfoType flag) const;
 
  void setValue(PetscInt row, PetscInt col, PetscScalar value);
 
  void fillWithRandoms();
 
  void setColumn(Vector &v, PetscInt col);
 
  std::pair<PetscInt, PetscInt> getSize() const;
 
  std::pair<PetscInt, PetscInt> getLocalSize() const;
 
  std::pair<PetscInt, PetscInt> ownerRange() const;
 
  std::pair<PetscInt, PetscInt> ownerRangeColumn() const;
 
  PetscInt blockSize() const;
 
  void write(const std::string &filename, VIEWERFORMAT format = ASCII) const;
 
  void read(const std::string &filename);
 
  void view() const;
 
  void viewDraw() const;
};
 
class KSPSolver {
public:
  KSP ksp;
 
 
  KSPSolver(const KSPSolver &) = delete;
  KSPSolver &operator=(const KSPSolver &) = delete;
 
  explicit KSPSolver(std::string name = "");
 
  KSPSolver(KSPSolver &&other) = default;
 
  ~KSPSolver();
 
  operator KSP &();
 
  void reset();
 
  enum struct SolverResult {
    Converged, 
    Stopped,   
    Diverged   
  };
 
  SolverResult getSolverResult();
 
  SolverResult solve(Vector &b, Vector &x);
 
  SolverResult solveTranspose(Vector &b, Vector &x);
 
  std::string summaryFor(Vector &b);
 
  PetscInt getIterationNumber();
 
  PetscReal getRealtiveTolerance();
 
  PetscReal getResidualNorm();
};
 
void destroy(KSP *ksp);
 
void destroy(ISLocalToGlobalMapping *IS);
 
void destroy(AO *ao);
 
} // namespace petsc
} // namespace utils
} // namespace precice
 
#endif // PRECICE_NO_PETSC