Merge remote-tracking branch 'remotes/totto82/frankenstein_mod' into frankenstein_merge_master

* remotes/totto82/frankenstein_mod:
  Fix seg-fault for cases without wells
  Some micro performance improvments and cleaning
  Add THP support in the denseAD well model
  Only solve the linear system when it is not converged.
  Revert changes to NewtonIterationBlackoilInterleaved.cpp
  add and use class wellModelMatrixAdapter
  Remove unused code and remove Eigen vectors
  New updateState
  Some cleaning and small changes

opm/autodiff/NonlinearSolver.hpp

@@ -24,6 +24,9 @@
 #include <opm/autodiff/AutoDiffBlock.hpp>
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/simulator/SimulatorTimerInterface.hpp>
+#include <opm/autodiff/DuneMatrix.hpp>
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <memory>
 
 namespace Opm {
@@ -40,6 +43,10 @@ namespace Opm {
         typedef ADB::V V;
         typedef ADB::M M;
 
+        typedef double Scalar;
+        typedef Dune::FieldVector<Scalar, 3    >       VectorBlockType;
+        typedef Dune::BlockVector<VectorBlockType>      BVector;
+
         // Available relaxation scheme types.
         enum RelaxType { DAMPEN, SOR };
 
@@ -150,6 +157,8 @@ namespace Opm {
         /// Apply a stabilization to dx, depending on dxOld and relaxation parameters.
         void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const;
 
+        void stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const;
+
         /// The greatest relaxation factor (i.e. smallest factor) allowed.
         double relaxMax() const          { return param_.relax_max_; }
 

opm/autodiff/NonlinearSolver_impl.hpp

@@ -282,6 +282,38 @@ namespace Opm
         return;
     }
 
+    template <class PhysicalModel>
+    void
+    NonlinearSolver<PhysicalModel>::stabilizeNonlinearUpdate(BVector& dx, BVector& dxOld, const double omega) const
+    {
+        // The dxOld is updated with dx.
+        // If omega is equal to 1., no relaxtion will be appiled.
+
+        BVector tempDxOld = dxOld;
+        dxOld = dx;
+
+        switch (relaxType()) {
+            case DAMPEN:
+                if (omega == 1.) {
+                    return;
+                }
+                dx *= omega;
+                return;
+            case SOR:
+                if (omega == 1.) {
+                    return;
+                }
+                dx *= omega;
+                tempDxOld *= (1.-omega);
+                dx += tempDxOld;
+                return;
+            default:
+                OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type.");
+        }
+
+        return;
+    }
+
 
 } // namespace Opm
 

opm/autodiff/StandardWells_impl.hpp

@@ -672,17 +672,18 @@ namespace Opm
                         int table_id = well_controls_iget_vfp(wc, ctrl_index);
 
                         const WellType& well_type = wells().type[w];
+                        const int perf = wells().well_connpos[w]; //first perforation.
                         if (well_type == INJECTOR) {
                             double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getInj()->getTable(table_id)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
 
                             well_state.thp()[w] = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp);
                         }
                         else if (well_type == PRODUCER) {
                             double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getProd()->getTable(table_id)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
 
                             well_state.thp()[w] = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq);
                         }
@@ -786,17 +787,18 @@ namespace Opm
                 //Set *BHP* target by calculating bhp from THP
                 const WellType& well_type = wells().type[w];
 
+                const int perf = wells().well_connpos[w]; // first perforation.
                 if (well_type == INJECTOR) {
                     double dp = wellhelpers::computeHydrostaticCorrection(
                             wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
-                            wellPerforationDensities(), gravity_);
+                            wellPerforationDensities()[perf], gravity_);
 
                     xw.bhp()[w] = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
                 }
                 else if (well_type == PRODUCER) {
                     double dp = wellhelpers::computeHydrostaticCorrection(
                             wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
-                            wellPerforationDensities(), gravity_);
+                            wellPerforationDensities()[perf], gravity_);
 
                     xw.bhp()[w] = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
                 }
@@ -1096,11 +1098,12 @@ namespace Opm
 
                     //Set *BHP* target by calculating bhp from THP
                     const WellType& well_type = wells().type[w];
