working with small increase in performance

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -57,13 +57,14 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     using EqVector = GetPropType<TypeTag, Properties::EqVector>;
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
-    
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
 
     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 };
@@ -195,27 +196,117 @@ public:
         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);
 
-        const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        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 = elemCtx.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;
+
+        //
+        //const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
         unsigned focusDofIdx = elemCtx.focusDofIndex();
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx))
                 continue;
+            // darcy flux calculation
+            short dnIdx;
+            //
+            short upIdx;
+            Evaluation pressureDifference;
+            ExtensiveQuantities::calculatePhasePressureDiff_(upIdx,
+                                                             dnIdx,
+                                                             pressureDifference,
+                                                             intQuantsIn,
+                                                             intQuantsEx,
+                                                             scvfIdx,//input
+                                                             timeIdx,//input
+                                                             phaseIdx,//input
+                                                             interiorDofIdx,//input
+                                                             exteriorDofIdx,//intput
+                                                             Vin,
+                                                             Vex,
+                                                             globalIndexIn,
+                                                             globalIndexEx,
+                                                             distZ*g,
+                                                             thpres);
+            
+            
+            
+            const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
+                unsigned globalIndex;
+                if(upIdx == interiorDofIdx){
+                    //up = intQuantsIn;
+                    globalIndex = globalIndexIn;
+                }else{
+                    //up = intQuantsEx;
+                    globalIndex = globalIndexEx;
+                }
+                // TODO: should the rock compaction transmissibility multiplier be upstreamed
+                // or averaged? all fluids should see the same compaction?!
+                //const auto& globalIndex = stencil.globalSpaceIndex(upstreamIdx);
+                const Evaluation& transMult =
+                    problem.template rockCompTransMultiplier<Evaluation>(up, globalIndex);
+                Evaluation darcyFlux;
+                if(pressureDifference == 0){
+                    darcyFlux = 0.0; //NB maybe we could drop calculations
+                }else{
+                    if (upIdx == interiorDofIdx)
+                        darcyFlux =
+                            pressureDifference*up.mobility(phaseIdx)*transMult*(-trans/faceArea);
+                    else
+                        darcyFlux =
+                            pressureDifference*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea));
+                }
+                //const auto& darcyFlux = extQuants.volumeFlux(phaseIdx);                       
+                //unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
+                               
+                //const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
+                unsigned pvtRegionIdx = up.pvtRegionIndex();
+                using FluidState = typename IntensiveQuantities::FluidState;
+                if (upIdx == focusDofIdx){
+                    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());
+                }
+            
             
-            const auto& darcyFlux = extQuants.volumeFlux(phaseIdx);                       
-            unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
-            const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
-            unsigned pvtRegionIdx = up.pvtRegionIndex();
-            using FluidState = typename IntensiveQuantities::FluidState;
-            if (upIdx == focusDofIdx){
-                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)

opm/models/discretization/common/fvbaseadlocallinearizertpfa.hh

