Merge branch 'NumLibMathLibCleanups' into 'master'

NumLib and MathLib PETSc related cleanups

See merge request ogs/ogs!4072

MathLib/LinAlg/LinAlg.cpp

@@ -29,7 +29,7 @@ void setLocalAccessibleVector(PETScVector const& x)
     x.setLocalAccessibleVector();
 }
 
-void set(PETScVector& x, double const a)
+void set(PETScVector& x, PetscScalar const a)
 {
     VecSet(x.getRawVector(), a);
 }
@@ -41,27 +41,28 @@ void copy(PETScVector const& x, PETScVector& y)
     VecCopy(x.getRawVector(), y.getRawVector());
 }
 
-void scale(PETScVector& x, double const a)
+void scale(PETScVector& x, PetscScalar const a)
 {
     VecScale(x.getRawVector(), a);
 }
 
 // y = a*y + X
-void aypx(PETScVector& y, double const a, PETScVector const& x)
+void aypx(PETScVector& y, PetscScalar const a, PETScVector const& x)
 {
     // TODO check sizes
     VecAYPX(y.getRawVector(), a, x.getRawVector());
 }
 
 // y = a*x + y
-void axpy(PETScVector& y, double const a, PETScVector const& x)
+void axpy(PETScVector& y, PetscScalar const a, PETScVector const& x)
 {
     // TODO check sizes
     VecAXPY(y.getRawVector(), a, x.getRawVector());
 }
 
 // y = a*x + b*y
-void axpby(PETScVector& y, double const a, double const b, PETScVector const& x)
+void axpby(PETScVector& y, PetscScalar const a, PetscScalar const b,
+           PETScVector const& x)
 {
     // TODO check sizes
     VecAXPBY(y.getRawVector(), a, b, x.getRawVector());
@@ -114,13 +115,13 @@ void copy(PETScMatrix const& A, PETScMatrix& B)
 }
 
 // A = a*A
-void scale(PETScMatrix& A, double const a)
+void scale(PETScMatrix& A, PetscScalar const a)
 {
     MatScale(A.getRawMatrix(), a);
 }
 
 // Y = a*Y + X
-void aypx(PETScMatrix& Y, double const a, PETScMatrix const& X)
+void aypx(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X)
 {
     // TODO check sizes
     // TODO sparsity pattern, currently they are assumed to be different (slow)
@@ -128,7 +129,7 @@ void aypx(PETScMatrix& Y, double const a, PETScMatrix const& X)
 }
 
 // Y = a*X + Y
-void axpy(PETScMatrix& Y, double const a, PETScMatrix const& X)
+void axpy(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X)
 {
     // TODO check sizes
     // TODO sparsity pattern, currently they are assumed to be different (slow)

MathLib/LinAlg/LinAlg.h

@@ -14,6 +14,10 @@
 #include "BaseLib/Error.h"
 #include "LinAlgEnums.h"
 
+#ifdef USE_PETSC
+#include <petscsystypes.h>
+#endif
+
 namespace MathLib
 {
 
@@ -40,29 +44,30 @@ void copy(MatrixOrVector const& x, MatrixOrVector& y)
 }
 
 //! Scales \c x by \c a
-template<typename MatrixOrVector>
+template <typename MatrixOrVector>
 void scale(MatrixOrVector& x, double const a)
 {
     x *= a;
 }
 
 //! Computes \f$ y = a \cdot y + x \f$.
-template<typename MatrixOrVector>
+template <typename MatrixOrVector>
 void aypx(MatrixOrVector& y, double const a, MatrixOrVector const& x)
 {
     y = a*y + x;
 }
 
 //! Computes \f$ y = a \cdot x + y \f$.
-template<typename MatrixOrVector>
+template <typename MatrixOrVector>
 void axpy(MatrixOrVector& y, double const a, MatrixOrVector const& x)
 {
     y += a*x;
 }
 
 //! Computes \f$ y = a \cdot x + b \cdot y \f$.
-template<typename MatrixOrVector>
-void axpby(MatrixOrVector& y, double const a, double const b, MatrixOrVector const& x)
+template <typename MatrixOrVector>
+void axpby(MatrixOrVector& y, double const a, double const b,
+           MatrixOrVector const& x)
 {
     y = a*x + b*y;
 }
@@ -150,35 +155,34 @@ namespace LinAlg
 /// Call this before call operator[] or get(...) of x.
 void setLocalAccessibleVector(PETScVector const& x);
 
-void set(PETScVector& x, double const a);
+void set(PETScVector& x, PetscScalar const a);
 
 void copy(PETScVector const& x, PETScVector& y);
 
-void scale(PETScVector& x, double const a);
+void scale(PETScVector& x, PetscScalar const a);
 
 // y = a*y + X
-void aypx(PETScVector& y, double const a, PETScVector const& x);
+void aypx(PETScVector& y, PetscScalar const a, PETScVector const& x);
 
 // y = a*x + y
-void axpy(PETScVector& y, double const a, PETScVector const& x);
+void axpy(PETScVector& y, PetscScalar const a, PETScVector const& x);
 
 // y = a*x + y
-void axpby(PETScVector& y, double const a, double const b, PETScVector const& x);
-
+void axpby(PETScVector& y, PetscScalar const a, PetscScalar const b,
+           PETScVector const& x);
 
 // Matrix
 
 void copy(PETScMatrix const& A, PETScMatrix& B);
 
 // A = a*A
-void scale(PETScMatrix& A, double const a);
+void scale(PETScMatrix& A, PetscScalar const a);
 
 // Y = a*Y + X
-void aypx(PETScMatrix& Y, double const a, PETScMatrix const& X);
+void aypx(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X);
 
 // Y = a*X + Y
-void axpy(PETScMatrix& Y, double const a, PETScMatrix const& X);
-
+void axpy(PETScMatrix& Y, PetscScalar const a, PETScMatrix const& X);
 
 // Matrix and Vector
 

MathLib/LinAlg/PETSc/PETScVector.cpp

@@ -314,7 +314,7 @@ void PETScVector::shallowCopy(const PETScVector& v)
     has_ghost_id_ = v.has_ghost_id_;
     global_ids2local_ids_ghost_ = v.global_ids2local_ids_ghost_;
 
-    VecSetOption(v_, VEC_IGNORE_NEGATIVE_INDICES, PETSC_TRUE);
+    config();
 }
 
 void finalizeVectorAssembly(PETScVector& vec)

