implementing two apply() function in StandardWell_impl

The two functions will be essentially the same even for different types
of wells. Maybe later we should try to put them in WellInterface.

opm/autodiff/StandardWell.hpp

@@ -25,13 +25,6 @@
 
 #include <opm/autodiff/WellInterface.hpp>
 
-#include<dune/common/fmatrix.hh>
-#include<dune/istl/bcrsmatrix.hh>
-#include<dune/istl/matrixmatrix.hh>
-
-#include <opm/material/densead/Math.hpp>
-#include <opm/material/densead/Evaluation.hpp>
-
 namespace Opm
 {
 
@@ -64,20 +57,18 @@ namespace Opm
             SFrac = 3
         };
 
+        using WellInterface<TypeTag>::numEq;
+        using typename WellInterface<TypeTag>::VectorBlockType;
+        using typename WellInterface<TypeTag>::MatrixBlockType;
+        using typename WellInterface<TypeTag>::Mat;
+        using typename WellInterface<TypeTag>::BVector;
+        using typename WellInterface<TypeTag>::Eval;
 
-        typedef double Scalar;
-        // static const int numEq = BlackoilIndices::numEq;
-        static const int numEq = BlackoilIndices::numEq;
         static const int numWellEq = GET_PROP_VALUE(TypeTag, EnablePolymer)? 3:numEq; // //numEq; //number of wellEq is only for 3 for polymer
         static const int contiSolventEqIdx = BlackoilIndices::contiSolventEqIdx;
         static const int contiPolymerEqIdx = BlackoilIndices::contiPolymerEqIdx;
 
-        typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
-        typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
-        typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
-        typedef Dune::BlockVector<VectorBlockType> BVector;
         typedef DenseAd::Evaluation<double, /*size=*/numEq + numWellEq> EvalWell;
-        typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
 
         StandardWell(const Well* well, const int time_step, const Wells* wells);
 
@@ -145,6 +136,11 @@ namespace Opm
         virtual void computeWellConnectionPressures(const Simulator& ebosSimulator,
                                                     const WellState& xw);
 
+        // Ax = Ax - C D^-1 B x
+        virtual void apply(const BVector& x, BVector& Ax) const;
+        // r = r - C D^-1 Rw
+        virtual void apply(BVector& r) const;
+
         using WellInterface<TypeTag>::phaseUsage;
         using WellInterface<TypeTag>::active;
         using WellInterface<TypeTag>::numberOfPerforations;

opm/autodiff/StandardWell_impl.hpp

@@ -1864,4 +1864,44 @@ namespace Opm
         }
     }
 
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    apply(const BVector& x, BVector& Ax) const
+    {
+        assert( Bx_.size() == duneB_.N() );
+        assert( invDrw_.size() == invDuneD_.N() );
+
+        // Bx_ = duneB_ * x
+        duneB_.mv(x, Bx_);
+        // invDBx = invDuneD_ * Bx_
+        // TODO: with this, we modified the content of the invDrw_.
+        // Is it necessary to do this to save some memory?
+        BVector& invDBx = invDrw_;
+        invDuneD_.mv(Bx_, invDBx);
+
+        // Ax = Ax - duneC_^T * invDBx
+        duneC_.mmtv(invDBx,Ax);
+    }
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    apply(BVector& r) const
+    {
+        assert( invDrw_.size() == invDuneD_.N() );
+
+        // invDrw_ = invDuneD_ * resWell_
+        invDuneD_.mv(resWell_, invDrw_);
+        // r = r - duneC_^T * invDrw_
+        duneC_.mmtv(invDrw_, r);
+    }
+
 }

opm/autodiff/StandardWellsDense_impl.hpp