@@ -0,0 +1,308 @@
+// -*- 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::FvBaseAdLocalLinearizer
+ */
+#ifndef EWOMS_FV_BASE_AD_LOCAL_TPFA_LINEARIZER_HH
+#define EWOMS_FV_BASE_AD_LOCAL_TPFA_LINEARIZER_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/material/densead/Math.hpp>
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/matrix.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+namespace Opm {
+//forward declaration
+template<class TypeTag>
+class FvBaseAdLocalLinearizerTPFA;
+}
+
+namespace Opm::Properties {
+
+//declare the property tags required for the finite differences local linearizer
+
+namespace TTag {
+struct AutoDiffLocalLinearizerTPFA {};
+} // namespace TTag
+
+// set the properties to be spliced in
+template<class TypeTag>
+struct LocalLinearizer<TypeTag, TTag::AutoDiffLocalLinearizerTPFA>
+{ using type = FvBaseAdLocalLinearizerTPFA<TypeTag>; };
+
+//! Set the function evaluation w.r.t. the primary variables
+template<class TypeTag>
+struct Evaluation<TypeTag, TTag::AutoDiffLocalLinearizerTPFA>
+{
+private:
+    static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
+public:
+    using type = DenseAd::Evaluation<Scalar, numEq>;
+};
+
+} // namespace Opm::Properties
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Calculates the local residual and its Jacobian for a single element of the grid.
+ *
+ * This class uses automatic differentiation to calculate the partial derivatives (the
+ * alternative is finite differences).
+ */
+template<class TypeTag>
+class FvBaseAdLocalLinearizerTPFA
+{
+private:
+    using Implementation = GetPropType<TypeTag, Properties::LocalLinearizer>;
+    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using Model = GetPropType<TypeTag, Properties::Model>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+
+    using ScalarVectorBlock = Dune::FieldVector<Scalar, numEq>;
+    // extract local matrices from jacobian matrix for consistency
+    using ScalarMatrixBlock = typename GetPropType<TypeTag, Properties::SparseMatrixAdapter>::MatrixBlock;
+
+    using ScalarLocalBlockVector = Dune::BlockVector<ScalarVectorBlock>;
+    using ScalarLocalBlockMatrix = Dune::Matrix<ScalarMatrixBlock>;
+
+public:
+    FvBaseAdLocalLinearizerTPFA()
+        : internalElemContext_(0)
+    { }
+
+    // copying local linearizer objects around is a very bad idea, so we explicitly
+    // prevent it...
+    FvBaseAdLocalLinearizerTPFA(const FvBaseAdLocalLinearizerTPFA&) = delete;
+
+    ~FvBaseAdLocalLinearizerTPFA()
+    { delete internalElemContext_; }
+
+    /*!
+     * \brief Register all run-time parameters for the local jacobian.
+     */
+    static void registerParameters()
+    { }
+
+    /*!
+     * \brief Initialize the local Jacobian object.
+     *
+     * At this point we can assume that everything has been allocated,
+     * although some objects may not yet be completely initialized.
+     *
+     * \param simulator The simulator object of the simulation.
+     */
+    void init(Simulator& simulator)
+    {
+        simulatorPtr_ = &simulator;
+        delete internalElemContext_;
+        internalElemContext_ = new ElementContext(simulator);
+    }
+
+    /*!
+     * \brief Compute an element's local Jacobian matrix and evaluate its residual.
+     *
+     * The local Jacobian for a given context is defined as the derivatives of the
+     * residuals of all degrees of freedom featured by the stencil with regard to the
+     * primary variables of the stencil's "primary" degrees of freedom. Adding the local
+     * Jacobians for all elements in the grid will give the global Jacobian 'grad f(x)'.
+     *
+     * \param element The grid element for which the local residual and its local
+     *                Jacobian should be calculated.
+     */
+    void linearize(const Element& element)
+    {
+        linearize(*internalElemContext_, element);
+    }
+
+    /*!
+     * \brief Compute an element's local Jacobian matrix and evaluate its residual.
+     *
+     * The local Jacobian for a given context is defined as the derivatives of the
+     * residuals of all degrees of freedom featured by the stencil with regard to the
+     * primary variables of the stencil's "primary" degrees of freedom. Adding the local
+     * Jacobians for all elements in the grid will give the global Jacobian 'grad f(x)'.
+     *
+     * After calling this method the ElementContext is in an undefined state, so do not
+     * use it anymore!
+     *
+     * \param elemCtx The element execution context for which the local residual and its
+     *                local Jacobian should be calculated.
+     */
+    void linearize(ElementContext& elemCtx, const Element& elem)
+    {
+        elemCtx.updateStencil(elem);
+        elemCtx.updateAllIntensiveQuantities();
+
+        // update the weights of the primary variables for the context
+        model_().updatePVWeights(elemCtx);
+
+        // resize the internal arrays of the linearizer
+        resize_(elemCtx);
+
+        // compute the local residual and its Jacobian
+        unsigned numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+        for (unsigned focusDofIdx = 0; focusDofIdx < numPrimaryDof; focusDofIdx++) {
+            elemCtx.setFocusDofIndex(focusDofIdx);
+            //elemCtx.updateAllExtensiveQuantities();//NB should not be need anymore
+
+            // calculate the local residual
+            localResidual_.eval(elemCtx);
+
+            // convert the local Jacobian matrix and the right hand side from the data
+            // structures used by the automatic differentiation code to the conventional
+            // ones used by the linear solver.
+            updateLocalLinearization_(elemCtx, focusDofIdx);
+        }
+    }
+
+    /*!
+     * \brief Return reference to the local residual.
+     */
+    LocalResidual& localResidual()
+    { return localResidual_; }
+
+    /*!
+     * \brief Return reference to the local residual.
+     */
+    const LocalResidual& localResidual() const
+    { return localResidual_; }
+
+    /*!
+     * \brief Returns the local Jacobian matrix of the residual of a sub-control volume.
+     *
+     * \param domainScvIdx The local index of the sub control volume to which the primary
+     *                     variables are associated with
+     * \param rangeScvIdx The local index of the sub control volume which contains the
+     *                    local residual
+     */
+    const ScalarMatrixBlock& jacobian(unsigned domainScvIdx, unsigned rangeScvIdx) const
+    { return jacobian_[domainScvIdx][rangeScvIdx]; }
+
+    /*!
+     * \brief Returns the local residual of a sub-control volume.
+     *
+     * \param dofIdx The local index of the sub control volume
+     */
+    const ScalarVectorBlock& residual(unsigned dofIdx) const
+    { return residual_[dofIdx]; }
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+
+    const Simulator& simulator_() const
+    { return *simulatorPtr_; }
+    const Problem& problem_() const
+    { return simulatorPtr_->problem(); }
+    const Model& model_() const
+    { return simulatorPtr_->model(); }
+
+    /*!
+     * \brief Resize all internal attributes to the size of the
+     *        element.
+     */
+    void resize_(const ElementContext& elemCtx)
+    {
+        size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+        size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+
+        residual_.resize(numDof);
+        if (jacobian_.N() != numDof || jacobian_.M() != numPrimaryDof)
+          jacobian_.setSize(numDof, numPrimaryDof);
+    }
+
+    /*!
+     * \brief Reset the all relevant internal attributes to 0
+     */
+    void reset_(const ElementContext&)
+    {
+        residual_ = 0.0;
+        jacobian_ = 0.0;
+    }
+
+    /*!
+     * \brief Updates the current local Jacobian matrix with the partial derivatives of
+     *        all equations for the degree of freedom associated with 'focusDofIdx'.
+     */
+    void updateLocalLinearization_(const ElementContext& elemCtx,
+                                   unsigned focusDofIdx)
+    {
+        const auto& resid = localResidual_.residual();
+
+        for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++)
+            residual_[focusDofIdx][eqIdx] = resid[focusDofIdx][eqIdx].value();
+
+        size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+        for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
+            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'
+                    jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = resid[dofIdx][eqIdx].derivative(pvIdx);
+                    Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
+                }
+            }
+        }
+    }
+
+    Simulator *simulatorPtr_;
+    Model *modelPtr_;
+
+    ElementContext *internalElemContext_;
+
+    LocalResidual localResidual_;
+
+    ScalarLocalBlockVector residual_;
+    ScalarLocalBlockMatrix jacobian_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaselocalresidualtpfa.hh

