Merge pull request #703 from atgeirr/faster_assembly_minimal

Faster assembly minimal

examples/problems/fractureproblem.hh

@@ -50,7 +50,6 @@
 #include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
 #include <opm/material/components/SimpleH2O.hpp>
 #include <opm/material/components/Dnapl.hpp>
-#include <opm/material/common/Unused.hpp>
 
 #include <dune/common/version.hh>
 #include <dune/common/fmatrix.hh>

opm/models/blackoil/blackoilintensivequantities.hh

@@ -109,10 +109,20 @@ class BlackOilIntensiveQuantities
     using Toolbox = MathToolbox<Evaluation>;
     using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
     using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
-    using FluidState = BlackOilFluidState<Evaluation, FluidSystem, enableTemperature, enableEnergy, compositionSwitchEnabled,  enableEvaporation, enableBrine, enableSaltPrecipitation, Indices::numPhases >;
     using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>;
 
 public:
+    using FluidState = BlackOilFluidState<Evaluation,
+                                          FluidSystem,
+                                          enableTemperature,
+                                          enableEnergy,
+                                          compositionSwitchEnabled,
+                                          enableEvaporation,
+                                          enableBrine,
+                                          enableSaltPrecipitation,
+                                          Indices::numPhases>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+
     BlackOilIntensiveQuantities()
     {
         if (compositionSwitchEnabled) {
@@ -134,11 +144,18 @@ public:
 
         const auto& problem = elemCtx.problem();
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
-        const auto& linearizationType = elemCtx.linearizationType();
+
+        const auto& linearizationType = problem.model().linearizer().getLinearizationType();
+        unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
+        Scalar RvMax = FluidSystem::enableVaporizedOil()
+            ? problem.maxOilVaporizationFactor(timeIdx, globalSpaceIdx)
+            : 0.0;
+        Scalar RsMax = FluidSystem::enableDissolvedGas()
+            ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
+            : 0.0;
 
         asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
 
-        unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
         unsigned pvtRegionIdx = priVars.pvtRegionIndex();
         fluidState_.setPvtRegionIndex(pvtRegionIdx);
 
@@ -207,10 +224,10 @@ public:
 
         // now we compute all phase pressures
         Evaluation pC[numPhases];
-        const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
+        const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
         MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
 
-        //oil is the reference phase for pressure 
+        // oil is the reference phase for pressure
         if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv || priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
             const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
@@ -239,7 +256,7 @@ public:
         Evaluation SoMax = 0.0;
         if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
             SoMax = max(fluidState_.saturation(oilPhaseIdx),
-                        elemCtx.problem().maxOilSaturation(globalSpaceIdx));
+                        problem.maxOilSaturation(globalSpaceIdx));
         }
 
         // take the meaning of the switching primary variable into account for the gas