@@ -377,6 +377,13 @@ namespace Opm {
 
 
 
+    // applying the well residual to reservoir residuals
+    // r = r - duneC_^T * invDuneD_ * resWell_
+    // TODO: for this, we should calcuate the duneC_^T * invDuneD_ * resWell_ for each
+    // well, then sum them up and apply to r finally
+    // In a more general case, the number of the equations for reservoir and wells can be different,
+    // we need to think about the possible data types can be faced.
+    // we do not want to expose the some well related data type even inside the Well Model
     template<typename TypeTag>
     void
     StandardWellsDense<TypeTag>::
@@ -402,16 +409,23 @@ namespace Opm {
             return;
         }
 
-        assert( invDrw_.size() == invDuneD_.N() );
+        for (auto& well : well_container_) {
+            well->apply(r);
+        }
 
+        /* assert( invDrw_.size() == invDuneD_.N() );
+
+        // invDrw_ = invDuneD_ * resWell_
         invDuneD_.mv(resWell_,invDrw_);
-        duneC_.mmtv(invDrw_, r);
+        // r = r - duneC_^T * invDrw_
+        duneC_.mmtv(invDrw_, r); */
     }
 
 
 
 
 
+    // Ax = A x - C D^-1 B x
     template<typename TypeTag>
     void
     StandardWellsDense<TypeTag>::
@@ -421,20 +435,33 @@ namespace Opm {
             return;
         }
 
-        assert( Bx_.size() == duneB_.N() );
+        for (auto& well : well_container_) {
+            well->apply(x, Ax);
+        }
+
+        /* assert( Bx_.size() == duneB_.N() );
 
         BVector& invDBx = invDrw_;
         assert( invDBx.size() == invDuneD_.N());
 
+        // Bx_ = duneB_ * x
         duneB_.mv(x, Bx_);
+        // invDBx = invDuneD_ * Bx_
         invDuneD_.mv(Bx_, invDBx);
+        // Ax = Ax - duneC_^T * invDBx
         duneC_.mmtv(invDBx,Ax);
+        */
     }
 
 
 
 
 
+    // Ax = Ax - alpha * C D^-1 B x
+    // TODO: for the new Well Model, we will calcuate
+    // C D^-1 B for each well and sum it up
+    // while it can be implemented in the function apply()
+    // then this function does not need to change
     template<typename TypeTag>
     void
     StandardWellsDense<TypeTag>::
@@ -449,7 +476,9 @@ namespace Opm {
         }
 
         scaleAddRes_ = 0.0;
+        // scaleAddRes_  = - C D^-1 B x
         apply( x, scaleAddRes_ );
+        // Ax = Ax + alpha * scaleAddRes_
         Ax.axpy( alpha, scaleAddRes_ );
     }
 
@@ -457,6 +486,10 @@ namespace Opm {
 
 
 
+    // xw = D^-1(resWell - B * x)
+    // TODO: this function should be moved to StandardWell
+    // xw should not appear in StandardWellsDense to avoid
+    // the data type to hold different types of xw
     template<typename TypeTag>
     void
     StandardWellsDense<TypeTag>::
@@ -466,7 +499,9 @@ namespace Opm {
              return;
         }
         BVector resWell = resWell_;
+        // resWell = resWell - B * x
         duneB_.mmv(x, resWell);
+        // xw = D^-1 * resWell
         invDuneD_.mv(resWell, xw);
     }
 

opm/autodiff/WellInterface.hpp

@@ -39,6 +39,13 @@
 
 #include <opm/simulators/WellSwitchingLogger.hpp>
 
+#include<dune/common/fmatrix.hh>
+#include<dune/istl/bcrsmatrix.hh>
+#include<dune/istl/matrixmatrix.hh>
+
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+
 #include <string>
 #include <memory>
 #include <vector>
@@ -62,7 +69,15 @@ namespace Opm
         typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
         typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
 
-        static const int solventCompIdx = 3; //TODO get this from ebos
+        static const int numEq = BlackoilIndices::numEq;
+        typedef double Scalar;
+
+        typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
+        typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
+        typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
+        typedef Dune::BlockVector<VectorBlockType> BVector;
+        typedef DenseAd::Evaluation<double, /*size=*/numEq> Eval;
+
         static const bool has_solvent = GET_PROP_VALUE(TypeTag, EnableSolvent);
 
         /// Constructor
@@ -171,6 +186,12 @@ namespace Opm
         virtual void computeWellConnectionPressures(const Simulator& ebosSimulator,
                                                     const WellState& xw) = 0;
 
+        // Ax = Ax - C D^-1 B x
+        virtual void apply(const BVector& x, BVector& Ax) const = 0;
+
+        // r = r - C D^-1 Rw
+        virtual void apply(BVector& r) const = 0;
+
     protected:
         // TODO: some variables shared by all the wells should be made static
         // well name