+                    const int perf = wells().well_connpos[w]; //first perforation
 
                     if (well_type == INJECTOR) {
                         double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getInj()->getTable(vfp)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
                         const double bhp = vfp_properties_->getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
                         // apply the strictest of the bhp controlls i.e. smallest bhp for injectors
                         if ( bhp < bhps[w]) {
@@ -1110,7 +1113,7 @@ namespace Opm
                     else if (well_type == PRODUCER) {
                         double dp = wellhelpers::computeHydrostaticCorrection(
                                     wells(), w, vfp_properties_->getProd()->getTable(vfp)->getDatumDepth(),
-                                    wellPerforationDensities(), gravity_);
+                                    wellPerforationDensities()[perf], gravity_);
 
                         const double bhp = vfp_properties_->getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
                         // apply the strictest of the bhp controlls i.e. largest bhp for producers

opm/autodiff/VFPHelpers.hpp

@@ -25,7 +25,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
-
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 
 /**
  * This file contains a set of helper functions used by VFPProd / VFPInj.
@@ -35,14 +36,7 @@ namespace detail {
 
 
 typedef AutoDiffBlock<double> ADB;
-
+typedef DenseAd::Evaluation<double, /*size=*/6> EvalWell;
-
-
-
-
-
-
-
 
 
 /**
@@ -52,7 +46,12 @@ inline double zeroIfNan(const double& value) {
     return (std::isnan(value)) ? 0.0 : value;
 }
 
-
+/**
+ * Returns zero if input value is NaN
+ */
+inline double zeroIfNan(const EvalWell& value) {
+    return (std::isnan(value.value)) ? 0.0 : value.value;
+}
 
 
 

opm/autodiff/VFPInjProperties.cpp

@@ -89,7 +89,36 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector<int>& table_id,
 
 
 
+VFPInjProperties::EvalWell VFPInjProperties::bhp(const int table_id,
+                                                   const EvalWell& aqua,
+                                                   const EvalWell& liquid,
+                                                   const EvalWell& vapour,
+                                                   const double& thp) const {
 
+    //Get the table
+    const VFPInjTable* table = detail::getTable(m_tables, table_id);
+    EvalWell bhp = 0.0;
+
+    //Find interpolation variables
+    EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
+
+    //Compute the BHP for each well independently
+    if (table != nullptr) {
+        //First, find the values to interpolate between
+        //Value of FLO is negative in OPM for producers, but positive in VFP table
+        auto flo_i = detail::findInterpData(flo.value, table->getFloAxis());
+        auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant
+
+        detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i);
+
+        bhp = bhp_val.dflo * flo;
+        bhp.value = bhp_val.value; // thp is assumed constant i.e.
+    }
+    else {
+        bhp.value = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
+    }
+    return bhp;
+}
 
 
 

opm/autodiff/VFPInjProperties.hpp

@@ -23,6 +23,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
 #include <opm/core/wells.h>
 #include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 
 #include <vector>
 #include <map>
@@ -91,6 +93,13 @@ public:
             const ADB& vapour,
             const ADB& thp) const;
 
