Mixed-precision improvements

CMakeLists_files.cmake

@@ -278,6 +278,7 @@ list(APPEND PRIVATE_HEADER_FILES
 
 if (HAVE_AVX2_EXTENSION)
   set (AVX2_SOURCE_FILES
+    opm/simulators/linalg/mixed/dot.c
     opm/simulators/linalg/mixed/bsr.c
     opm/simulators/linalg/mixed/prec.c
     opm/simulators/linalg/mixed/bslv.c)

opm/simulators/linalg/FlexibleSolver_impl.hpp

@@ -219,13 +219,13 @@ namespace Dune
                                                                             maxiter, // maximum number of iterations
                                                                             verbosity);
 #if HAVE_AVX2_EXTENSION
-          } else if (solver_type == "mixed-bicgstab") {
+          } else if (solver_type == "mixed-legacy") {
               if constexpr (Opm::is_gpu_operator_v<Operator>) {
-                OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for GPU operators");
+                OPM_THROW(std::invalid_argument, "legacy mixed-bicgstab solver not supported for GPU operators");
             } else if constexpr (Opm::detail::is_multi_type_block_vector_v<VectorType>) {
-                OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for multi-type block vectors.");
+                OPM_THROW(std::invalid_argument, "legacy mixed-bicgstab solver not supported for multi-type block vectors.");
             } else if constexpr (std::is_same_v<typename VectorType::field_type, float>){
-                OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for single precision.");
+                OPM_THROW(std::invalid_argument, "legacy mixed-bicgstab solver not supported for single precision.");
             } else {
                 const std::string prec_type = prm.get<std::string>("preconditioner.type", "error");
                 bool use_mixed_dilu= (prec_type=="legacy-mixed-dilu");
@@ -238,13 +238,13 @@ namespace Dune
                                                                         );
             }
             // MixedBiCGSTABSolver starts here