@@ -248,7 +265,6 @@ public:
             // in the threephase case, gas and oil phases are potentially present, i.e.,
             // we use the compositions of the gas-saturated oil and oil-saturated gas.
             if (FluidSystem::enableDissolvedGas()) {
-                Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
                 const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             oilPhaseIdx,
@@ -260,7 +276,6 @@ public:
                 fluidState_.setRs(0.0);
 
             if (FluidSystem::enableVaporizedOil()) {
-                Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
                 const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             gasPhaseIdx,
@@ -284,7 +299,6 @@ public:
             fluidState_.setRvw(Rvw);
 
             if (FluidSystem::enableVaporizedOil()) {
-                Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
                 const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             gasPhaseIdx,
@@ -306,7 +320,6 @@ public:
             if (FluidSystem::enableDissolvedGas()) {
                 // the oil phase is not present, but we need to compute its "composition" for
                 // the gravity correction anyway
-                Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
                 const auto& RsSat = enableExtbo ? asImp_().rs() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             oilPhaseIdx,
@@ -321,8 +334,6 @@ public:
         }
         else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
             // if the switching variable is the mole fraction of the gas component in the
-            Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
-
             // oil phase, we can directly set the composition of the oil phase
             const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
             fluidState_.setRs(min(RsMax, Rs));
@@ -330,7 +341,6 @@ public:
             if (FluidSystem::enableVaporizedOil()) {
                 // the gas phase is not present, but we need to compute its "composition"
                 // for the gravity correction anyway
-                Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
                 const auto& RvSat = enableExtbo ? asImp_().rv() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             gasPhaseIdx,
@@ -356,7 +366,6 @@ public:
             if (FluidSystem::enableDissolvedGas()) {
                 // the oil phase is not present, but we need to compute its "composition" for
                 // the gravity correction anyway
-                Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
                 const auto& RsSat = enableExtbo ? asImp_().rs() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             oilPhaseIdx,
@@ -380,7 +389,7 @@ public:
 
         typename FluidSystem::template ParameterCache<Evaluation> paramCache;
         paramCache.setRegionIndex(pvtRegionIdx);
-        if(FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
+        if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
             paramCache.setMaxOilSat(SoMax);
         }
         paramCache.updateAll(fluidState_);
@@ -447,9 +456,9 @@ public:
 
         // the porosity must be modified by the compressibility of the
         // rock...
-        Scalar rockCompressibility = problem.rockCompressibility(elemCtx, dofIdx, timeIdx);
+        Scalar rockCompressibility = problem.rockCompressibility(globalSpaceIdx);
         if (rockCompressibility > 0.0) {
-            Scalar rockRefPressure = problem.rockReferencePressure(elemCtx, dofIdx, timeIdx);
+            Scalar rockRefPressure = problem.rockReferencePressure(globalSpaceIdx);
             Evaluation x;
             if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
                 x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
@@ -477,6 +486,8 @@ public:
             porosity_ *= (1.0 - Sp);
         }
 
+        rockCompTransMultiplier_ = problem.template rockCompTransMultiplier<Evaluation>(*this, globalSpaceIdx);
+
         asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
         asImp_().zPvtUpdate_();
         asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
@@ -527,6 +538,12 @@ public:
     const Evaluation& porosity() const
     { return porosity_; }
 
+    /*!
+     * The pressure-dependent transmissibility multiplier due to rock compressibility.
+     */
+    const Evaluation& rockCompTransMultiplier() const
+    { return rockCompTransMultiplier_; }
+
     /*!
      * \brief Returns the index of the PVT region used to calculate the thermodynamic
      *        quantities.
@@ -578,6 +595,7 @@ private:
     FluidState fluidState_;
     Scalar referencePorosity_;
     Evaluation porosity_;
+    Evaluation rockCompTransMultiplier_;
     Evaluation mobility_[numPhases];
 };
 

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -0,0 +1,549 @@
+// -*- 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::BlackOilLocalResidual
+ */
+#ifndef EWOMS_BLACK_OIL_LOCAL_TPFA_RESIDUAL_HH
+#define EWOMS_BLACK_OIL_LOCAL_TPFA_RESIDUAL_HH
+
+#include "blackoilproperties.hh"
+#include "blackoilsolventmodules.hh"
+#include "blackoilextbomodules.hh"
+#include "blackoilpolymermodules.hh"
+#include "blackoilenergymodules.hh"
+#include "blackoilfoammodules.hh"
+#include "blackoilbrinemodules.hh"
+#include "blackoildiffusionmodule.hh"
+#include "blackoilmicpmodules.hh"
+#include <opm/material/fluidstates/BlackOilFluidState.hpp>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOilModel
+ *
+ * \brief Calculates the local residual of the black oil model.
+ */
+template <class TypeTag>
+class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLocalResidual>
+{
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
+    enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
+
+    enum { dimWorld = GridView::dimensionworld };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+
+    enum { gasCompIdx = FluidSystem::gasCompIdx };
+    enum { oilCompIdx = FluidSystem::oilCompIdx };
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+    enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
+
+    static const bool waterEnabled = Indices::waterEnabled;
+    static const bool gasEnabled = Indices::gasEnabled;
+    static const bool oilEnabled = Indices::oilEnabled;
+    static const bool compositionSwitchEnabled = (compositionSwitchIdx >= 0);
+
+    static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
+
+    static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();
+    static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
+    static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
+    static constexpr bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
+    static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
+    static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
+    static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
+    static constexpr bool enableMICP = getPropValue<TypeTag, Properties::EnableMICP>();
+
+    using SolventModule = BlackOilSolventModule<TypeTag>;
+    using ExtboModule = BlackOilExtboModule<TypeTag>;
+    using PolymerModule = BlackOilPolymerModule<TypeTag>;
+    using EnergyModule = BlackOilEnergyModule<TypeTag>;
+    using FoamModule = BlackOilFoamModule<TypeTag>;
+    using BrineModule = BlackOilBrineModule<TypeTag>;
+    using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+    using MICPModule = BlackOilMICPModule<TypeTag>;
+
+    using Toolbox = MathToolbox<Evaluation>;
+
+public:
+    /*!
+     * \copydoc FvBaseLocalResidual::computeStorage
+     */
+    template <class LhsEval>
+    void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
+                        const ElementContext& elemCtx,
+                        unsigned dofIdx,
+                        unsigned timeIdx) const
+    {
+        const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
+        computeStorage(storage,
+                       intQuants);
+    }
+
+    template <class LhsEval>
+    static void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
+                               const IntensiveQuantities& intQuants)
+    {
+        // retrieve the intensive quantities for the SCV at the specified point in time
+        const auto& fs = intQuants.fluidState();
+        storage = 0.0;
+
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
+            LhsEval surfaceVolume =
+                Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
+                * Toolbox::template decay<LhsEval>(fs.invB(phaseIdx))
+                * Toolbox::template decay<LhsEval>(intQuants.porosity());
+
+            storage[conti0EqIdx + activeCompIdx] += surfaceVolume;
+
+            // account for dissolved gas
+            if (phaseIdx == oilPhaseIdx && FluidSystem::enableDissolvedGas()) {
+                unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
+                storage[conti0EqIdx + activeGasCompIdx] +=
+                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs())
+                    * surfaceVolume;
+            }
+
+            // account for vaporized oil
+            if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedOil()) {
+                unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
+                storage[conti0EqIdx + activeOilCompIdx] +=
+                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rv())
+                    * surfaceVolume;
+            }
+
+            // account for vaporized water
+            if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedWater()) {
+                unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
+                storage[conti0EqIdx + activeWaterCompIdx] +=
+                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rvw())
+                    * surfaceVolume;
+            }
+        }
+
+        adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex());
+
+        // deal with solvents (if present)
+        SolventModule::addStorage(storage, intQuants);
+
+        // deal with zFracton (if present)
+        ExtboModule::addStorage(storage, intQuants);
+
+        // deal with polymer (if present)
+        PolymerModule::addStorage(storage, intQuants);
+
+        // deal with energy (if present)
+        EnergyModule::addStorage(storage, intQuants);
+
+        // deal with foam (if present)
+        FoamModule::addStorage(storage, intQuants);
+
+        // deal with salt (if present)
+        BrineModule::addStorage(storage, intQuants);
+
+        // deal with micp (if present)
+        MICPModule::addStorage(storage, intQuants);
+    }
+
+    /*!
+     * This function works like the ElementContext-based version with
+     * one main difference: The darcy flux is calculated here, not
+     * read from the extensive quantities of the element context.
+     */
+    static void computeFlux(RateVector& flux,
+                            const Problem& problem,
+                            const unsigned globalIndexIn,
+                            const unsigned globalIndexEx,
+                            const IntensiveQuantities& intQuantsIn,
+                            const IntensiveQuantities& intQuantsEx,
+                            const Scalar trans,
+                            const Scalar faceArea)
+    {
+        flux = 0.0;
+        Scalar Vin = problem.model().dofTotalVolume(globalIndexIn);
+        Scalar Vex = problem.model().dofTotalVolume(globalIndexEx);
+
+        Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx);
+
+        // estimate the gravity correction: for performance reasons we use a simplified
+        // approach for this flux module that assumes that gravity is constant and always
+        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
+        Scalar g = problem.gravity()[dimWorld - 1];
+
+        // this is quite hacky because the dune grid interface does not provide a
+        // cellCenterDepth() method (so we ask the problem to provide it). The "good"
+        // solution would be to take the Z coordinate of the element centroids, but since
+        // ECL seems to like to be inconsistent on that front, it needs to be done like
+        // here...
+        Scalar zIn = problem.dofCenterDepth(globalIndexIn);
+        Scalar zEx = problem.dofCenterDepth(globalIndexEx);
+
+        // the distances from the DOF's depths. (i.e., the additional depth of the
+        // exterior DOF)
+        Scalar distZ = zIn - zEx; // NB could be precalculated
+
+        calculateFluxes_(flux,
+                         intQuantsIn,
+                         intQuantsEx,
+                         Vin,
+                         Vex,
+                         globalIndexIn,
+                         globalIndexEx,
+                         distZ * g,
+                         thpres,
+                         trans,
+                         faceArea);
+    }
+
+    // This function demonstrates compatibility with the ElementContext-based interface.
+    // Actually using it will lead to double work since the element context already contains
+    // fluxes through its stored ExtensiveQuantities.
+    static void computeFlux(RateVector& flux,
+                            const ElementContext& elemCtx,
+                            unsigned scvfIdx,
+                            unsigned timeIdx)
+    {
+        assert(timeIdx == 0);
+
+        flux = 0.0;
+        // need for dary flux calculation
+        const auto& problem = elemCtx.problem();
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& scvf = stencil.interiorFace(scvfIdx);
+
+        unsigned interiorDofIdx = scvf.interiorIndex();
+        unsigned exteriorDofIdx = scvf.exteriorIndex();
+        assert(interiorDofIdx != exteriorDofIdx);
+
+        // unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
+        // unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
+        Scalar Vin = elemCtx.dofVolume(interiorDofIdx, /*timeIdx=*/0);
+        Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, /*timeIdx=*/0);
+        const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx);
+        const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx);
+        Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
+        Scalar faceArea = scvf.area();
+        Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx);
+
+        // estimate the gravity correction: for performance reasons we use a simplified
+        // approach for this flux module that assumes that gravity is constant and always
+        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
+        Scalar g = problem.gravity()[dimWorld - 1];
+        const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
+        const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
+
+        // this is quite hacky because the dune grid interface does not provide a
+        // cellCenterDepth() method (so we ask the problem to provide it). The "good"
+        // solution would be to take the Z coordinate of the element centroids, but since
+        // ECL seems to like to be inconsistent on that front, it needs to be done like
+        // here...
+        Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
+        Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
+
+        // the distances from the DOF's depths. (i.e., the additional depth of the
+        // exterior DOF)
+        Scalar distZ = zIn - zEx;
+
+        calculateFluxes_(flux,
+                         intQuantsIn,
+                         intQuantsEx,
+                         Vin,
+                         Vex,
+                         globalIndexIn,
+                         globalIndexEx,
+                         distZ * g,
+                         thpres,
+                         trans,
+                         faceArea);
+    }
+
+    static void calculateFluxes_(RateVector& flux,
+                                 const IntensiveQuantities& intQuantsIn,
+                                 const IntensiveQuantities& intQuantsEx,
+                                 const Scalar& Vin,
+                                 const Scalar& Vex,
+                                 const unsigned& globalIndexIn,
+                                 const unsigned& globalIndexEx,
+                                 const Scalar& distZg,
+                                 const Scalar& thpres,
+                                 const Scalar& trans,
+                                 const Scalar& faceArea)
+    {
+        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            if (!FluidSystem::phaseIsActive(phaseIdx))
+                continue;
+            // darcy flux calculation
+            short dnIdx;
+            //
+            short upIdx;
+            // fake intices should only be used to get upwind anc compatibility with old functions
+            short interiorDofIdx = 0; // NB
+            short exteriorDofIdx = 1; // NB
+            Evaluation pressureDifference;
+            ExtensiveQuantities::calculatePhasePressureDiff_(upIdx,
+                                                             dnIdx,
+                                                             pressureDifference,
+                                                             intQuantsIn,
+                                                             intQuantsEx,
+                                                             phaseIdx, // input
+                                                             interiorDofIdx, // input
+                                                             exteriorDofIdx, // intput
+                                                             Vin,
+                                                             Vex,
+                                                             globalIndexIn,
+                                                             globalIndexEx,
+                                                             distZg,
+                                                             thpres);
+
+
+
+            const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
+            unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
+            const Evaluation& transMult = up.rockCompTransMultiplier();
+            Evaluation darcyFlux;
+            if (pressureDifference == 0) {
+                darcyFlux = 0.0; // NB maybe we could drop calculations
+            } else {
+                if (globalUpIndex == globalIndexIn)
+                    darcyFlux = pressureDifference * up.mobility(phaseIdx) * transMult * (-trans / faceArea);
+                else
+                    darcyFlux = pressureDifference * (Toolbox::value(up.mobility(phaseIdx)) * Toolbox::value(transMult) * (-trans / faceArea));
+            }
+
+            unsigned pvtRegionIdx = up.pvtRegionIndex();
+            using FluidState = typename IntensiveQuantities::FluidState;
+            // if (upIdx == globalFocusDofIdx){
+            if (globalUpIndex == globalIndexIn) {
+                const auto& invB
+                    = getInvB_<FluidSystem, FluidState, Evaluation>(up.fluidState(), phaseIdx, pvtRegionIdx);
+                const auto& surfaceVolumeFlux = invB * darcyFlux;
+                evalPhaseFluxes_<Evaluation, Evaluation, FluidState>(
+                    flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
+            } else {
+                const auto& invB = getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), phaseIdx, pvtRegionIdx);
+                const auto& surfaceVolumeFlux = invB * darcyFlux;
+                evalPhaseFluxes_<Scalar, Evaluation, FluidState>(
+                    flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState());
+            }
+        }
+
+        // deal with solvents (if present)
+        static_assert(!enableSolvent, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with zFracton (if present)
+        static_assert(!enableExtbo, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with polymer (if present)
+        static_assert(!enablePolymer, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with energy (if present)
+        static_assert(!enableEnergy, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with foam (if present)
+        static_assert(!enableFoam, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with salt (if present)
+        static_assert(!enableBrine, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with diffusion (if present)
+        static_assert(!enableDiffusion, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with micp (if present)
+        static_assert(!enableMICP, "Relevant computeFlux() method must be implemented for this module before enabling.");
+        // MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+    }
+
+    static void computeSource(RateVector& source,
+                              const Problem& problem,
+                              unsigned globalSpaceIdex,
+                              unsigned timeIdx)
+    {
+        // retrieve the source term intrinsic to the problem
+        problem.source(source, globalSpaceIdex, timeIdx);
+
+        // deal with MICP (if present)
+        // deal with micp (if present)
+        static_assert(!enableMICP, "Relevant addSource() method must be implemented for this module before enabling.");
+        // MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
+
+        // scale the source term of the energy equation
+        if (enableEnergy)
+            source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+    }
+
+    /*!
+     * \copydoc FvBaseLocalResidual::computeSource
+     */
+    void computeSource(RateVector& source,
+                       const ElementContext& elemCtx,
+                       unsigned dofIdx,
+                       unsigned timeIdx) const
+    {
+        // retrieve the source term intrinsic to the problem
+        elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
+
+        // deal with MICP (if present)
+        MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
+
+        // scale the source term of the energy equation
+        if (enableEnergy)
+            source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+    }
+
+    template <class UpEval, class FluidState>
+    static void evalPhaseFluxes_(RateVector& flux,
+                                 unsigned phaseIdx,
+                                 unsigned pvtRegionIdx,
+                                 const ExtensiveQuantities& extQuants,
+                                 const FluidState& upFs)
+    {
+
+        const auto& invB = getInvB_<FluidSystem, FluidState, UpEval>(upFs, phaseIdx, pvtRegionIdx);
+        const auto& surfaceVolumeFlux = invB * extQuants.volumeFlux(phaseIdx);
+        evalPhaseFluxes_<UpEval>(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, upFs);
+    }
+
+    /*!
+     * \brief Helper function to calculate the flux of mass in terms of conservation
+     *        quantities via specific fluid phase over a face.
+     */
+    template <class UpEval, class Eval,class FluidState>
+    static void evalPhaseFluxes_(RateVector& flux,
+                                 unsigned phaseIdx,
+                                 unsigned pvtRegionIdx,
+                                 const Eval& surfaceVolumeFlux,
+                                 const FluidState& upFs)
+    {
+        unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx));
+
+        if (blackoilConserveSurfaceVolume)
+            flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux;
+        else
+            flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux*FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
+
+        if (phaseIdx == oilPhaseIdx) {
+            // dissolved gas (in the oil phase).
+            if (FluidSystem::enableDissolvedGas()) {
+                const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
+
+                unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
+                if (blackoilConserveSurfaceVolume)
+                    flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux;
+                else
+                    flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+            }
+        }
+        else if (phaseIdx == gasPhaseIdx) {
+            // vaporized oil (in the gas phase).
+            if (FluidSystem::enableVaporizedOil()) {
+                const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
+
+                unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
+                if (blackoilConserveSurfaceVolume)
+                    flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux;
+                else
+                    flux[conti0EqIdx + activeOilCompIdx] += Rv*surfaceVolumeFlux*FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+            }
+             // vaporized water (in the gas phase).
+            if (FluidSystem::enableVaporizedWater()) {
+                const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);
+
+                unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
+                if (blackoilConserveSurfaceVolume)
+                    flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux;
+                else
+                    flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux*FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+            }
+        }
+    }
+
+    /*!
+     * \brief Helper function to convert the mass-related parts of a Dune::FieldVector
+     *        that stores conservation quantities in terms of "surface-volume" to the
+     *        conservation quantities used by the model.
+     *
+     * Depending on the value of the BlackoilConserveSurfaceVolume property, the model
+     * either conserves mass by means of "surface volume" of the components or mass
+     * directly. In the former case, this method is a no-op; in the latter, the values
+     * passed are multiplied by their respective pure component's density at surface
+     * conditions.
+     */
+    template <class Scalar>
+    static void adaptMassConservationQuantities_(Dune::FieldVector<Scalar, numEq>& container, unsigned pvtRegionIdx)
+    {
+        if (blackoilConserveSurfaceVolume)
+            return;
+
+        // convert "surface volume" to mass. this is complicated a bit by the fact that
+        // not all phases are necessarily enabled. (we here assume that if a fluid phase
+        // is disabled, its respective "main" component is not considered as well.)
+
+        if (waterEnabled) {
+            unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx);
+            container[conti0EqIdx + activeWaterCompIdx] *=
+                FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+        }
+
+        if (gasEnabled) {
+            unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx);
+            container[conti0EqIdx + activeGasCompIdx] *=
+                FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+        }
+
+        if (oilEnabled) {
+            unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx);
+            container[conti0EqIdx + activeOilCompIdx] *=
+                FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
+        }
+    }
+};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilmodel.hh

