Add TPFA-specific linearizer variant.

opm/models/discretization/common/fvbasediscretization.hh

@@ -321,6 +321,10 @@ struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization> { static const
 template<class TypeTag>
 struct EnableExperiments<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = true; };
 
+//! Default to not using the specialized TPFA linearizer.
+template<class TypeTag>
+struct UseTpfaLinearizer<TypeTag, TTag::FvBaseDiscretization> { static constexpr bool value = false; };
+
 } // namespace Opm::Properties
 
 namespace Opm {

opm/models/discretization/common/fvbaselinearizer.hh

@@ -37,6 +37,7 @@
 #include <opm/models/discretization/common/baseauxiliarymodule.hh>
 
 #include <opm/material/common/Exceptions.hpp>
+#include <opm/grid/utility/SparseTable.hpp>
 
 #include <dune/common/version.hh>
 #include <dune/common/fvector.hh>
@@ -86,6 +87,8 @@ class FvBaseLinearizer
     using Constraints = GetPropType<TypeTag, Properties::Constraints>;
     using Stencil = GetPropType<TypeTag, Properties::Stencil>;
     using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
+    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
 
     using GridCommHandleFactory = GetPropType<TypeTag, Properties::GridCommHandleFactory>;
 
@@ -103,8 +106,11 @@ class FvBaseLinearizer
 
     using MatrixBlock = typename SparseMatrixAdapter::MatrixBlock;
     using VectorBlock = Dune::FieldVector<Scalar, numEq>;
+    using ADVectorBlock = GetPropType<TypeTag, Properties::RateVector>;
+
 
     static const bool linearizeNonLocalElements = getPropValue<TypeTag, Properties::LinearizeNonLocalElements>();
+    static const bool useTpfaLinearizer = getPropValue<TypeTag, Properties::UseTpfaLinearizer>();
 
     // copying the linearizer is not a good idea
     FvBaseLinearizer(const FvBaseLinearizer&);
@@ -368,6 +374,26 @@ private:
 
         // create matrix structure based on sparsity pattern
         jacobian_->reserve(sparsityPattern);
+
+        for (unsigned globI = 0; globI < model.numTotalDof(); globI++) {
+            sparsityPattern[globI].erase(globI);
+        }
+        unsigned numCells = model.numTotalDof();
+        neighbours_.reserve(numCells, 6 * numCells);
+        trans_.reserve(numCells, 6 * numCells);
+        std::vector<double> loctrans;
+        for (unsigned globI = 0; globI < numCells; globI++) {
+            const auto& cells = sparsityPattern[globI];
+            neighbours_.appendRow(cells.begin(), cells.end());
+            unsigned n = cells.size();
+            loctrans.resize(n);
+            short loc = 0;
+            for (const int& cell : cells) {
+                loctrans[loc] = problem_().transmissibility(globI, cell);
+                loc++;
+            }
+            trans_.appendRow(loctrans.begin(), loctrans.end());
+        }
     }
 
     // reset the global linear system of equations.
@@ -424,9 +450,103 @@ private:
         }
     }
 