+    typedef DenseAd::Evaluation<double, /*size=*/6> EvalWell;
+    EvalWell bhp(const int table_id,
+            const EvalWell& aqua,
+            const EvalWell& liquid,
+            const EvalWell& vapour,
+            const double& thp) const;
+
     /**
      * Linear interpolation of bhp as a function of the input parameters
      * @param table_id Table number to use

opm/autodiff/VFPProdProperties.cpp

@@ -24,7 +24,8 @@
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/common/ErrorMacros.hpp>
 #include <opm/autodiff/AutoDiffHelpers.hpp>
-
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 #include <opm/autodiff/VFPHelpers.hpp>
 
 
@@ -74,7 +75,42 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector<int>& table_id,
 }
 
 
+VFPProdProperties::EvalWell VFPProdProperties::bhp(const int table_id,
+                                                   const EvalWell& aqua,
+                                                   const EvalWell& liquid,
+                                                   const EvalWell& vapour,
+                                                   const double& thp,
+                                                   const double& alq) const {
 
+    //Get the table
+    const VFPProdTable* table = detail::getTable(m_tables, table_id);
+    EvalWell bhp = 0.0;
+
+    //Find interpolation variables
+    EvalWell flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
+    EvalWell wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
+    EvalWell gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
+
+    //Compute the BHP for each well independently
+    if (table != nullptr) {
+        //First, find the values to interpolate between
+        //Value of FLO is negative in OPM for producers, but positive in VFP table
+        auto flo_i = detail::findInterpData(-flo.value, table->getFloAxis());
+        auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); // assume constant
+        auto wfr_i = detail::findInterpData( wfr.value, table->getWFRAxis());
+        auto gfr_i = detail::findInterpData( gfr.value, table->getGFRAxis());
+        auto alq_i = detail::findInterpData( alq, table->getALQAxis()); //assume constant
+
+        detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
+
+        bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (bhp_val.dflo * flo);
+        bhp.value = bhp_val.value;
+    }
+    else {
+        bhp.value = -1e100; //Signal that this value has not been calculated properly, due to "missing" table
+    }
+    return bhp;
+}
 
 
 

opm/autodiff/VFPProdProperties.hpp

@@ -23,6 +23,8 @@
 #include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
 #include <opm/core/wells.h>
 #include <opm/autodiff/AutoDiffBlock.hpp>
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/densead/Evaluation.hpp>
 
 #include <vector>
 #include <map>
@@ -99,6 +101,15 @@ public:
             const ADB& thp,
             const ADB& alq) const;
 
+
+    typedef DenseAd::Evaluation<double, /*size=*/6> EvalWell;
+    EvalWell bhp(const int table_id,
+            const EvalWell& aqua,
+            const EvalWell& liquid,
+            const EvalWell& vapour,
+            const double& thp,
+            const double& alq) const;
+
     /**
      * Linear interpolation of bhp as a function of the input parameters
      * @param table_id Table number to use

opm/autodiff/WellHelpers.hpp

@@ -131,7 +131,7 @@ namespace Opm {
          */
         inline
         double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth,
-                                            const Vector& well_perforation_densities, const double gravity) {
+                                            const double& rho, const double gravity) {
             if ( wells.well_connpos[w] == wells.well_connpos[w+1] )
             {
                 // This is a well with no perforations.
@@ -142,13 +142,12 @@ namespace Opm {
             }
             const double well_ref_depth = wells.depth_ref[w];
             const double dh = vfp_ref_depth - well_ref_depth;
-            const int perf = wells.well_connpos[w];
-            const double rho = well_perforation_densities[perf];
             const double dp = rho*gravity*dh;
 
             return dp;
         }
 
+
         inline
         Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth,
                 const Vector& well_perforation_densities, const double gravity) {
@@ -157,7 +156,23 @@ namespace Opm {
 
 #pragma omp parallel for schedule(static)
             for (int i=0; i<nw; ++i) {
-                retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities, gravity);
+                const int perf = wells.well_connpos[i];
+                retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
+            }
+
+            return retval;
+        }
+
+        inline
+        std::vector<double> computeHydrostaticCorrection(const Wells& wells, const std::vector<double>& vfp_ref_depth,
+                                                         const std::vector<double>& well_perforation_densities, const double gravity) {
+            const int nw = wells.number_of_wells;
+            std::vector<double> retval(nw,0.0);
+
+#pragma omp parallel for schedule(static)
+            for (int i=0; i<nw; ++i) {
+                const int perf = wells.well_connpos[i];
+                retval[i] = computeHydrostaticCorrection(wells, i, vfp_ref_depth[i], well_perforation_densities[perf], gravity);
             }
 
             return retval;

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -119,6 +119,7 @@ namespace Opm
                 const WellType& well_type = wells->type[w];
 
                 switch (well_controls_iget_type(wc, current)) {
+                case THP: // Intentional fall-through
                 case BHP:
                 {
                     if (well_type == INJECTOR) {
@@ -148,13 +149,13 @@ namespace Opm
                     total_rates += g[p] * wellRates()[np*w + p];
                 }
 
-                //if(std::abs(total_rates) > 0) {
+                if(std::abs(total_rates) > 0) {
-                //    wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water;
+                    wellSolutions()[nw + w] = g[Water] * wellRates()[np*w + Water] / total_rates; //wells->comp_frac[np*w + Water]; // Water;
-                //    wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas
+                    wellSolutions()[2*nw + w] = g[Gas] * wellRates()[np*w + Gas] / total_rates ; //wells->comp_frac[np*w + Gas]; //Gas
-                //} else {
+                } else {
                     wellSolutions()[nw + w] = wells->comp_frac[np*w + Water];
                     wellSolutions()[2*nw + w] = wells->comp_frac[np*w + Gas];
-                //}
+                }
 
 
             }