@@ -298,6 +298,9 @@ class BlackOilModel
     using MICPModule = BlackOilMICPModule<TypeTag>;
 
 public:
+
+    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
+
     BlackOilModel(Simulator& simulator)
         : ParentType(simulator)
     {}

opm/models/blackoil/blackoilproblem.hh

@@ -130,6 +130,12 @@ public:
                                unsigned) const
     { return 0.0; }
 
+    /*!
+     * \brief Returns the compressibility of the porous medium of a cell
+     */
+    Scalar rockCompressibility(unsigned) const
+    { return 0.0; }
+
     /*!
      * \brief Returns the reference pressure for rock the compressibility of a cell
      */
@@ -139,6 +145,12 @@ public:
                                  unsigned) const
     { return 1e5; }
 
+    /*!
+     * \brief Returns the reference pressure for rock the compressibility of a cell
+     */
+    Scalar rockReferencePressure(unsigned) const
+    { return 1e5; }
+
     /*!
      * \brief Returns the reference temperature
      *
@@ -159,6 +171,11 @@ public:
                                   unsigned) const
     { return 1.0; }
 
+    template <class LhsEval>
+    LhsEval rockCompTransMultiplier(const IntensiveQuantities&,
+                                    unsigned) const
+    { return 1.0; }
+
 private:
     //! Returns the implementation of the problem (i.e. static polymorphism)
     Implementation& asImp_()

opm/models/discretefracture/discretefractureproblem.hh

@@ -33,7 +33,6 @@
 #include <opm/models/common/multiphasebaseproblem.hh>
 
 #include <opm/material/common/Means.hpp>
-#include <opm/material/common/Unused.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>

opm/models/discretization/common/fvbasediscretization.hh

@@ -1411,7 +1411,7 @@ public:
         // previous time step so that we can start the next
         // update at a physically meaningful solution.
         solution(/*timeIdx=*/0) = solution(/*timeIdx=*/1);