NumLib/ODESolver/TimeDiscretization.cpp

@@ -47,4 +47,28 @@ double computeRelativeChangeFromPreviousTimestep(GlobalVector const& x,
     // Only norm_x is close to zero
     return norm_dx / std::numeric_limits<double>::epsilon();
 }
+
+void TimeDiscretization::getXdot(GlobalVector const& x_at_new_timestep,
+                                 GlobalVector const& x_old,
+                                 GlobalVector& xdot) const
+{
+    namespace LinAlg = MathLib::LinAlg;
+
+    double const dt = getCurrentTimeIncrement();
+
+    // xdot = 1/dt * x_at_new_timestep - x_old
+    getWeightedOldX(xdot, x_old);
+    LinAlg::axpby(xdot, 1. / dt, -1.0, x_at_new_timestep);
+}
+
+void BackwardEuler::getWeightedOldX(GlobalVector& y,
+                                    GlobalVector const& x_old) const
+{
+    namespace LinAlg = MathLib::LinAlg;
+
+    // y = x_old / delta_t
+    LinAlg::copy(x_old, y);
+    LinAlg::scale(y, 1.0 / _delta_t);
+}
+
 }  // end of namespace NumLib

NumLib/ODESolver/TimeDiscretization.h

@@ -132,16 +132,7 @@ public:
     //! \f$.
     void getXdot(GlobalVector const& x_at_new_timestep,
                  GlobalVector const& x_old,
-                 GlobalVector& xdot) const
-    {
-        namespace LinAlg = MathLib::LinAlg;
-
-        double const dt = getCurrentTimeIncrement();
-
-        // xdot = 1/dt * x_at_new_timestep - x_old
-        getWeightedOldX(xdot, x_old);
-        LinAlg::axpby(xdot, 1. / dt, -1.0, x_at_new_timestep);
-    }
+                 GlobalVector& xdot) const;
 
     //! Returns \f$ x_O \f$.
     virtual void getWeightedOldX(
@@ -155,7 +146,6 @@ class BackwardEuler final : public TimeDiscretization
 {
 public:
     void setInitialState(const double t0) override { _t = t0; }
-
     void nextTimestep(const double t, const double delta_t) override
     {
         _t = t;
@@ -165,14 +155,7 @@ public:
     double getCurrentTime() const override { return _t; }
     double getCurrentTimeIncrement() const override { return _delta_t; }
     void getWeightedOldX(GlobalVector& y,
-                         GlobalVector const& x_old) const override
-    {
-        namespace LinAlg = MathLib::LinAlg;
-
-        // y = x_old / delta_t
-        LinAlg::copy(x_old, y);
-        LinAlg::scale(y, 1.0 / _delta_t);
-    }
+                         GlobalVector const& x_old) const override;
 
 private:
     double _t = std::numeric_limits<double>::quiet_NaN();  //!< \f$ t_C \f$