@@ -0,0 +1,648 @@
+// -*- 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::FvBaseLocalResidual
+ */
+#ifndef EWOMS_FV_BASE_LOCAL_TPFA_RESIDUAL_HH
+#define EWOMS_FV_BASE_LOCAL_TPFA_RESIDUAL_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/models/utils/parametersystem.hh>
+#include <opm/models/utils/alignedallocator.hh>
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+
+#include <dune/istl/bvector.hh>
+#include <dune/grid/common/geometry.hh>
+
+#include <dune/common/fvector.hh>
+
+#include <dune/common/classname.hh>
+
+#include <cmath>
+
+namespace Opm {
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Element-wise caculation of the residual matrix for models based on a finite
+ *        volume spatial discretization.
+ *
+ * \copydetails Doxygen::typeTagTParam
+ */
+template<class TypeTag>
+class FvBaseLocalResidualTPFA
+{
+private:
+    using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
+
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Element = typename GridView::template Codim<0>::Entity;
+
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using BoundaryContext = GetPropType<TypeTag, Properties::BoundaryContext>;
+
+    static constexpr bool useVolumetricResidual = getPropValue<TypeTag, Properties::UseVolumetricResidual>();
+
+    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
+    enum { extensiveStorageTerm = getPropValue<TypeTag, Properties::ExtensiveStorageTerm>() };
+
+    using Toolbox = MathToolbox<Evaluation>;
+    using EvalVector = Dune::FieldVector<Evaluation, numEq>;
+
+    // copying the local residual class is not a good idea
+    FvBaseLocalResidualTPFA(const FvBaseLocalResidualTPFA& )
+    {}
+
+public:
+    using LocalEvalBlockVector = Dune::BlockVector<EvalVector, aligned_allocator<EvalVector, alignof(EvalVector)> >;
+
+    FvBaseLocalResidualTPFA()
+    { }
+
+    ~FvBaseLocalResidualTPFA()
+    { }
+
+    /*!
+     * \brief Register all run-time parameters for the local residual.
+     */
+    static void registerParameters()
+    { }
+
+    /*!
+     * \brief Return the result of the eval() call using internal
+     *        storage.
+     */
+    const LocalEvalBlockVector& residual() const
+    { return internalResidual_; }
+
+    /*!
+     * \brief Return the result of the eval() call using internal
+     *        storage.
+     *
+     * \copydetails Doxygen::ecfvScvIdxParam
+     */
+    const EvalVector& residual(unsigned dofIdx) const
+    { return internalResidual_[dofIdx]; }
+
+    /*!
+     * \brief Compute the local residual, i.e. the deviation of the
+     *        conservation equations from zero and store the results
+     *        internally.
+     *
+     * The results can be requested afterwards using the residual() method.
+     *
+     * \copydetails Doxygen::problemParam
+     * \copydetails Doxygen::elementParam
+     */
+    void eval(const Problem& problem, const Element& element)
+    {
+        ElementContext elemCtx(problem);
+        elemCtx.updateAll(element);
+        eval(elemCtx);
+    }
+
+    /*!
+     * \brief Compute the local residual, i.e. the deviation of the
+     *        conservation equations from zero and store the results
+     *        internally.
+     *
+     * The results can be requested afterwards using the residual() method.
+     *
+     * \copydetails Doxygen::ecfvElemCtxParam
+     */
+    void eval(ElementContext& elemCtx)
+    {
+        size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+        internalResidual_.resize(numDof);
+        asImp_().eval(internalResidual_, elemCtx);
+    }
+
+    /*!
+     * \brief Compute the local residual, i.e. the deviation of the
+     *        conservation equations from zero.
+     *
+     * \copydetails Doxygen::residualParam
+     * \copydetails Doxygen::ecfvElemCtxParam
+     */
+    void eval(LocalEvalBlockVector& residual,
+              ElementContext& elemCtx) const
+    {
+        assert(residual.size() == elemCtx.numDof(/*timeIdx=*/0));
+
+        residual = 0.0;
+
+        // evaluate the flux terms
+        asImp_().evalFluxes(residual, elemCtx, /*timeIdx=*/0);
+
+        // evaluate the storage and the source terms
+        asImp_().evalVolumeTerms_(residual, elemCtx);
+
+        // evaluate the boundary conditions
+        //asImp_().evalBoundary_(residual, elemCtx, /*timeIdx=*/0);
+
+        if (useVolumetricResidual) {
+            // make the residual volume specific (i.e., make it incorrect mass per cubic
+            // meter instead of total mass)
+            size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+            for (unsigned dofIdx=0; dofIdx < numDof; ++dofIdx) {
+                if (elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0) > 0.0) {
+                    // interior DOF
+                    Scalar dofVolume = elemCtx.dofTotalVolume(dofIdx, /*timeIdx=*/0);
+
+                    assert(std::isfinite(dofVolume));
+                    Valgrind::CheckDefined(dofVolume);
+
+                    for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
+                        residual[dofIdx][eqIdx] /= dofVolume;
+                }
+            }
+        }
+    }
+
+    /*!
+     * \brief Calculate the amount of all conservation quantities stored in all element's
+     *        sub-control volumes for a given history index.
+     *
+     * This is used to figure out how much of each conservation quantity is inside the
+     * element.
+     *
+     * \copydetails Doxygen::storageParam
+     * \copydetails Doxygen::ecfvElemCtxParam
+     * \copydetails Doxygen::timeIdxParam
+     */
+    void evalStorage(LocalEvalBlockVector& storage,
+                     const ElementContext& elemCtx,
+                     unsigned timeIdx) const
+    {
+        // the derivative of the storage term depends on the current primary variables;
+        // for time indices != 0, the storage term is constant (because these solutions
+        // are not changed by the Newton method!)
+        if (timeIdx == 0) {
+            // calculate the amount of conservation each quantity inside
+            // all primary sub control volumes
+            size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+            for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
+                storage[dofIdx] = 0.0;
+
+                // the volume of the associated DOF
+                Scalar alpha =
+                    elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume();
+                    //* elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
+
+                // If the degree of freedom which we currently look at is the one at the
+                // center of attention, we need to consider the derivatives for the
+                // storage term, else the storage term is constant w.r.t. the primary
+                // variables of the focused DOF.
+                if (dofIdx == elemCtx.focusDofIndex()) {
+                    asImp_().computeStorage(storage[dofIdx],
+                                            elemCtx,
+                                            dofIdx,
+                                            timeIdx);
+
+                    for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
+                        storage[dofIdx][eqIdx] *= alpha;
+                }
+                else {
+                    Dune::FieldVector<Scalar, numEq> tmp;
+                    asImp_().computeStorage(tmp,
+                                            elemCtx,
+                                            dofIdx,
+                                            timeIdx);
+
+                    for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
+                        storage[dofIdx][eqIdx] = tmp[eqIdx]*alpha;
+                }
+            }
+        }
+        else {
+            // for all previous solutions, the storage term does _not_ depend on the
+            // current primary variables, so we use scalars to store it.
+            if (elemCtx.enableStorageCache()) {
+                size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
+                for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
+                    unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
+                    const auto& cachedStorage = elemCtx.model().cachedStorage(globalDofIdx, timeIdx);
+                    for (unsigned eqIdx=0; eqIdx < numEq; eqIdx++)
+                        storage[dofIdx][eqIdx] = cachedStorage[eqIdx];
+                }
+            }
+            else {
+                // calculate the amount of conservation each quantity inside
+                // all primary sub control volumes
+                Dune::FieldVector<Scalar, numEq> tmp;
+                size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+                for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
+                    tmp = 0.0;
+                    asImp_().computeStorage(tmp,
+                                            elemCtx,
+                                            dofIdx,
+                                            timeIdx);
+                    tmp *=
+                        elemCtx.stencil(timeIdx).subControlVolume(dofIdx).volume()
+                        * elemCtx.intensiveQuantities(dofIdx, timeIdx).extrusionFactor();
+
+                    for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                        storage[dofIdx][eqIdx] = tmp[eqIdx];
+                }
+            }
+        }
+
+#ifndef NDEBUG
+        size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+        for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                Valgrind::CheckDefined(storage[dofIdx][eqIdx]);
+                assert(isfinite(storage[dofIdx][eqIdx]));
+            }
+        }
+#endif
+    }
+
+    /*!
+     * \brief Add the flux term to a local residual.
+     *
+     * \copydetails Doxygen::residualParam
+     * \copydetails Doxygen::ecfvElemCtxParam
+     * \copydetails Doxygen::timeIdxParam
+     */
+    void evalFluxes(LocalEvalBlockVector& residual,
+                    const ElementContext& elemCtx,
+                    unsigned timeIdx) const
+    {
+        RateVector flux;
+
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        // calculate the mass flux over the sub-control volume faces
+        size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
+        for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
+            const auto& face = stencil.interiorFace(scvfIdx);
+            unsigned i = face.interiorIndex();
+            unsigned j = face.exteriorIndex();
+
+            Valgrind::SetUndefined(flux);
+            asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx);
+            Valgrind::CheckDefined(flux);
+#ifndef NDEBUG
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                assert(isfinite(flux[eqIdx]));
+#endif
+
+            // Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor();
+            Scalar alpha = face.area();
+            // Valgrind::CheckDefined(alpha);
+            // assert(alpha > 0.0);
+            // assert(isfinite(alpha));
+
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
+                 flux[eqIdx] *= alpha;
+
+            // The balance equation for a finite volume is given by
+            //
+            // dStorage/dt + Flux = Source
+            //
+            // where the 'Flux' and the 'Source' terms represent the
+            // mass per second which leaves the finite
+            // volume. Re-arranging this, we get
+            //
+            // dStorage/dt + Flux - Source = 0
+            //
+            // Since the mass flux as calculated by computeFlux() goes out of sub-control
+            // volume i and into sub-control volume j, we need to add the flux to finite
+            // volume i and subtract it from finite volume j
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                assert(isfinite(flux[eqIdx]));
+                residual[i][eqIdx] += flux[eqIdx];
+                residual[j][eqIdx] -= flux[eqIdx];
+            }
+        }
+
+#if !defined NDEBUG
+        // in debug mode, ensure that the residual is well-defined
+        size_t numDof = elemCtx.numDof(timeIdx);
+        for (unsigned i=0; i < numDof; i++) {
+            for (unsigned j = 0; j < numEq; ++ j) {
+                assert(isfinite(residual[i][j]));
+                Valgrind::CheckDefined(residual[i][j]);
+            }
+        }
+#endif
+
+    }
+
+    /////////////////////////////
+    // The following methods _must_ be overloaded by the actual flow
+    // models!
+    /////////////////////////////
+
+    /*!
+     * \brief Evaluate the amount all conservation quantities
+     *        (e.g. phase mass) within a finite sub-control volume.
+     *
+     * \copydetails Doxygen::storageParam
+     * \copydetails Doxygen::ecfvScvCtxParams
+     */
+    void computeStorage(EqVector&,
+                        const ElementContext&,
+                        unsigned,
+                        unsigned) const
+    {
+        throw std::logic_error("Not implemented: The local residual "+Dune::className<Implementation>()
+                               +" does not implement the required method 'computeStorage()'");
+    }
+
+    /*!
+     * \brief Evaluates the total mass flux of all conservation
+     *        quantities over a face of a sub-control volume.
+     *
+     * \copydetails Doxygen::areaFluxParam
+     * \copydetails Doxygen::ecfvScvfCtxParams
+     */
+    void computeFlux(RateVector&,
+                     const ElementContext&,
+                     unsigned,
+                     unsigned) const
+    {
+        throw std::logic_error("Not implemented: The local residual "+Dune::className<Implementation>()
+                               +" does not implement the required method 'computeFlux()'");
+    }
+
+    /*!
+     * \brief Calculate the source term of the equation
+     *
+     * \copydoc Doxygen::sourceParam
+     * \copydoc Doxygen::ecfvScvCtxParams
+     */
+    void computeSource(RateVector&,
+                       const ElementContext&,
+                       unsigned,
+                       unsigned) const
+    {
+        throw std::logic_error("Not implemented: The local residual "+Dune::className<Implementation>()
+                               +" does not implement the required method 'computeSource()'");
+    }
+
+protected:
+    /*!
+     * \brief Evaluate the boundary conditions of an element.
+     */
+    void evalBoundary_(LocalEvalBlockVector& residual,
+                       const ElementContext& elemCtx,
+                       unsigned timeIdx) const
+    {
+        
+        if (!elemCtx.onBoundary())
+            return;
+        throw std::logic_error("Not implemented: Boundary??? "+Dune::className<Implementation>()
+                               +" does not implement the required method 'computeSource()'");
+
+        BoundaryContext boundaryCtx(elemCtx);
+        // move the iterator to the first boundary
+        if(boundaryCtx.intersection(0).neighbor())
+            boundaryCtx.increment();
+
+        // evaluate the boundary for all boundary faces of the current context
+        size_t numBoundaryFaces = boundaryCtx.numBoundaryFaces(/*timeIdx=*/0);
+        for (unsigned faceIdx = 0; faceIdx < numBoundaryFaces; ++faceIdx, boundaryCtx.increment()) {
+            // add the residual of all vertices of the boundary
+            // segment
+            evalBoundarySegment_(residual,
+                                 boundaryCtx,
+                                 faceIdx,
+                                 timeIdx);
+        }
+
+#if !defined NDEBUG
+        // in debug mode, ensure that the residual and the storage terms are well-defined
+        size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+        for (unsigned i=0; i < numDof; i++) {
+            for (unsigned j = 0; j < numEq; ++ j) {
+                assert(isfinite(residual[i][j]));
+                Valgrind::CheckDefined(residual[i][j]);
+            }
+        }
+#endif
+
+    }
+
+    /*!
+     * \brief Evaluate all boundary conditions for a single
+     *        sub-control volume face to the local residual.
+     */
+    void evalBoundarySegment_(LocalEvalBlockVector& residual,
+                              const BoundaryContext& boundaryCtx,
+                              unsigned boundaryFaceIdx,
+                              unsigned timeIdx) const
+    {
+        throw std::logic_error("Not implemented: Boundary??? "+Dune::className<Implementation>()
+                               +" does not implement the required method 'computeSource()'");
+
+        BoundaryRateVector values;
+
+        Valgrind::SetUndefined(values);
+        boundaryCtx.problem().boundary(values, boundaryCtx, boundaryFaceIdx, timeIdx);
+        Valgrind::CheckDefined(values);
+
+        const auto& stencil = boundaryCtx.stencil(timeIdx);
+        unsigned dofIdx = stencil.boundaryFace(boundaryFaceIdx).interiorIndex();
+        const auto& insideIntQuants = boundaryCtx.elementContext().intensiveQuantities(dofIdx, timeIdx);
+        for (unsigned eqIdx = 0; eqIdx < values.size(); ++eqIdx)  {
+            values[eqIdx] *=
+                stencil.boundaryFace(boundaryFaceIdx).area()
+                * insideIntQuants.extrusionFactor();
+
+            Valgrind::CheckDefined(values[eqIdx]);
+            assert(isfinite(values[eqIdx]));
+        }
+
+        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+            residual[dofIdx][eqIdx] += values[eqIdx];
+    }
+
+    /*!
+     * \brief Add the change in the storage terms and the source term
+     *        to the local residual of all sub-control volumes of the
+     *        current element.
+     */
+    void evalVolumeTerms_(LocalEvalBlockVector& residual,
+                          ElementContext& elemCtx) const
+    {
+        EvalVector tmp;
+        EqVector tmp2;
+        RateVector sourceRate;
+
+        tmp = 0.0;
+        tmp2 = 0.0;
+
+        // evaluate the volumetric terms (storage + source terms)
+        size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+        for (unsigned dofIdx=0; dofIdx < numPrimaryDof; dofIdx++) {
+            // Scalar extrusionFactor =
+            //     elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).extrusionFactor();
+            // Valgrind::CheckDefined(extrusionFactor);
+            // assert(isfinite(extrusionFactor));
+            // assert(extrusionFactor > 0.0);
+            Scalar scvVolume =
+                elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume();// * extrusionFactor;
+            Valgrind::CheckDefined(scvVolume);
+            assert(isfinite(scvVolume));
+            assert(scvVolume > 0.0);
+
+            // if the model uses extensive quantities in its storage term, and we use
+            // automatic differention and current DOF is also not the one we currently
+            // focus on, the storage term does not need any derivatives!
+            if (!extensiveStorageTerm &&
+                !std::is_same<Scalar, Evaluation>::value &&
+                dofIdx != elemCtx.focusDofIndex())
+            {
+                asImp_().computeStorage(tmp2, elemCtx, dofIdx, /*timeIdx=*/0);
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                    tmp[eqIdx] = tmp2[eqIdx];
+            }
+            else
+                asImp_().computeStorage(tmp, elemCtx, dofIdx, /*timeIdx=*/0);
+
+#ifndef NDEBUG
+            Valgrind::CheckDefined(tmp);
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                assert(isfinite(tmp[eqIdx]));
+#endif
+
+            if (elemCtx.enableStorageCache()) {
+                const auto& model = elemCtx.model();
+                unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
+                if (model.newtonMethod().numIterations() == 0 &&
+                    !elemCtx.haveStashedIntensiveQuantities())
+                {
+                    if (!elemCtx.problem().recycleFirstIterationStorage()) {
+                        // we re-calculate the storage term for the solution of the
+                        // previous time step from scratch instead of using the one of
+                        // the first iteration of the current time step.
+                        tmp2 = 0.0;
+                        elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/1);
+                        asImp_().computeStorage(tmp2, elemCtx,  dofIdx, /*timeIdx=*/1);
+                    }
+                    else {
+                        // if the storage term is cached and we're in the first iteration
+                        // of the time step, use the storage term of the first iteration
+                        // as the one as the solution of the last time step (this assumes
+                        // that the initial guess for the solution at the end of the time
+                        // step is the same as the solution at the beginning of the time
+                        // step. This is usually true, but some fancy preprocessing
+                        // scheme might invalidate that assumption.)
+                        for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
+                            tmp2[eqIdx] = Toolbox::value(tmp[eqIdx]);
+                    }
+
+                    Valgrind::CheckDefined(tmp2);
+
+                    model.updateCachedStorage(globalDofIdx, /*timeIdx=*/1, tmp2);
+                }
+                else {
+                    // if the mass storage at the beginning of the time step is not cached,
+                    // if the storage term is cached and we're not looking at the first
+                    // iteration of the time step, we take the cached data.
+                    tmp2 = model.cachedStorage(globalDofIdx, /*timeIdx=*/1);
+                    Valgrind::CheckDefined(tmp2);
+                }
+            }
+            else {
+                // if the mass storage at the beginning of the time step is not cached,
+                // we re-calculate it from scratch.
+                tmp2 = 0.0;
+                asImp_().computeStorage(tmp2, elemCtx,  dofIdx, /*timeIdx=*/1);
+                Valgrind::CheckDefined(tmp2);
+            }
+
+            // Use the implicit Euler time discretization
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                double dt = elemCtx.simulator().timeStepSize();
+                assert(dt > 0);
+                tmp[eqIdx] -= tmp2[eqIdx];
+                tmp[eqIdx] *= scvVolume / dt;
+
+                residual[dofIdx][eqIdx] += tmp[eqIdx];
+            }
+
+            Valgrind::CheckDefined(residual[dofIdx]);
+
+            // deal with the source term
+            asImp_().computeSource(sourceRate, elemCtx, dofIdx, /*timeIdx=*/0);
+
+            // if the model uses extensive quantities in its storage term, and we use
+            // automatic differention and current DOF is also not the one we currently
+            // focus on, the storage term does not need any derivatives!
+            if (!extensiveStorageTerm &&
+                !std::is_same<Scalar, Evaluation>::value &&
+                dofIdx != elemCtx.focusDofIndex())
+            {
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                    residual[dofIdx][eqIdx] -= scalarValue(sourceRate[eqIdx])*scvVolume;
+            }
+            else {
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                    sourceRate[eqIdx] *= scvVolume;
+                    residual[dofIdx][eqIdx] -= sourceRate[eqIdx];
+                }
+            }
+
+            Valgrind::CheckDefined(residual[dofIdx]);
+        }
+
+#if !defined NDEBUG
+        // in debug mode, ensure that the residual is well-defined
+        size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+        for (unsigned i=0; i < numDof; i++) {
+            for (unsigned j = 0; j < numEq; ++ j) {
+                assert(isfinite(residual[i][j]));
+                Valgrind::CheckDefined(residual[i][j]);
+            }
+        }
+#endif
+    }
+
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+
+    LocalEvalBlockVector internalResidual_;
+};
+
+} // namespace Opm
+
+#endif