-        invalidateIntensiveQuantitiesCache(/*timeIdx=*/0);
+        invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
 
 #ifndef NDEBUG
         for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {

opm/models/discretization/common/tpfalinearizer.hh

@@ -0,0 +1,476 @@
+// -*- 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 TPFA_LINEARIZER_HH
+#define TPFA_LINEARIZER_HH
+
+#include "fvbaseproperties.hh"
+#include "linearizationtype.hh"
+
+#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>
+#include <dune/common/fmatrix.hh>
+
+#include <type_traits>
+#include <iostream>
+#include <vector>
+#include <thread>
+#include <set>
+#include <exception>   // current_exception, rethrow_exception
+#include <mutex>
+
+
+namespace Opm {
+// forward declarations
+template<class TypeTag>
+class EcfvDiscretization;
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief The common code for the linearizers of non-linear systems of equations
+ *
+ * This class assumes that these system of equations to be linearized are stemming from
+ * models that use an finite volume scheme for spatial discretization and an Euler
+ * scheme for time discretization.
+ */
+template<class TypeTag>
+class TpfaLinearizer
+{
+//! \cond SKIP_THIS
+    using Model = GetPropType<TypeTag, Properties::Model>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+
+    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+    using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
+    using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
+    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
+    using Constraints = GetPropType<TypeTag, Properties::Constraints>;
+    using Stencil = GetPropType<TypeTag, Properties::Stencil>;
+    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+
+    using Element = typename GridView::template Codim<0>::Entity;
+    using ElementIterator = typename GridView::template Codim<0>::Iterator;
+
+    using Vector = GlobalEqVector;
+
+    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+    enum { historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>() };
+
+    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>();
+
+    // copying the linearizer is not a good idea
+    TpfaLinearizer(const TpfaLinearizer&);
+//! \endcond
+
+public:
+    TpfaLinearizer()
+        : jacobian_()
+    {
+        simulatorPtr_ = 0;
+    }
+
+    ~TpfaLinearizer()
+    {
+    }
+
+    /*!
+     * \brief Register all run-time parameters for the Jacobian linearizer.
+     */
+    static void registerParameters()
+    { }
+
+    /*!
+     * \brief Initialize the linearizer.
+     *
+     * At this point we can assume that all objects in the simulator
+     * have been allocated. We cannot assume that they are fully
+     * initialized, though.
+     *
+     * \copydetails Doxygen::simulatorParam
+     */
+    void init(Simulator& simulator)
+    {
+        simulatorPtr_ = &simulator;
+        eraseMatrix();
+    }
+
+    /*!
+     * \brief Causes the Jacobian matrix to be recreated from scratch before the next
+     *        iteration.
+     *
+     * This method is usally called if the sparsity pattern has changed for some
+     * reason. (e.g. by modifications of the grid or changes of the auxiliary equations.)
+     */
+    void eraseMatrix()
+    {
+        jacobian_.reset();
+    }
+
+    /*!
+     * \brief Linearize the full system of non-linear 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()
+    {
+        linearizeDomain();
+        linearizeAuxiliaryEquations();
+    }
+
+    /*!
+     * \brief Linearize the part of the non-linear system of equations that is associated
+     *        with the spatial domain.
+     *
+     * That means that the global Jacobian of the residual is assembled and the residual
+     * is evaluated for the current solution.
+     *
+     * The current state of affairs (esp. the previous and the current solutions) is
+     * represented by the model object.
+     */
+    void linearizeDomain()
+    {
+        // we defer the initialization of the Jacobian matrix until here because the
+        // auxiliary modules usually assume the problem, model and grid to be fully
+        // initialized...
+        if (!jacobian_)
+            initFirstIteration_();
+
+        int succeeded;
+        try {
+            linearize_();
+            succeeded = 1;
+        }
+        catch (const std::exception& e)
+        {
+            std::cout << "rank " << simulator_().gridView().comm().rank()
+                      << " caught an exception while linearizing:" << e.what()
+                      << "\n"  << std::flush;
+            succeeded = 0;
+        }
+        catch (...)
+        {
+            std::cout << "rank " << simulator_().gridView().comm().rank()
+                      << " caught an exception while linearizing"
+                      << "\n"  << std::flush;
+            succeeded = 0;
+        }
+        succeeded = gridView_().comm().min(succeeded);
+
+        if (!succeeded)
+            throw NumericalIssue("A process did not succeed in linearizing the system");
+    }
+
+    void finalize()
+    { jacobian_->finalize(); }
+
+    /*!
+     * \brief Linearize the part of the non-linear system of equations that is associated
+     *        with the spatial domain.
+     */
+    void linearizeAuxiliaryEquations()
+    {
+        // flush possible local caches into matrix structure
+        jacobian_->commit();
+
+        auto& model = model_();
+        const auto& comm = simulator_().gridView().comm();
+        for (unsigned auxModIdx = 0; auxModIdx < model.numAuxiliaryModules(); ++auxModIdx) {
+            bool succeeded = true;
+            try {
+                model.auxiliaryModule(auxModIdx)->linearize(*jacobian_, residual_);
+            }
+            catch (const std::exception& e) {
+                succeeded = false;
+
+                std::cout << "rank " << simulator_().gridView().comm().rank()
+                          << " caught an exception while linearizing:" << e.what()
+                          << "\n"  << std::flush;
+            }
+
+            succeeded = comm.min(succeeded);
+
+            if (!succeeded)
+                throw NumericalIssue("linearization of an auxiliary equation failed");
+        }
+    }
+
+    /*!
+     * \brief Return constant reference to global Jacobian matrix backend.
+     */
+    const SparseMatrixAdapter& jacobian() const
+    { return *jacobian_; }
+
+    SparseMatrixAdapter& jacobian()
+    { return *jacobian_; }
+
+    /*!
+     * \brief Return constant reference to global residual vector.
+     */
+    const GlobalEqVector& residual() const
+    { return residual_; }
+
+    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.
+     *
+     * (This object is only non-empty if the EnableConstraints property is true.)
+     */
+    const std::map<unsigned, Constraints> constraintsMap() const
+    { return {}; }
+
+private:
+    Simulator& simulator_()
+    { return *simulatorPtr_; }
+    const Simulator& simulator_() const
+    { return *simulatorPtr_; }
+
+    Problem& problem_()
+    { return simulator_().problem(); }
+    const Problem& problem_() const
+    { return simulator_().problem(); }
+
+    Model& model_()
+    { return simulator_().model(); }
+    const Model& model_() const
+    { return simulator_().model(); }
+
+    const GridView& gridView_() const
+    { return problem_().gridView(); }
+
+    void initFirstIteration_()
+    {
+        // initialize the BCRS matrix for the Jacobian of the residual function
+        createMatrix_();
+
+        // initialize the Jacobian matrix and the vector for the residual function
+        residual_.resize(model_().numTotalDof());
+        resetSystem_();
+    }
+
+    // Construct the BCRS matrix for the Jacobian of the residual function
+    void createMatrix_()
+    {
+        const auto& model = model_();
+        Stencil stencil(gridView_(), model_().dofMapper());
+
+        // for the main model, find out the global indices of the neighboring degrees of
+        // freedom of each primary degree of freedom
+        using NeighborSet = std::set< unsigned >;
+        std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
+
+        unsigned numCells = model.numTotalDof();
+        neighborInfo_.reserve(numCells, 6 * numCells);
+        std::vector<NeighborInfo> loc_nbinfo;
+
+        ElementIterator elemIt = gridView_().template begin<0>();
+        const ElementIterator elemEndIt = gridView_().template end<0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const Element& elem = *elemIt;
+            stencil.update(elem);
+
+            for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx) {
+                unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
+                loc_nbinfo.resize(stencil.numDof() - 1); // Do not include the primary dof in neighborInfo_
+
+                for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
+                    unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
+                    sparsityPattern[myIdx].insert(neighborIdx);
+                    if (dofIdx > 0) {
+                        const double trans = problem_().transmissibility(myIdx, neighborIdx);
+                        const double area = stencil.interiorFace(dofIdx - 1).area();
+                        loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area};
+                    }
+                }
+                neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
+            }
+        }
+
+        // add the additional neighbors and degrees of freedom caused by the auxiliary
+        // equations
+        size_t numAuxMod = model.numAuxiliaryModules();
+        for (unsigned auxModIdx = 0; auxModIdx < numAuxMod; ++auxModIdx)
+            model.auxiliaryModule(auxModIdx)->addNeighbors(sparsityPattern);
+
+        // allocate raw matrix
+        jacobian_.reset(new SparseMatrixAdapter(simulator_()));
+
+        // create matrix structure based on sparsity pattern
+        jacobian_->reserve(sparsityPattern);
+    }
+
+    // reset the global linear system of equations.
+    void resetSystem_()
+    {
+        residual_ = 0.0;
+        // zero all matrix entries
+        jacobian_->clear();
+    }
+
+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 linearize_()
+    {
+        const bool well_local = true;
+        resetSystem_();
+        unsigned numCells = model_().numTotalDof();
+#ifdef _OPENMP
+#pragma omp parallel for
+#endif
+        for (unsigned globI = 0; globI < numCells; globI++) {
+            const auto& nbInfos = neighborInfo_[globI]; // this is a set but should maybe be changed
+            VectorBlock res(0.0);
+            MatrixBlock bMat(0.0);
+            ADVectorBlock adres(0.0);
+            const IntensiveQuantities* intQuantsInP = model_().cachedIntensiveQuantities(globI, /*timeIdx*/ 0);
+            if (intQuantsInP == nullptr) {
+                throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globI));
+            }
+            const IntensiveQuantities& intQuantsIn = *intQuantsInP;
+
+            // Flux term.
+            short loc = 0;
+            for (const auto& nbInfo : nbInfos) {
+                unsigned globJ = nbInfo.neighbor;
+                assert(globJ != globI);
+                res = 0.0;
+                bMat = 0.0;
+                adres = 0.0;
+                const IntensiveQuantities* intQuantsExP = model_().cachedIntensiveQuantities(globJ, /*timeIdx*/ 0);
+                if (intQuantsExP == nullptr) {
+                    throw std::logic_error("Missing updated intensive quantities for cell " + std::to_string(globJ) + " when assembling fluxes for cell " + std::to_string(globI));
+                }
+                const IntensiveQuantities& intQuantsEx = *intQuantsExP;
+                LocalResidual::computeFlux(
+                    adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, nbInfo.trans, nbInfo.faceArea);
+                adres *= nbInfo.faceArea;
+                setResAndJacobi(res, bMat, adres);
+                residual_[globI] += res;
+                jacobian_->addToBlock(globI, globI, bMat);
+                bMat *= -1.0;
+                jacobian_->addToBlock(globJ, globI, bMat);
+                ++loc;
+            }
+
+            // Accumulation term.
+            double dt = simulator_().timeStepSize();
+            double volume = model_().dofTotalVolume(globI);
+            Scalar storefac = volume / dt;
+            adres = 0.0;
+            LocalResidual::computeStorage(adres, intQuantsIn);
+            setResAndJacobi(res, bMat, adres);
+            // TODO: check recycleFirst etc.
+            // first we use it as storage cache
+            if (model_().newtonMethod().numIterations() == 0) {
+                model_().updateCachedStorage(globI, /*timeIdx=*/1, res);
+            }
+            res -= model_().cachedStorage(globI, 1);
+            res *= storefac;
+            bMat *= storefac;
+            // 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);
+            }
+        }
+    }
+
+    Simulator *simulatorPtr_;
+
+    // the jacobian matrix
+    std::unique_ptr<SparseMatrixAdapter> jacobian_;
+
+    // the right-hand side
+    GlobalEqVector residual_;
+
+    LinearizationType linearizationType_;
+
+    struct NeighborInfo
+    {
+        unsigned int neighbor;
+        double trans;
+        double faceArea;
+    };
+    SparseTable<NeighborInfo> neighborInfo_;
+};
+
+} // namespace Opm
+
+#endif // TPFA_LINEARIZER