+public:
+    void setResAndJacobi(VectorBlock& res, MatrixBlock& bMat, const ADVectorBlock& resid) const
+    {
+        for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++)
+            res[eqIdx] = resid[eqIdx].value();
+
+        for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
+            for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
+                // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of
+                // the residual function 'eqIdx' for the degree of freedom 'dofIdx'
+                // with regard to the focus variable 'pvIdx' of the degree of freedom
+                // 'focusDofIdx'
+                bMat[eqIdx][pvIdx] = resid[eqIdx].derivative(pvIdx);
+            }
+        }
+    }
+
+private:
+    void linearizeGlobalTPFA_()
+    {
+        const bool well_local = false;
+        resetSystem_();
+        unsigned numCells = model_().numTotalDof();
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+        for (unsigned globI = 0; globI < numCells; globI++) {
+            const auto& neighbours = neighbours_[globI]; // this is a set but should maybe be changed
+            // accumulation term
+            double dt = simulator_().timeStepSize();
+            double volume = model_().dofTotalVolume(globI);
+            Scalar storefac = volume / dt;
+            ADVectorBlock adres(0.0);
+            const IntensiveQuantities* intQuantsInP = model_().cachedIntensiveQuantities(globI, /*timeIdx*/ 0);
+            assert(intQuantsInP);
+            const IntensiveQuantities& intQuantsIn = *intQuantsInP;
+            // intensiveQuantity(globI, 0);
+            LocalResidual::computeStorage(adres, intQuantsIn, 0);
+            adres *= storefac;
+            VectorBlock res(0.0);
+            MatrixBlock bMat(0.0);
+            setResAndJacobi(res, bMat, adres);
+            // first we use it as storage cache
+            if (model_().newtonMethod().numIterations() == 0) {
+                model_().updateCachedStorage(globI, /*timeIdx=*/1, res);
+            }
+            residual_[globI] -= model_().cachedStorage(globI, 1); //*storefac;
+            residual_[globI] += res;
+            jacobian_->addToBlock(globI, globI, bMat);
+            // wells sources for now (should be moved out)
+            if (well_local) {
+                res = 0.0;
+                bMat = 0.0;
+                adres = 0.0;
+                LocalResidual::computeSource(adres, problem_(), globI, 0);
+                adres *= -volume;
+                setResAndJacobi(res, bMat, adres);
+                residual_[globI] += res;
+                jacobian_->addToBlock(globI, globI, bMat);
+            }
+            short loc = 0;
+            for (const auto& globJ : neighbours) {
+                assert(globJ != globI);
+                res = 0.0;
+                bMat = 0.0;
+                adres = 0.0;
+                const IntensiveQuantities* intQuantsExP = model_().cachedIntensiveQuantities(globJ, /*timeIdx*/ 0);
+                assert(intQuantsExP);
+                const IntensiveQuantities& intQuantsEx = *intQuantsExP;
+                unsigned globalFocusDofIdx = globI;
+                LocalResidual::computeFlux(
+                    adres, problem_(), globalFocusDofIdx, globI, globJ, intQuantsIn, intQuantsEx, 0);
+                adres *= trans_[globI][loc];
+                setResAndJacobi(res, bMat, adres);
+                residual_[globI] += res;
+                jacobian_->addToBlock(globI, globI, bMat);
+                bMat *= -1.0;
+                jacobian_->addToBlock(globJ, globI, bMat);
+                loc++;
+            }
+        }
+        if (not(well_local)) {
+            problem_().wellModel().addReseroirSourceTerms(residual_, *jacobian_);
+        }
+        // before the first iteration of each time step, we need to update the
+        // constraints. (i.e., we assume that constraints can be time dependent, but they
+        // can't depend on the solution.)
+    }
+
     // linearize the whole system
     void linearize_()
     {
+        if (useTpfaLinearizer) {
+            linearizeGlobalTPFA_();
+            return;
+        }
+
         resetSystem_();
 
         // before the first iteration of each time step, we need to update the
@@ -595,6 +715,9 @@ private:
     LinearizationType linearizationType_;
 
     std::mutex globalMatrixMutex_;
+
+    SparseTable<unsigned> neighbours_;
+    SparseTable<double> trans_;
 };
 
 } // namespace Opm

opm/models/discretization/common/fvbaseproperties.hh

@@ -113,6 +113,9 @@ struct LocalLinearizer { using type = UndefinedProperty; };
 //! skipped
 template<class TypeTag, class MyTypeTag>
 struct LinearizeNonLocalElements { using type = UndefinedProperty; };
+//! Specify if the specialized TPFA linearizer should be used.
+template<class TypeTag, class MyTypeTag>
+struct UseTpfaLinearizer { using type = UndefinedProperty; };
 
 //! Linearizes the global non-linear system of equations
 template<class TypeTag, class MyTypeTag>