Add LinearizationType, which will be used for sequential methods.

opm/models/blackoil/blackoilenergymodules.hh

@@ -358,7 +358,7 @@ public:
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
 
         // set temperature
-        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx));
+        fs.setTemperature(priVars.makeEvaluation(temperatureIdx, timeIdx, elemCtx.linearizationType()));
     }
 
     /*!

opm/models/blackoil/blackoilpolymermodules.hh

@@ -884,6 +884,7 @@ public:
         // Use log(v0) as initial value for u
         auto u = v0AbsLog;
         bool converged = false;
+        // TODO make this into parameters
         for (int i = 0; i < 20; ++i) {
             auto f = F(u);
             auto df = dF(u);
@@ -1044,10 +1045,11 @@ public:
                                   unsigned dofIdx,
                                   unsigned timeIdx)
     {
+        const auto linearizationType = elemCtx.linearizationType();
         const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
-        polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx);
+        polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx, linearizationType);
         if (enablePolymerMolarWeight) {
-            polymerMoleWeight_ = priVars.makeEvaluation(polymerMoleWeightIdx, timeIdx);
+            polymerMoleWeight_ = priVars.makeEvaluation(polymerMoleWeightIdx, timeIdx, linearizationType);
         }
         const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx);
 

opm/models/blackoil/blackoilsolventmodules.hh

@@ -902,7 +902,7 @@ public:
     {
         const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
         auto& fs = asImp_().fluidState_;
-        solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx);
+        solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType());
         hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
 
         // apply a cut-off. Don't waste calculations if no solvent
@@ -951,10 +951,11 @@ public:
             MaterialLaw::capillaryPressures(pC, materialParams, fs);
 
             //oil is the reference phase for pressure
+            const auto linearizationType = elemCtx.linearizationType();
             if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv)
-                pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+                pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, linearizationType);
             else {
-                const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
+                const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx, linearizationType);
                 pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
             }
 

opm/models/discretization/common/fvbasediscretization.hh

@@ -806,6 +806,15 @@ public:
     bool enableStorageCache() const
     { return enableStorageCache_; }
 
+    /*!
+     * \brief Set the value of enable storage cache
+     *
+     * Be aware that calling the *CachedStorage() methods if the storage cache is
+     * disabled will crash the program.
+     */
+    void setEnableStorageCache(bool enableStorageCache)
+    { enableStorageCache_= enableStorageCache; }
+
     /*!
      * \brief Retrieve an entry of the cache for the storage term.
      *

opm/models/discretization/common/fvbaseelementcontext.hh

@@ -30,6 +30,7 @@
 
 #include "fvbaseproperties.hh"
 
+#include <opm/models/discretization/common/linearizationtype.hh>
 #include <opm/models/utils/alignedallocator.hh>
 
 #include <opm/material/common/Unused.hpp>
@@ -268,6 +269,14 @@ public:
     unsigned focusDofIndex() const
     { return focusDofIdx_; }
 
+    /*!
+     * \brief Returns the linearization type.
+     *
+     * \copydetails setLinearizationType()
+     */
+    LinearizationType linearizationType() const
+    { return this->model().linearizer().getLinearizationType(); }
+
     /*!
      * \brief Return a reference to the simulator.
      */

opm/models/discretization/common/fvbaselinearizer.hh

@@ -29,6 +29,7 @@
 #define EWOMS_FV_BASE_LINEARIZER_HH
 
 #include "fvbaseproperties.hh"
+#include "linearizationtype.hh"
 
 #include <opm/models/parallel/gridcommhandles.hh>
 #include <opm/models/parallel/threadmanager.hh>
@@ -143,6 +144,12 @@ public:
     {
         simulatorPtr_ = &simulator;
         eraseMatrix();
+        auto it = elementCtx_.begin();
+        const auto& endIt = elementCtx_.end();
+        for (; it != endIt; ++it){
+            delete *it;
+        }
+        elementCtx_.resize(0);
     }
 
     /*!
@@ -160,7 +167,11 @@ public:
     /*!
      * \brief Linearize the full system of non-linear equations.
      *
-     * This means the spatial domain plus all auxiliary equations.
+     * The linearizationType() controls the scheme used and the focus
+     * time index. The default is fully implicit scheme, and focus index
+     * equal to 0, i.e. current time (end of step).
+     *
+     * This linearizes the spatial domain and all auxiliary equations.
      */
     void linearize()
     {
@@ -263,6 +274,14 @@ public:
     GlobalEqVector& residual()
     { return residual_; }
 
+    void setLinearizationType(LinearizationType linearizationType){
+        linearizationType_ = linearizationType;
+    };
+
+    const LinearizationType& getLinearizationType() const{
+        return linearizationType_;
+    };
+
     /*!
      * \brief Returns the map of constraint degrees of freedom.
      *
@@ -573,6 +592,7 @@ private:
     // the right-hand side
     GlobalEqVector residual_;
 
+    LinearizationType linearizationType_;
 
     std::mutex globalMatrixMutex_;
 };

opm/models/discretization/common/fvbaseprimaryvariables.hh

@@ -31,7 +31,7 @@
 #include <type_traits>
 
 #include "fvbaseproperties.hh"
-
+#include "linearizationtype.hh"
 #include <opm/material/common/Valgrind.hpp>
 #include <opm/material/common/Unused.hpp>
 #include <opm/material/common/Exceptions.hpp>
@@ -87,13 +87,13 @@ public:
      * it represents the a constant f = x_i. (the difference is that in the first case,
      * the derivative w.r.t. x_i is 1, while it is 0 in the second case.
      */
-    Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx) const
+    Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx, LinearizationType linearizationType = LinearizationType()) const
     {
         if (std::is_same<Evaluation, Scalar>::value)
             return (*this)[varIdx]; // finite differences
         else {
             // automatic differentiation
-            if (timeIdx == 0)
+            if (timeIdx == linearizationType.time)
                 return Toolbox::createVariable((*this)[varIdx], varIdx);
             else
                 return Toolbox::createConstant((*this)[varIdx]);

opm/models/discretization/common/linearizationtype.hh

@@ -0,0 +1,44 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::FvBaseLinearizer
+ */
+#ifndef EWOMS_FV_BASE_LINEARIZATIONTYPE_HH
+#define EWOMS_FV_BASE_LINEARIZATIONTYPE_HH
+
+namespace Opm
+{
+
+struct LinearizationType
+{
+    enum VarType { implicit, pressure, seqtransport };
+    VarType type = implicit;
+    unsigned int time = 0;
+};
+
+} // namespace Opm
+
+
+#endif