-          } else if (solver_type == "mixed-precision") {
+          } else if (solver_type == "mixed-bicgstab") {
               if constexpr (Opm::is_gpu_operator_v<Operator>) {
-                OPM_THROW(std::invalid_argument, "mixed-precision solver not supported for GPU operators");
+                OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for GPU operators");
             } else if constexpr (Opm::detail::is_multi_type_block_vector_v<VectorType>) {
                 OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for multi-type block vectors.");
             } else if constexpr (std::is_same_v<typename VectorType::field_type, float>){
-                OPM_THROW(std::invalid_argument, "mixed-precision solver not supported for single precision.");
+                OPM_THROW(std::invalid_argument, "mixed-bicgstab solver not supported for single precision.");
             } else {
                 linsolver_ = std::make_shared<Dune::MixedBiCGSTABSolver<Comm,Operator,VectorType>>(linearoperator_for_solver_,
                                                                             scalarproduct_,

opm/simulators/linalg/ScalarProducts.hpp

@@ -0,0 +1,194 @@
+#ifndef OPM_SCALAR_PRODUCTS_HEADER_INCLUDED
+#define OPM_SCALAR_PRODUCTS_HEADER_INCLUDED
+
+#include <opm/simulators/linalg/mixed/dot.h>
+
+namespace Dune
+{
+
+/// A parallel scalar product that takes advantage of the fact that all
+/// elements associated with ghost cells are located at the end of the
+/// vector. This allows us to ignore the block structure of the vector
+/// and eliminate the use of a mask to exclude ghost entries from being
+/// included in the scalar product
+template<class Vector, class Comm>
+class GhostLastScalarProduct : public ScalarProduct<Vector>
+{
+    public:
+
+    ///Exctract block size from vector type
+    static constexpr auto block_size = Vector::block_type::dimension;
+
+    /*! \brief constructor
+    * \param com The communication object for syncing overlap and copy
+    * data points.
+    * \param cat parallel solver category (nonoverlapping or overlapping)
+    */
+    GhostLastScalarProduct (std::shared_ptr<const Comm> com, SolverCategory::Category cat)
+        : _communication(com), _category(cat)
+    {
+        count_ = getLocalCount();        // number or local cells
+        int verify = verifyLocalCount(); // redundant check on numbef of local cells
+        if (count_ != verify) OPM_THROW(std::runtime_error, "Inconsistent local node count!!\n");
+    }
+
+    /*! \brief constructor
+    * \param com The communication object for syncing overlap and copy
+    * data points.
+    * \param cat parallel solver category (nonoverlapping or overlapping)
+    * \note if you use this constructor you have to make sure com stays alive
+    */
+    GhostLastScalarProduct (const Comm& com, SolverCategory::Category cat)
+        : GhostLastScalarProduct(stackobject_to_shared_ptr(com), cat)
+    {}
+
+    /*! \brief Dot product of two vectors.
+    * \param vx first input vector
+    * \param vy second input vector
+    */
+    virtual double dot (const Vector& vx, const Vector& vy) const override
+    {
+
+        // access underlying data
+        double const *x = &vx[0][0];
+        double const *y = &vy[0][0];
+
+        // total array length
+        int NN = block_size*count_;
+#if 0
+        // unroll loop in multiples of 8
+        int n=NN/8;
+        int N=8*n;
+        double agg[8];
+        for(int i=0;i<8;i++) agg[i]=0.0;
+        for(int i=0;i<N;i+=8) for(int j=0;j<8;j++) agg[j]+=x[i+j]*y[i+j];
+        for(int j=0;j<4;j++) agg[j]+=agg[j+4];
+        for(int j=0;j<2;j++) agg[j]+=agg[j+2];
+        for(int j=0;j<1;j++) agg[j]+=agg[j+1];
+
+        // loop-peeling of trailing end
+        for(int j=N;j<NN;j++) agg[0]+=x[j]*y[j];
+
+        // Global summation
+        auto cc = _communication->communicator();
+        double result = cc.sum(agg[0]);
+        return result;
+#else
+        auto cc = _communication->communicator();
+        return cc.sum(vec_dot(x,y,NN));
+#endif
+        //return cc.sum(vec_dot(x,y,NN));
+    }
+
+    /*! \brief Vector L2-norm.
+    * \param vx input vector
+    */
+    virtual double norm (const Vector& vx) const override
+    {
+        return sqrt(dot(vx,vx));
+    }
+
+    //! Category of the scalar product (see SolverCategory::Category)
+    virtual SolverCategory::Category category() const override
+    {
+        return _category;
+    }
+
+    private:
+    std::shared_ptr<const Comm> _communication;
+    SolverCategory::Category _category;
+    int count_;
+
+    /*! \brief Count number of local cells.
+    */
+    int getLocalCount() const
+    {
+        int count = 0;
+        // Loop over index set
+        auto indexSet = _communication->indexSet();
+        for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) {
+            if (idx->local().attribute()==1) count++; // count non-local indices
+        }
+        return count;
+    }
+
+    /*! \brief Infer number of local cells from largest local index.
+    */
+    int verifyLocalCount() const
+    {
+        auto indexSet = _communication->indexSet();
+
+        size_t is = 0;
+        // Loop over index set
+        for (auto idx = indexSet.begin(); idx!=indexSet.end(); ++idx) {
+            //Only take "owner" indices
+            if (idx->local().attribute()==1) {
+                //get local index
+                auto loc = idx->local().local();
+                // if loc is higher than "old interior size", update it
+                if (loc > is) {
+                    is = loc;
+                }
+            }
+        }
+        return is + 1; //size is plus 1 since we start at 0
+    }
+
+};
+
+
+
+/// A sequential scalar product that ignores block structure of the vector
+/// to facilitate well-known optimization techniques
+template<class Vector>
+class SeqOptmizedProduct : public Dune::SeqScalarProduct<Vector>
+{
+public:
+
+    // extract block size
+    static constexpr auto block_size = Vector::block_type::dimension;
+
+    /*! \brief Dot product of two vectors.
+    * \param vx first input vector
+    * \param vy second input vector
+    */
+    virtual double dot(const Vector& vx, const Vector& vy) const override
+    {
+        // access underlying data
+        double const *x = &vx[0][0];
+        double const *y = &vy[0][0];
+
+        // total array length
+        int NN = block_size*vx.N();
+#if 1
+        // unroll loop in multiples of 8
+        int n=NN/8;
+        int N=8*n;
+        double agg[8];
+        for(int i=0;i<8;i++) agg[i]=0.0;
+        for(int i=0;i<N;i+=8) for(int j=0;j<8;j++) agg[j]+=x[i+j]*y[i+j];
+        for(int j=0;j<4;j++) agg[j]+=agg[j+4];
+        for(int j=0;j<2;j++) agg[j]+=agg[j+2];
+        for(int j=0;j<1;j++) agg[j]+=agg[j+1];
+
+        // loop-peeling of trailing end
+        for(int j=N;j<NN;j++) agg[0]+=x[j]*y[j];
+
+        return agg[0];
+#else
+        return vec_dot(x,y,NN);
+#endif
+    }
+
+    /*! \brief Vector L2-norm.
+    * \param vx input vector
+    */
+    virtual double norm(const Vector& vx) const override {
+        return std::sqrt(this->dot(vx, vx));
+    }
+};
+
+}
+
+#endif //OPM_SCALAR_PRODUCTS_HEADER_INCLUDED
+

opm/simulators/linalg/WellOperators.hpp

@@ -262,6 +262,8 @@ class WellModelMatrixAdapter : public Dune::AssembledLinearOperator<M,X,Y>
 
     const matrix_type& getmat() const override { return A_; }
 
+    const LinearOperatorExtra<X, Y>& getwellOper() const { return wellOper_; }
+
     void addWellPressureEquations(PressureMatrix& jacobian,
                                   const X& weights,
                                   const bool use_well_weights) const
@@ -357,6 +359,8 @@ class WellModelGhostLastMatrixAdapter : public Dune::AssembledLinearOperator<M,X
 
     const matrix_type& getmat() const override { return A_; }
 
+    const LinearOperatorExtra<X, Y>& getwellOper() const { return wellOper_; }
+
     void addWellPressureEquations(PressureMatrix& jacobian,
                                   const X& weights,
                                   const bool use_well_weights) const