changed: ewoms/disc -> opm/models/discretization

examples/co2injection_flash_ecfv.cpp

@@ -34,7 +34,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/flash/flashmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionflash.hh"
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_flash_ni_ecfv.cpp

@@ -31,7 +31,7 @@
 #include <opm/material/common/quad.hpp>
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/flash/flashmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionflash.hh"
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_flash_ni_vcfv.cpp

@@ -31,7 +31,7 @@
 #include <opm/material/common/quad.hpp>
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/flash/flashmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/co2injectionflash.hh"
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_flash_vcfv.cpp

@@ -34,7 +34,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/flash/flashmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/co2injectionflash.hh"
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_immiscible_ecfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/immiscible/immisciblemodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_immiscible_ni_ecfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/immiscible/immisciblemodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_immiscible_ni_vcfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/immiscible/immisciblemodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_immiscible_vcfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/immiscible/immisciblemodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_ncp_ecfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/ncp/ncpmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_ncp_ni_ecfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/ncp/ncpmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_ncp_ni_vcfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/ncp/ncpmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_ncp_vcfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/ncp/ncpmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 
 #include "problems/co2injectionproblem.hh"
 

examples/co2injection_pvs_ecfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/pvs/pvsmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_pvs_ni_ecfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/pvs/pvsmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_pvs_ni_vcfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/pvs/pvsmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/co2injectionproblem.hh"
 
 BEGIN_PROPERTIES

examples/co2injection_pvs_vcfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/pvs/pvsmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 
 #include "problems/co2injectionproblem.hh"
 

examples/finger_immiscible_ecfv.cpp

@@ -29,7 +29,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/immiscible/immisciblemodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/fingerproblem.hh"
 
 BEGIN_PROPERTIES

examples/finger_immiscible_vcfv.cpp

@@ -29,7 +29,7 @@
 
 #include <ewoms/common/start.hh>
 #include <ewoms/models/immiscible/immisciblemodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/fingerproblem.hh"
 
 BEGIN_PROPERTIES

examples/lens_immiscible_ecfv_ad.hh

@@ -30,7 +30,7 @@
 #define EWOMS_LENS_IMMISCIBLE_ECFV_AD_HH
 
 #include <ewoms/models/immiscible/immisciblemodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/lensproblem.hh"
 
 BEGIN_PROPERTIES

examples/lens_richards_ecfv.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 
 #include "problems/richardslensproblem.hh"
 

examples/lens_richards_vcfv.cpp

@@ -28,7 +28,7 @@
 #include "config.h"
 
 #include <opm/models/utils/start.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 
 #include "problems/richardslensproblem.hh"
 

examples/problems/fingerproblem.hh

@@ -42,7 +42,7 @@
 #include <opm/material/components/Air.hpp>
 
 #include <ewoms/models/immiscible/immiscibleproperties.hh>
-#include <ewoms/disc/common/restrictprolong.hh>
+#include <opm/models/discretization/common/restrictprolong.hh>
 
 #if HAVE_DUNE_ALUGRID
 #include <dune/alugrid/grid.hh>

examples/problems/lensproblem.hh

@@ -30,8 +30,8 @@
 
 #include <ewoms/io/structuredgridvanguard.hh>
 #include <ewoms/models/immiscible/immiscibleproperties.hh>
-#include <ewoms/disc/common/fvbaseadlocallinearizer.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/common/fvbaseadlocallinearizer.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 
 #include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>

examples/reservoir_blackoil_ecfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <opm/models/blackoil/blackoilmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/reservoirproblem.hh"
 
 BEGIN_PROPERTIES

examples/reservoir_blackoil_vcfv.cpp

@@ -29,7 +29,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <opm/models/blackoil/blackoilmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/reservoirproblem.hh"
 
 BEGIN_PROPERTIES

examples/reservoir_ncp_ecfv.cpp

@@ -29,7 +29,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/ncp/ncpmodel.hh>
-#include <ewoms/disc/ecfv/ecfvdiscretization.hh>
+#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
 #include "problems/reservoirproblem.hh"
 
 BEGIN_PROPERTIES

examples/reservoir_ncp_vcfv.cpp

@@ -30,7 +30,7 @@
 
 #include <opm/models/utils/start.hh>
 #include <ewoms/models/ncp/ncpmodel.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 #include "problems/reservoirproblem.hh"
 
 BEGIN_PROPERTIES

examples/tutorial1problem.hh

@@ -32,7 +32,7 @@
 #include <ewoms/models/immiscible/immisciblemodel.hh>
 
 // The spatial discretization (VCFV == Vertex-Centered Finite Volumes)
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>  /*@\label{tutorial1:include-discretization}@*/
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>  /*@\label{tutorial1:include-discretization}@*/
 
 // The chemical species that are used
 #include <opm/material/components/SimpleH2O.hpp>

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -34,7 +34,7 @@
 #include "blackoilenergymodules.hh"
 #include "blackoilfoammodules.hh"
 
-#include <ewoms/disc/common/fvbaseprimaryvariables.hh>
+#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
 
 #include <dune/common/fvector.hh>
 

opm/models/common/diffusionmodule.hh

@@ -28,7 +28,7 @@
 #ifndef EWOMS_DIFFUSION_MODULE_HH
 #define EWOMS_DIFFUSION_MODULE_HH
 
-#include <ewoms/disc/common/fvbaseproperties.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
 #include <opm/models/common/quantitycallbacks.hh>
 
 #include <opm/material/common/Valgrind.hpp>

opm/models/common/energymodule.hh

@@ -29,7 +29,7 @@
 #ifndef EWOMS_ENERGY_MODULE_HH
 #define EWOMS_ENERGY_MODULE_HH
 
-#include <ewoms/disc/common/fvbaseproperties.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
 #include <opm/models/common/quantitycallbacks.hh>
 
 #include <opm/material/common/Valgrind.hpp>

opm/models/common/forchheimerfluxmodule.hh

@@ -32,7 +32,7 @@
 
 #include "darcyfluxmodule.hh"
 
-#include <ewoms/disc/common/fvbaseproperties.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
 
 #include <opm/material/common/Valgrind.hpp>
 #include <opm/material/common/Unused.hpp>

opm/models/common/multiphasebaseextensivequantities.hh

@@ -31,7 +31,7 @@
 #include "multiphasebaseproperties.hh"
 
 #include <opm/models/common/quantitycallbacks.hh>
-#include <ewoms/disc/common/fvbaseextensivequantities.hh>
+#include <opm/models/discretization/common/fvbaseextensivequantities.hh>
 #include <opm/models/utils/parametersystem.hh>
 
 #include <opm/material/common/Valgrind.hpp>

opm/models/common/multiphasebasemodel.hh

@@ -35,7 +35,7 @@
 #include "multiphasebaseextensivequantities.hh"
 
 #include <opm/models/common/flux.hh>
-#include <ewoms/disc/vcfv/vcfvdiscretization.hh>
+#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
 
 #include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>

opm/models/common/multiphasebaseproblem.hh

@@ -30,8 +30,8 @@
 
 #include "multiphasebaseproperties.hh"
 
-#include <ewoms/disc/common/fvbaseproblem.hh>
-#include <ewoms/disc/common/fvbaseproperties.hh>
+#include <opm/models/discretization/common/fvbaseproblem.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
 
 #include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
 #include <opm/material/common/Means.hpp>

opm/models/common/multiphasebaseproperties.hh

@@ -30,7 +30,7 @@
 #ifndef EWOMS_MULTI_PHASE_BASE_PROPERTIES_HH
 #define EWOMS_MULTI_PHASE_BASE_PROPERTIES_HH
 
-#include <ewoms/disc/common/fvbaseproperties.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
 #include <ewoms/io/vtkmultiphasemodule.hh>
 #include <ewoms/io/vtktemperaturemodule.hh>
 

opm/models/common/quantitycallbacks.hh

@@ -29,7 +29,7 @@
 #ifndef EWOMS_QUANTITY_CALLBACKS_HH
 #define EWOMS_QUANTITY_CALLBACKS_HH
 
-#include <ewoms/disc/common/fvbaseproperties.hh>
+#include <opm/models/discretization/common/fvbaseproperties.hh>
 
 #include <opm/material/common/MathToolbox.hpp>
 #include <opm/material/common/Valgrind.hpp>

opm/models/discretization/common/baseauxiliarymodule.hh

@@ -0,0 +1,138 @@
+// -*- 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::BaseAuxiliaryModule
+ */
+#ifndef EWOMS_BASE_AUXILIARY_MODULE_HH
+#define EWOMS_BASE_AUXILIARY_MODULE_HH
+
+#include <opm/models/utils/propertysystem.hh>
+
+#include <opm/models/discretization/common/fvbaseproperties.hh>
+
+#include <set>
+#include <vector>
+
+BEGIN_PROPERTIES
+
+NEW_TYPE_TAG(AuxModule);
+
+// declare the properties required by the for the ecl grid manager
+NEW_PROP_TAG(Grid);
+NEW_PROP_TAG(GridView);
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(DofMapper);
+NEW_PROP_TAG(GlobalEqVector);
+NEW_PROP_TAG(SparseMatrixAdapter);
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup ModelModules
+ *
+ * \brief Base class for specifying auxiliary equations.
+ *
+ * For example, these equations can be wells, non-neighboring connections, interfaces
+ * between model domains, etc.
+ */
+template <class TypeTag>
+class BaseAuxiliaryModule
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
+
+protected:
+    typedef std::set<unsigned> NeighborSet;
+
+public:
+    virtual ~BaseAuxiliaryModule()
+    {}
+
+    /*!
+     * \brief Returns the number of additional degrees of freedom required for the
+     *        auxiliary module.
+     */
+    virtual unsigned numDofs() const = 0;
+
+    /*!
+     * \brief Set the offset in the global system of equations for the first degree of
+     *        freedom of this auxiliary module.
+     */
+    void setDofOffset(int value)
+    { dofOffset_ = value; }
+
+    /*!
+     * \brief Return the offset in the global system of equations for the first degree of
+     *        freedom of this auxiliary module.
+     */
+    int dofOffset()
+    { return dofOffset_; }
+
+    /*!
+     * \brief Given a degree of freedom relative to the current auxiliary equation,
+     *        return the corresponding index in the global system of equations.
+     */
+    int localToGlobalDof(unsigned localDofIdx) const
+    {
+        assert(0 <= localDofIdx);
+        assert(localDofIdx < numDofs());
+        return dofOffset_ + localDofIdx;
+    }
+
+    /*!
+     * \brief Specify the additional neighboring correlations caused by the auxiliary
+     *        module.
+     */
+    virtual void addNeighbors(std::vector<NeighborSet>& neighbors) const = 0;
+
+    /*!
+     * \brief Set the initial condition of the auxiliary module in the solution vector.
+     */
+    virtual void applyInitial() = 0;
+
+    /*!
+     * \brief Linearize the auxiliary equation.
+     */
+    virtual void linearize(SparseMatrixAdapter& matrix, GlobalEqVector& residual) = 0;
+
+    /*!
+     * \brief This method is called after the linear solver has been called but before
+     *        the solution is updated for the next iteration.
+     *
+     * It is intended to implement stuff like Schur complements.
+     */
+    virtual void postSolve(GlobalEqVector& residual OPM_UNUSED)
+    {};
+
+private:
+    int dofOffset_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseadlocallinearizer.hh

@@ -0,0 +1,316 @@
+// -*- 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_LINEARIZER_HH
+#define EWOMS_FV_BASE_AD_LOCAL_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 FvBaseAdLocalLinearizer;
+}
+
+BEGIN_PROPERTIES
+
+// declare the property tags required for the finite differences local linearizer
+NEW_TYPE_TAG(AutoDiffLocalLinearizer);
+
+NEW_PROP_TAG(LocalLinearizer);
+NEW_PROP_TAG(Evaluation);
+
+NEW_PROP_TAG(LocalResidual);
+NEW_PROP_TAG(Simulator);
+NEW_PROP_TAG(Problem);
+NEW_PROP_TAG(Model);
+NEW_PROP_TAG(PrimaryVariables);
+NEW_PROP_TAG(ElementContext);
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(Evaluation);
+NEW_PROP_TAG(GridView);
+
+// set the properties to be spliced in
+SET_TYPE_PROP(AutoDiffLocalLinearizer, LocalLinearizer,
+              Opm::FvBaseAdLocalLinearizer<TypeTag>);
+
+//! Set the function evaluation w.r.t. the primary variables
+SET_PROP(AutoDiffLocalLinearizer, Evaluation)
+{
+private:
+    static const unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+public:
+    typedef Opm::DenseAd::Evaluation<Scalar, numEq> type;
+};
+
+END_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 FvBaseAdLocalLinearizer
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizer) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+    typedef Dune::FieldVector<Scalar, numEq> ScalarVectorBlock;
+    // extract local matrices from jacobian matrix for consistency
+    typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter)::MatrixBlock ScalarMatrixBlock;
+
+    typedef Dune::BlockVector<ScalarVectorBlock> ScalarLocalBlockVector;
+    typedef Dune::Matrix<ScalarMatrixBlock> ScalarLocalBlockMatrix;
+
+public:
+    FvBaseAdLocalLinearizer()
+        : internalElemContext_(0)
+    { }
+
+    // copying local linearizer objects around is a very bad idea, so we explicitly
+    // prevent it...
+    FvBaseAdLocalLinearizer(const FvBaseAdLocalLinearizer&) = delete;
+
+    ~FvBaseAdLocalLinearizer()
+    { 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);
+        reset_(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();
+
+            // 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);
+        jacobian_.setSize(numDof, numPrimaryDof);
+    }
+
+    /*!
+     * \brief Reset the all relevant internal attributes to 0
+     */
+    void reset_(const ElementContext& elemCtx OPM_UNUSED)
+    {
+        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);
+                    Opm::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/fvbaseboundarycontext.hh

@@ -0,0 +1,269 @@
+// -*- 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::FvBaseBoundaryContext
+ */
+#ifndef EWOMS_FV_BASE_BOUNDARY_CONTEXT_HH
+#define EWOMS_FV_BASE_BOUNDARY_CONTEXT_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/material/common/Unused.hpp>
+
+#include <dune/common/fvector.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Represents all quantities which available on boundary segments
+ */
+template<class TypeTag>
+class FvBaseBoundaryContext
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, GradientCalculator) GradientCalculator;
+
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::IntersectionIterator IntersectionIterator;
+    typedef typename GridView::Intersection Intersection;
+
+    enum { dim = GridView::dimension };
+    enum { dimWorld = GridView::dimensionworld };
+
+    typedef typename GridView::ctype CoordScalar;
+    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+    typedef Dune::FieldVector<Scalar, dimWorld> Vector;
+
+public:
+    /*!
+     * \brief The constructor.
+     */
+    explicit FvBaseBoundaryContext(const ElementContext& elemCtx)
+        : elemCtx_(elemCtx)
+        , intersectionIt_(gridView().ibegin(element()))
+    { }
+
+    void increment()
+    {
+        const auto& iend = gridView().iend(element());
+        while (intersectionIt_ != iend && !intersectionIt_->boundary())
+            ++ intersectionIt_;
+    }
+
+    /*!
+     * \copydoc Opm::ElementContext::problem()
+     */
+    const Problem& problem() const
+    { return elemCtx_.problem(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::model()
+     */
+    const Model& model() const
+    { return elemCtx_.model(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::gridView()
+     */
+    const GridView& gridView() const
+    { return elemCtx_.gridView(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::element()
+     */
+    const Element& element() const
+    { return elemCtx_.element(); }
+
+    /*!
+     * \brief Returns a reference to the element context object.
+     */
+    const ElementContext& elementContext() const
+    { return elemCtx_; }
+
+    /*!
+     * \brief Returns a reference to the current gradient calculator.
+     */
+    const GradientCalculator& gradientCalculator() const
+    { return elemCtx_.gradientCalculator(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::numDof()
+     */
+    size_t numDof(unsigned timeIdx) const
+    { return elemCtx_.numDof(timeIdx); }
+
+    /*!
+     * \copydoc Opm::ElementContext::numPrimaryDof()
+     */
+    size_t numPrimaryDof(unsigned timeIdx) const
+    { return elemCtx_.numPrimaryDof(timeIdx); }
+
+    /*!
+     * \copydoc Opm::ElementContext::numInteriorFaces()
+     */
+    size_t numInteriorFaces(unsigned timeIdx) const
+    { return elemCtx_.numInteriorFaces(timeIdx); }
+
+    /*!
+     * \brief Return the number of boundary segments of the current element
+     */
+    size_t numBoundaryFaces(unsigned timeIdx) const
+    { return elemCtx_.stencil(timeIdx).numBoundaryFaces(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::stencil()
+     */
+    const Stencil& stencil(unsigned timeIdx) const
+    { return elemCtx_.stencil(timeIdx); }
+
+    /*!
+     * \brief Returns the outer unit normal of the boundary segment
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     * \param timeIdx The index of the solution used by the time discretization
+     */
+    Vector normal(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    {
+        auto tmp = stencil(timeIdx).boundaryFace[boundaryFaceIdx].normal;
+        tmp /= tmp.two_norm();
+        return tmp;
+    }
+
+    /*!
+     * \brief Returns the area [m^2] of a given boudary segment.
+     */
+    Scalar boundarySegmentArea(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    { return elemCtx_.stencil(timeIdx).boundaryFace(boundaryFaceIdx).area(); }
+
+    /*!
+     * \brief Return the position of a local entity in global coordinates.
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     * \param timeIdx The index of the solution used by the time discretization
+     */
+    const GlobalPosition& pos(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    { return stencil(timeIdx).boundaryFace(boundaryFaceIdx).integrationPos(); }
+
+    /*!
+     * \brief Return the position of a control volume's center in global coordinates.
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     * \param timeIdx The index of the solution used by the time discretization
+     */
+    const GlobalPosition& cvCenter(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    {
+        unsigned scvIdx = stencil(timeIdx).boundaryFace(boundaryFaceIdx).interiorIndex();
+        return stencil(timeIdx).subControlVolume(scvIdx).globalPos();
+    }
+
+    /*!
+     * \brief Return the local sub-control volume index upon which the linearization is
+     *        currently focused.
+     */
+    unsigned focusDofIndex() const
+    { return elemCtx_.focusDofIndex(); }
+
+    /*!
+     * \brief Return the local sub-control volume index of the
+     *        interior of a boundary segment
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     * \param timeIdx The index of the solution used by the time discretization
+     */
+    unsigned interiorScvIndex(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    { return stencil(timeIdx).boundaryFace(boundaryFaceIdx).interiorIndex(); }
+
+    /*!
+     * \brief Return the global space index of the sub-control volume
+     *        at the interior of a boundary segment
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     * \param timeIdx The index of the solution used by the time discretization
+     */
+    unsigned globalSpaceIndex(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    { return elemCtx_.globalSpaceIndex(interiorScvIndex(boundaryFaceIdx, timeIdx), timeIdx); }
+
+    /*!
+     * \brief Return the intensive quantities for the finite volume in the
+     *        interiour of a boundary segment
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     * \param timeIdx The index of the solution used by the time discretization
+     */
+    const IntensiveQuantities& intensiveQuantities(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    {
+        unsigned interiorScvIdx = this->interiorScvIndex(boundaryFaceIdx, timeIdx);
+        return elemCtx_.intensiveQuantities(interiorScvIdx, timeIdx);
+    }
+
+    /*!
+     * \brief Return the extensive quantities for a given boundary face.
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     * \param timeIdx The index of the solution used by the time discretization
+     */
+    const ExtensiveQuantities& extensiveQuantities(unsigned boundaryFaceIdx, unsigned timeIdx) const
+    { return elemCtx_.boundaryExtensiveQuantities(boundaryFaceIdx, timeIdx); }
+
+    /*!
+     * \brief Return the intersection for the neumann segment
+     *
+     * TODO/HACK: The intersection should take a local index as an
+     * argument. since that's not supported efficiently by the DUNE
+     * grid interface, we just ignore the index argument here!
+     *
+     * \param boundaryFaceIdx The local index of the boundary segment
+     */
+    const Intersection& intersection(unsigned boundaryFaceIdx OPM_UNUSED) const
+    { return *intersectionIt_; }
+
+    /*!
+     * \brief Return the intersection for the neumann segment
+     *
+     * TODO/HACK: the intersection iterator can basically be
+     * considered as an index which is manipulated externally, but
+     * context classes should not store any indices. it is done this
+     * way for performance reasons
+     */
+    IntersectionIterator& intersectionIt()
+    { return intersectionIt_; }
+
+protected:
+    const ElementContext& elemCtx_;
+    IntersectionIterator intersectionIt_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseconstraints.hh

@@ -0,0 +1,86 @@
+// -*- 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::FvBaseConstraints
+ */
+#ifndef EWOMS_FV_BASE_CONSTRAINTS_HH
+#define EWOMS_FV_BASE_CONSTRAINTS_HH
+
+#include <opm/models/utils/propertysystem.hh>
+#include <opm/material/common/Valgrind.hpp>
+
+BEGIN_PROPERTIES
+
+NEW_PROP_TAG(PrimaryVariables);
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Class to specify constraints for a finite volume spatial discretization.
+ *
+ * Note that eWoms does not support "partial" constraints. (I.e., either all equations
+ * must be constraint or none.)
+ */
+template <class TypeTag>
+class FvBaseConstraints : public GET_PROP_TYPE(TypeTag, PrimaryVariables)
+{
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) ParentType;
+
+public:
+    FvBaseConstraints()
+    { setActive(false); }
+
+    FvBaseConstraints(const FvBaseConstraints&) = default;
+
+//! \cond SKIP
+    /*!
+     * \brief Use the default assignment operator
+     */
+    FvBaseConstraints &operator=(const FvBaseConstraints&) = default;
+//! \endcond
+
+    /*!
+     * \brief Returns true if the constraints are enforced.
+     */
+    bool isActive() const
+    { return isActive_; }
+
+    /*!
+     * \brief Specify whether the constraints should be enforced or not
+     */
+    void setActive(bool yesno)
+    { isActive_ = yesno; }
+
+private:
+    bool isActive_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseconstraintscontext.hh

@@ -0,0 +1,118 @@
+// -*- 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::FvBaseConstraintsContext
+ */
+#ifndef EWOMS_FV_BASE_CONSTRAINTS_CONTEXT_HH
+#define EWOMS_FV_BASE_CONSTRAINTS_CONTEXT_HH
+
+#include "fvbaseproperties.hh"
+
+#include <dune/common/fvector.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Represents all quantities which available for calculating constraints
+ */
+template<class TypeTag>
+class FvBaseConstraintsContext
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    enum { dimWorld = GridView::dimensionworld };
+
+    typedef typename GridView::ctype CoordScalar;
+    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+
+public:
+    /*!
+     * \brief The constructor.
+     */
+    explicit FvBaseConstraintsContext(const ElementContext& elemCtx)
+        : elemCtx_(elemCtx)
+    { }
+
+    /*!
+     * \copydoc Opm::ElementContext::problem()
+     */
+    const Problem& problem() const
+    { return elemCtx_.problem(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::model()
+     */
+    const Model& model() const
+    { return elemCtx_.model(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::gridView()
+     */
+    const GridView& gridView() const
+    { return elemCtx_.gridView(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::element()
+     */
+    const Element& element() const
+    { return elemCtx_.element(); }
+
+    /*!
+     * \copydoc Opm::ElementContext::numDof()
+     */
+    int numDof(int timeIdx) const
+    { return elemCtx_.numDof(timeIdx); }
+
+    /*!
+     * \copydoc Opm::ElementContext::numInteriorFaces()
+     */
+    int numInteriorFaces(int timeIdx) const
+    { return elemCtx_.numInteriorFaces(timeIdx); }
+
+    /*!
+     * \copydoc Opm::ElementContext::globalSpaceIndex
+     */
+    int globalSpaceIndex(int dofIdx, int timeIdx) const
+    { return elemCtx_.globalSpaceIndex(dofIdx, timeIdx); }
+
+    /*!
+     * \copydoc Opm::ElementContext::pos
+     */
+    GlobalPosition pos(int dofIdx, int timeIdx) const
+    { return elemCtx_.pos(dofIdx, timeIdx); }
+
+protected:
+    const ElementContext& elemCtx_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseelementcontext.hh

@@ -0,0 +1,607 @@
+// -*- 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::FvBaseElementContext
+ */
+#ifndef EWOMS_FV_BASE_ELEMENT_CONTEXT_HH
+#define EWOMS_FV_BASE_ELEMENT_CONTEXT_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/models/utils/alignedallocator.hh>
+
+#include <opm/material/common/Unused.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <vector>
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief This class stores an array of IntensiveQuantities objects, one
+ *        intensive quantities object for each of the element's vertices
+ */
+template<class TypeTag>
+class FvBaseElementContext
+{
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
+    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
+
+    // the history size of the time discretization in number of steps
+    enum { timeDiscHistorySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) };
+
+    struct DofStore_ {
+        IntensiveQuantities intensiveQuantities[timeDiscHistorySize];
+        PrimaryVariables priVars[timeDiscHistorySize];
+        const IntensiveQuantities *thermodynamicHint[timeDiscHistorySize];
+    };
+    typedef std::vector<DofStore_> DofVarsVector;
+    typedef std::vector<ExtensiveQuantities> ExtensiveQuantitiesVector;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
+    typedef typename GET_PROP_TYPE(TypeTag, GradientCalculator) GradientCalculator;
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    static const unsigned dimWorld = GridView::dimensionworld;
+    static const unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
+
+    typedef typename GridView::ctype CoordScalar;
+    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+
+    // we don't allow copies of element contexts!
+    FvBaseElementContext(const FvBaseElementContext& ) = delete;
+
+public:
+    /*!
+     * \brief The constructor.
+     */
+    explicit FvBaseElementContext(const Simulator& simulator)
+        : gridView_(simulator.gridView())
+        , stencil_(gridView_, simulator.model().dofMapper() )
+    {
+        // remember the simulator object
+        simulatorPtr_ = &simulator;
+        enableStorageCache_ = EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache);
+        stashedDofIdx_ = -1;
+        focusDofIdx_ = -1;
+    }
+
+    static void *operator new(size_t size)
+    { return Opm::aligned_alloc(alignof(FvBaseElementContext), size); }
+
+    static void operator delete(void *ptr)
+    { Opm::aligned_free(ptr); }
+
+    /*!
+     * \brief Construct all volume and extensive quantities of an element
+     *        from scratch.
+     *
+     * \param elem The DUNE Codim<0> entity for which the volume
+     *             variables ought to be calculated
+     */
+    void updateAll(const Element& elem)
+    {
+        asImp_().updateStencil(elem);
+        asImp_().updateAllIntensiveQuantities();
+        asImp_().updateAllExtensiveQuantities();
+    }
+
+    /*!
+     * \brief Compute the finite volume geometry for an element.
+     *
+     * \param elem The grid element for which the finite volume geometry ought to be
+     *             computed.
+     */
+    void updateStencil(const Element& elem)
+    {
+        // remember the current element
+        elemPtr_ = &elem;
+
+        // update the stencil. the center gradients are quite expensive to calculate and
+        // most models don't need them, so that we only do this if the model explicitly
+        // enables them
+        stencil_.update(elem);
+
+        // resize the arrays containing the flux and the volume variables
+        dofVars_.resize(stencil_.numDof());
+        extensiveQuantities_.resize(stencil_.numInteriorFaces());
+    }
+
+    /*!
+     * \brief Update the primary topological part of the stencil, but nothing else.
+     *
+     * \param elem The grid element for which the finite volume geometry ought to be
+     *             computed.
+     */
+    void updatePrimaryStencil(const Element& elem)
+    {
+        // remember the current element
+        elemPtr_ = &elem;
+
+        // update the finite element geometry
+        stencil_.updatePrimaryTopology(elem);
+
+        dofVars_.resize(stencil_.numPrimaryDof());
+    }
+
+    /*!
+     * \brief Update the topological part of the stencil, but nothing else.
+     *
+     * \param elem The grid element for which the finite volume geometry ought to be
+     *             computed.
+     */
+    void updateStencilTopology(const Element& elem)
+    {
+        // remember the current element
+        elemPtr_ = &elem;
+
+        // update the finite element geometry
+        stencil_.updateTopology(elem);
+    }
+
+    /*!
+     * \brief Compute the intensive quantities of all sub-control volumes of the current
+     *        element for all time indices.
+     */
+    void updateAllIntensiveQuantities()
+    {
+        if (!enableStorageCache_) {
+            // if the storage cache is disabled, we need to calculate the storage term
+            // from scratch, i.e. we need the intensive quantities of all of the history.
+            for (unsigned timeIdx = 0; timeIdx < timeDiscHistorySize; ++ timeIdx)
+                asImp_().updateIntensiveQuantities(timeIdx);
+        }
+        else
+            // if the storage cache is enabled, we only need to recalculate the storage
+            // term for the most recent point of history (i.e., for the current iterative
+            // solution)
+            asImp_().updateIntensiveQuantities(/*timeIdx=*/0);
+    }
+
+    /*!
+     * \brief Compute the intensive quantities of all sub-control volumes of the current
+     *        element for a single time index.
+     *
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    void updateIntensiveQuantities(unsigned timeIdx)
+    { updateIntensiveQuantities_(timeIdx, numDof(timeIdx)); }
+
+    /*!
+     * \brief Compute the intensive quantities of all sub-control volumes of the current
+     *        element for a single time index.
+     *
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    void updatePrimaryIntensiveQuantities(unsigned timeIdx)
+    { updateIntensiveQuantities_(timeIdx, numPrimaryDof(timeIdx)); }
+
+    /*!
+     * \brief Compute the intensive quantities of a single sub-control volume of the
+     *        current element for a single time index.
+     *
+     * \param priVars The PrimaryVariables which should be used to calculate the
+     *                intensive quantities.
+     * \param dofIdx The local index in the current element of the sub-control volume
+     *               which should be updated.
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    void updateIntensiveQuantities(const PrimaryVariables& priVars, unsigned dofIdx, unsigned timeIdx)
+    { asImp_().updateSingleIntQuants_(priVars, dofIdx, timeIdx); }
+
+    /*!
+     * \brief Compute the extensive quantities of all sub-control volume
+     *        faces of the current element for all time indices.
+     */
+    void updateAllExtensiveQuantities()
+    { asImp_().updateExtensiveQuantities(/*timeIdx=*/0); }
+
+    /*!
+     * \brief Compute the extensive quantities of all sub-control volume
+     *        faces of the current element for a single time index.
+     *
+     * \param timeIdx The index of the solution vector used by the
+     *                time discretization.
+     */
+    void updateExtensiveQuantities(unsigned timeIdx)
+    {
+        gradientCalculator_.prepare(/*context=*/asImp_(), timeIdx);
+
+        for (unsigned fluxIdx = 0; fluxIdx < numInteriorFaces(timeIdx); fluxIdx++) {
+            extensiveQuantities_[fluxIdx].update(/*context=*/asImp_(),
+                                                 /*localIndex=*/fluxIdx,
+                                                 timeIdx);
+        }
+    }
+
+    /*!
+     * \brief Sets the degree of freedom on which the simulator is currently "focused" on
+     *
+     * I.e., in the case of automatic differentiation, all derivatives are with regard to
+     * the primary variables of that degree of freedom. Only "primary" DOFs can be
+     * focused on.
+     */
+    void setFocusDofIndex(unsigned dofIdx)
+    { focusDofIdx_ = dofIdx; }
+
+    /*!
+     * \brief Returns the degree of freedom on which the simulator is currently "focused" on
+     *
+     * \copydetails setFocusDof()
+     */
+    unsigned focusDofIndex() const
+    { return focusDofIdx_; }
+
+    /*!
+     * \brief Return a reference to the simulator.
+     */
+    const Simulator& simulator() const
+    { return *simulatorPtr_; }
+
+    /*!
+     * \brief Return a reference to the problem.
+     */
+    const Problem& problem() const
+    { return simulatorPtr_->problem(); }
+
+    /*!
+     * \brief Return a reference to the model.
+     */
+    const Model& model() const
+    { return simulatorPtr_->model(); }
+
+    /*!
+     * \brief Return a reference to the grid view.
+     */
+    const GridView& gridView() const
+    { return gridView_; }
+
+    /*!
+     * \brief Return the current element.
+     */
+    const Element& element() const
+    { return *elemPtr_; }
+
+    /*!
+     * \brief Return the number of sub-control volumes of the current element.
+     */
+    size_t numDof(unsigned timeIdx) const
+    { return stencil(timeIdx).numDof(); }
+
+    /*!
+     * \brief Return the number of primary degrees of freedom of the current element.
+     */
+    size_t numPrimaryDof(unsigned timeIdx) const
+    { return stencil(timeIdx).numPrimaryDof(); }
+
+    /*!
+     * \brief Return the number of non-boundary faces which need to be
+     *        considered for the flux apporixmation.
+     */
+    size_t numInteriorFaces(unsigned timeIdx) const
+    { return stencil(timeIdx).numInteriorFaces(); }
+
+    /*!
+     * \brief Return the number of boundary faces which need to be
+     *        considered for the flux apporixmation.
+     */
+    size_t numBoundaryFaces(unsigned timeIdx) const
+    { return stencil(timeIdx).numBoundaryFaces(); }
+
+    /*!
+     * \brief Return the current finite element geometry.
+     *
+     * \param timeIdx The index of the solution vector used by the
+     *                time discretization.
+     */
+    const Stencil& stencil(unsigned timeIdx OPM_UNUSED) const
+    { return stencil_; }
+
+    /*!
+     * \brief Return the position of a local entities in global coordinates
+     *
+     * \param dofIdx The local index of the degree of freedom
+     *               in the current element.
+     * \param timeIdx The index of the solution vector used by the
+     *                time discretization.
+     */
+    const GlobalPosition& pos(unsigned dofIdx, unsigned timeIdx OPM_UNUSED) const
+    { return stencil_.subControlVolume(dofIdx).globalPos(); }
+
+    /*!
+     * \brief Return the global spatial index for a sub-control volume
+     *
+     * \param dofIdx The local index of the degree of freedom
+     *               in the current element.
+     * \param timeIdx The index of the solution vector used by the
+     *                time discretization.
+     */
+    unsigned globalSpaceIndex(unsigned dofIdx, unsigned timeIdx) const
+    { return stencil(timeIdx).globalSpaceIndex(dofIdx); }
+
+
+    /*!
+     * \brief Return the element-local volume associated with a degree of freedom
+     *
+     * In the case of the vertex-centered finite volume method, this is different from
+     * the total volume because a finite volume usually spans multiple elements...
+     *
+     * \param dofIdx The local index of the degree of freedom in the current element.
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    Scalar dofVolume(unsigned dofIdx, unsigned timeIdx) const
+    { return stencil(timeIdx).subControlVolume(dofIdx).volume();  }
+
+    /*!
+     * \brief Return the total volume associated with a degree of freedom
+     *
+     * "Total" means the volume controlled by a degree of freedom disregarding the
+     * element. (For example in the vertex-centered finite volume method, a control
+     * volume typically encompasses parts of multiple elements.)
+     *
+     * \param dofIdx The local index of the degree of freedom in the current element.
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    Scalar dofTotalVolume(unsigned dofIdx, unsigned timeIdx) const
+    { return model().dofTotalVolume(globalSpaceIndex(dofIdx, timeIdx)); }
+
+    /*!
+     * \brief Returns whether the current element is on the domain's
+     *        boundary.
+     */
+    bool onBoundary() const
+    { return element().hasBoundaryIntersections(); }
+
+    /*!
+     * \brief Return a reference to the intensive quantities of a
+     *        sub-control volume at a given time.
+     *
+     * If the time step index is not given, return the volume
+     * variables for the current time.
+     *
+     * \param dofIdx The local index of the degree of freedom in the current element.
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    const IntensiveQuantities& intensiveQuantities(unsigned dofIdx, unsigned timeIdx) const
+    {
+#ifndef NDEBUG
+        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
+
+        if (enableStorageCache_ && timeIdx != 0 && problem().recycleFirstIterationStorage())
+            throw std::logic_error("If caching of the storage term is enabled, only the intensive quantities "
+                                   "for the most-recent substep (i.e. time index 0) are available!");
+#endif
+
+        return dofVars_[dofIdx].intensiveQuantities[timeIdx];
+    }
+
+    /*!
+     * \brief Return the thermodynamic hint for a given local index.
+     *
+     * \sa Discretization::thermodynamicHint(int, int)
+     *
+     * \param dofIdx The local index of the degree of freedom in the current element.
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    const IntensiveQuantities *thermodynamicHint(unsigned dofIdx, unsigned timeIdx) const
+    {
+        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
+        return dofVars_[dofIdx].thermodynamicHint[timeIdx];
+    }
+    /*!
+     * \copydoc intensiveQuantities()
+     */
+    IntensiveQuantities& intensiveQuantities(unsigned dofIdx, unsigned timeIdx)
+    {
+        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
+        return dofVars_[dofIdx].intensiveQuantities[timeIdx];
+    }
+
+    /*!
+     * \brief Return the primary variables for a given local index.
+     *
+     * \param dofIdx The local index of the degree of freedom
+     *               in the current element.
+     * \param timeIdx The index of the solution vector used by the
+     *                time discretization.
+     */
+    PrimaryVariables& primaryVars(unsigned dofIdx, unsigned timeIdx)
+    {
+        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
+        return dofVars_[dofIdx].priVars[timeIdx];
+    }
+    /*!
+     * \copydoc primaryVars()
+     */
+    const PrimaryVariables& primaryVars(unsigned dofIdx, unsigned timeIdx) const
+    {
+        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
+        return dofVars_[dofIdx].priVars[timeIdx];
+    }
+
+    /*!
+     * \brief Returns true if no intensive quanties are stashed
+     *
+     * In most cases quantities are stashed only if a partial derivative is to be
+     * calculated via finite difference methods.
+     */
+    bool haveStashedIntensiveQuantities() const
+    { return stashedDofIdx_ != -1; }
+
+    /*!
+     * \brief Return the (local) index of the DOF for which the primary variables were
+     *        stashed
+     *
+     * If none, then this returns -1.
+     */
+    int stashedDofIdx() const
+    { return stashedDofIdx_; }
+
+    /*!
+     * \brief Stash the intensive quantities for a degree of freedom on internal memory.
+     *
+     * \param dofIdx The local index of the degree of freedom in the current element.
+     */
+    void stashIntensiveQuantities(unsigned dofIdx)
+    {
+        assert(0 <= dofIdx && dofIdx < numDof(/*timeIdx=*/0));
+
+        intensiveQuantitiesStashed_ = dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0];
+        priVarsStashed_ = dofVars_[dofIdx].priVars[/*timeIdx=*/0];
+        stashedDofIdx_ = static_cast<int>(dofIdx);
+    }
+
+    /*!
+     * \brief Restores the intensive quantities for a degree of freedom from internal memory.
+     *
+     * \param dofIdx The local index of the degree of freedom in the current element.
+     */
+    void restoreIntensiveQuantities(unsigned dofIdx)
+    {
+        dofVars_[dofIdx].priVars[/*timeIdx=*/0] = priVarsStashed_;
+        dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0] = intensiveQuantitiesStashed_;
+        stashedDofIdx_ = -1;
+    }
+
+    /*!
+     * \brief Return a reference to the gradient calculation class of
+     *        the chosen spatial discretization.
+     */
+    const GradientCalculator& gradientCalculator() const
+    { return gradientCalculator_; }
+
+    /*!
+     * \brief Return a reference to the extensive quantities of a
+     *        sub-control volume face.
+     *
+     * \param fluxIdx The local index of the sub-control volume face for which the
+     *               extensive quantities are requested
+     * \param timeIdx The index of the solution vector used by the time discretization.
+     */
+    const ExtensiveQuantities& extensiveQuantities(unsigned fluxIdx, unsigned timeIdx OPM_UNUSED) const
+    { return extensiveQuantities_[fluxIdx]; }
+
+    /*!
+     * \brief Returns true iff the cache for the storage term ought to be used for this context.
+     *
+     * If it is used, intensive quantities can only be accessed for the most recent time
+     * index. (time index 0.)
+     */
+    bool enableStorageCache() const
+    { return enableStorageCache_; }
+
+    /*!
+     * \brief Specifies if the cache for the storage term ought to be used for this context.
+     */
+    void setEnableStorageCache(bool yesno)
+    { enableStorageCache_ = yesno; }
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+
+protected:
+    /*!
+     * \brief Update the first 'n' intensive quantities objects from the primary variables.
+     *
+     * This method considers the intensive quantities cache.
+     */
+    void updateIntensiveQuantities_(unsigned timeIdx, size_t numDof)
+    {
+        // update the intensive quantities for the whole history
+        const SolutionVector& globalSol = model().solution(timeIdx);
+
+        // update the non-gradient quantities
+        for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
+            unsigned globalIdx = globalSpaceIndex(dofIdx, timeIdx);
+            const PrimaryVariables& dofSol = globalSol[globalIdx];
+            dofVars_[dofIdx].priVars[timeIdx] = dofSol;
+
+            dofVars_[dofIdx].thermodynamicHint[timeIdx] =
+                model().thermodynamicHint(globalIdx, timeIdx);
+
+            const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
+            if (cachedIntQuants) {
+                dofVars_[dofIdx].intensiveQuantities[timeIdx] = *cachedIntQuants;
+            }
+            else {
+                updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
+                model().updateCachedIntensiveQuantities(dofVars_[dofIdx].intensiveQuantities[timeIdx],
+                                                        globalIdx,
+                                                        timeIdx);
+            }
+        }
+    }
+
+    void updateSingleIntQuants_(const PrimaryVariables& priVars, unsigned dofIdx, unsigned timeIdx)
+    {
+#ifndef NDEBUG
+        if (enableStorageCache_ && timeIdx != 0 && problem().recycleFirstIterationStorage())
+            throw std::logic_error("If caching of the storage term is enabled, only the intensive quantities "
+                                   "for the most-recent substep (i.e. time index 0) are available!");
+#endif
+
+        dofVars_[dofIdx].priVars[timeIdx] = priVars;
+        dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/asImp_(), dofIdx, timeIdx);
+    }
+
+    IntensiveQuantities intensiveQuantitiesStashed_;
+    PrimaryVariables priVarsStashed_;
+
+    GradientCalculator gradientCalculator_;
+
+    std::vector<DofStore_, Opm::aligned_allocator<DofStore_, alignof(DofStore_)> > dofVars_;
+    std::vector<ExtensiveQuantities, Opm::aligned_allocator<ExtensiveQuantities, alignof(ExtensiveQuantities)> > extensiveQuantities_;
+
+    const Simulator *simulatorPtr_;
+    const Element *elemPtr_;
+    const GridView gridView_;
+    Stencil stencil_;
+
+    int stashedDofIdx_;
+    int focusDofIdx_;
+    bool enableStorageCache_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseextensivequantities.hh

@@ -0,0 +1,134 @@
+// -*- 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::FvBaseExtensiveQuantities
+ */
+#ifndef EWOMS_FV_BASE_EXTENSIVE_QUANTITIES_HH
+#define EWOMS_FV_BASE_EXTENSIVE_QUANTITIES_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+
+namespace Opm {
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Provide the properties at a face which make sense indepentently
+ *        of the conserved quantities.
+ */
+template <class TypeTag>
+class FvBaseExtensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+    /*!
+     * \brief Register all run-time parameters for the extensive quantities.
+     */
+    static void registerParameters()
+    { }
+
+    /*!
+     * \brief Update the extensive quantities for a given sub-control-volume face.
+     *
+     * \param elemCtx Reference to the current element context.
+     * \param scvfIdx The local index of the sub-control-volume face for which the
+     *                extensive quantities should be calculated.
+     * \param timeIdx The index used by the time discretization.
+     */
+    void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
+    {
+        const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
+        interiorScvIdx_ = scvf.interiorIndex();
+        exteriorScvIdx_ = scvf.exteriorIndex();
+
+        extrusionFactor_ =
+            (elemCtx.intensiveQuantities(interiorScvIdx_, timeIdx).extrusionFactor()
+             + elemCtx.intensiveQuantities(exteriorScvIdx_, timeIdx).extrusionFactor()) / 2;
+        Opm::Valgrind::CheckDefined(extrusionFactor_);
+        assert(extrusionFactor_ > 0);
+    }
+
+
+    /*!
+     * \brief Update the extensive quantities for a given boundary face.
+     *
+     * \param context Reference to the current execution context.
+     * \param bfIdx The local index of the boundary face for which
+     *              the extensive quantities should be calculated.
+     * \param timeIdx The index used by the time discretization.
+     * \param fluidState The FluidState on the domain boundary.
+     * \param paramCache The FluidSystem's parameter cache.
+     */
+    template <class Context, class FluidState>
+    void updateBoundary(const Context& context,
+                        unsigned bfIdx,
+                        unsigned timeIdx,
+                        const FluidState& fluidState OPM_UNUSED)
+    {
+        unsigned dofIdx = context.interiorScvIndex(bfIdx, timeIdx);
+        interiorScvIdx_ = static_cast<unsigned short>(dofIdx);
+        exteriorScvIdx_ = static_cast<unsigned short>(dofIdx);
+
+        extrusionFactor_ = context.intensiveQuantities(bfIdx, timeIdx).extrusionFactor();
+        Opm::Valgrind::CheckDefined(extrusionFactor_);
+        assert(extrusionFactor_ > 0);
+    }
+
+    /*!
+     * \brief Returns th extrusion factor for the sub-control-volume face
+     */
+    Scalar extrusionFactor() const
+    { return extrusionFactor_; }
+
+    /*!
+     * \brief Return the local index of the control volume which is on the "interior" of
+     *        the sub-control volume face.
+     */
+    unsigned short interiorIndex() const
+    { return interiorScvIdx_; }
+
+    /*!
+     * \brief Return the local index of the control volume which is on the "exterior" of
+     *        the sub-control volume face.
+     */
+    unsigned short exteriorIndex() const
+    { return exteriorScvIdx_; }
+
+private:
+    // local indices of the interior and the exterior sub-control-volumes
+    unsigned short interiorScvIdx_;
+    unsigned short exteriorScvIdx_;
+
+    Scalar extrusionFactor_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbasefdlocallinearizer.hh

@@ -0,0 +1,520 @@
+// -*- 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::FvBaseFdLocalLinearizer
+ */
+#ifndef EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
+#define EWOMS_FV_BASE_FD_LOCAL_LINEARIZER_HH
+
+#include <opm/models/utils/propertysystem.hh>
+#include <opm/models/utils/parametersystem.hh>
+
+#include <opm/material/common/MathToolbox.hpp>
+#include <opm/material/common/Valgrind.hpp>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/matrix.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <limits>
+
+namespace Opm {
+// forward declaration
+template<class TypeTag>
+class FvBaseFdLocalLinearizer;
+
+} // namespace Opm
+
+BEGIN_PROPERTIES
+
+// declare the property tags required for the finite differences local linearizer
+NEW_TYPE_TAG(FiniteDifferenceLocalLinearizer);
+
+NEW_PROP_TAG(LocalLinearizer);
+NEW_PROP_TAG(Evaluation);
+NEW_PROP_TAG(NumericDifferenceMethod);
+NEW_PROP_TAG(BaseEpsilon);
+NEW_PROP_TAG(SparseMatrixAdapter);
+NEW_PROP_TAG(LocalResidual);
+NEW_PROP_TAG(Simulator);
+NEW_PROP_TAG(Problem);
+NEW_PROP_TAG(Model);
+NEW_PROP_TAG(PrimaryVariables);
+NEW_PROP_TAG(ElementContext);
+NEW_PROP_TAG(Scalar);
+NEW_PROP_TAG(Evaluation);
+NEW_PROP_TAG(GridView);
+NEW_PROP_TAG(NumEq);
+
+// set the properties to be spliced in
+SET_TYPE_PROP(FiniteDifferenceLocalLinearizer, LocalLinearizer,
+              Opm::FvBaseFdLocalLinearizer<TypeTag>);
+
+SET_TYPE_PROP(FiniteDifferenceLocalLinearizer, Evaluation,
+              typename GET_PROP_TYPE(TypeTag, Scalar));
+
+/*!
+ * \brief Specify which kind of method should be used to numerically
+ * calculate the partial derivatives of the residual.
+ *
+ * -1 means backward differences, 0 means central differences, 1 means
+ * forward differences. By default we use central differences.
+ */
+SET_INT_PROP(FiniteDifferenceLocalLinearizer, NumericDifferenceMethod, +1);
+
+//! The base epsilon value for finite difference calculations
+SET_SCALAR_PROP(FiniteDifferenceLocalLinearizer,
+                BaseEpsilon,
+                std::max<Scalar>(0.9123e-10, std::numeric_limits<Scalar>::epsilon()*1.23e3));
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Calculates the Jacobian of the local residual for finite volume spatial
+ *        discretizations using a finite difference method
+ *
+ * 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.
+ *
+ * This class implements numeric differentiation using finite difference methods, i.e.
+ * forward or backward differences (2nd order), or central differences (3rd order). The
+ * method used is determined by the "NumericDifferenceMethod" property:
+ *
+ * - If the value of this property is smaller than 0, backward differences are used,
+ *   i.e.:
+ *   \f[
+ *     \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
+ *   \f]
+ *
+ * - If the value of this property is 0, central differences are used, i.e.:
+ *   \f[
+ *     \frac{\partial f(x)}{\partial x} \approx
+ *          \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
+ *   \f]
+ *
+ * - if the value of this property is larger than 0, forward differences are used, i.e.:
+ *   \f[
+ *     \frac{\partial f(x)}{\partial x} \approx
+ *          \frac{f(x + \epsilon) - f(x)}{\epsilon}
+ *   \f]
+ *
+ * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$ is the value of a
+ * sub-control volume's primary variable at the evaluation point and \f$\epsilon\f$ is a
+ * small scalar value larger than 0.
+ */
+template<class TypeTag>
+class FvBaseFdLocalLinearizer
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, LocalLinearizer) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) LocalResidual;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+    // extract local matrices from jacobian matrix for consistency
+    typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter)::MatrixBlock ScalarMatrixBlock;
+    typedef Dune::FieldVector<Scalar, numEq> ScalarVectorBlock;
+
+    typedef Dune::BlockVector<ScalarVectorBlock> ScalarLocalBlockVector;
+    typedef Dune::Matrix<ScalarMatrixBlock> ScalarLocalBlockMatrix;
+
+    typedef typename LocalResidual::LocalEvalBlockVector LocalEvalBlockVector;
+
+#if __GNUC__ == 4 && __GNUC_MINOR__ <= 6
+public:
+    // make older GCCs happy by providing a public copy constructor (this is necessary
+    // for their implementation of std::vector, although the method is never called...)
+    FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&)
+        : internalElemContext_(0)
+    {}
+
+#else
+    // copying local residual objects around is a very bad idea, so we explicitly prevent
+    // it...
+    FvBaseFdLocalLinearizer(const FvBaseFdLocalLinearizer&) = delete;
+#endif
+public:
+    FvBaseFdLocalLinearizer()
+        : internalElemContext_(0)
+    { }
+
+    ~FvBaseFdLocalLinearizer()
+    { delete internalElemContext_; }
+
+    /*!
+     * \brief Register all run-time parameters for the local jacobian.
+     */
+    static void registerParameters()
+    {
+        EWOMS_REGISTER_PARAM(TypeTag, int, NumericDifferenceMethod,
+                             "The method used for numeric differentiation (-1: backward "
+                             "differences, 0: central differences, 1: forward differences)");
+    }
+
+    /*!
+     * \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.updateAll(elem);
+
+        // update the weights of the primary variables for the context
+        model_().updatePVWeights(elemCtx);
+
+        resize_(elemCtx);
+        reset_(elemCtx);
+
+        // calculate the local residual
+        localResidual_.eval(residual_, elemCtx);
+
+        // calculate the local jacobian matrix
+        size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+        for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; dofIdx++) {
+            for (unsigned pvIdx = 0; pvIdx < numEq; pvIdx++) {
+                asImp_().evalPartialDerivative_(elemCtx, dofIdx, pvIdx);
+
+                // incorporate the partial derivatives into the local Jacobian matrix
+                updateLocalJacobian_(elemCtx, dofIdx, pvIdx);
+            }
+        }
+    }
+
+    /*!
+     * \brief Returns the unweighted epsilon value used to calculate
+     *        the local derivatives
+     */
+    static Scalar baseEpsilon()
+    { return GET_PROP_VALUE(TypeTag, BaseEpsilon); }
+
+    /*!
+     * \brief Returns the epsilon value which is added and removed
+     *        from the current solution.
+     *
+     * \param elemCtx The element execution context for which the
+     *                local residual and its gradient should be
+     *                calculated.
+     * \param dofIdx     The local index of the element's vertex for
+     *                   which the local derivative ought to be calculated.
+     * \param pvIdx      The index of the primary variable which gets varied
+     */
+    Scalar numericEpsilon(const ElementContext& elemCtx,
+                          unsigned dofIdx,
+                          unsigned pvIdx) const
+    {
+        unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
+        Scalar pvWeight = elemCtx.model().primaryVarWeight(globalIdx, pvIdx);
+        assert(pvWeight > 0 && std::isfinite(pvWeight));
+        Opm::Valgrind::CheckDefined(pvWeight);
+
+        return baseEpsilon()/pvWeight;
+    }
+
+    /*!
+     * \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
+     *                     which contains the independents
+     * \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 Returns the numeric difference method which is applied.
+     */
+    static int numericDifferenceMethod_()
+    { return EWOMS_GET_PARAM(TypeTag, int, NumericDifferenceMethod); }
+
+    /*!
+     * \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);
+        jacobian_.setSize(numDof, numPrimaryDof);
+
+        derivResidual_.resize(numDof);
+    }
+
+    /*!
+     * \brief Reset the all relevant internal attributes to 0
+     */
+    void reset_(const ElementContext& elemCtx)
+    {
+        size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+        size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
+        for (unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx)
+            for (unsigned dof2Idx = 0; dof2Idx < numDof; ++ dof2Idx)
+                jacobian_[dof2Idx][primaryDofIdx] = 0.0;
+
+        for (unsigned primaryDofIdx = 0; primaryDofIdx < numDof; ++ primaryDofIdx)
+            residual_[primaryDofIdx] = 0.0;
+    }
+
+    /*!
+     * \brief Compute the partial derivatives of a context's residual functions
+     *
+     * This method can be overwritten by the implementation if a better scheme than
+     * numerical differentiation is available.
+     *
+     * The default implementation of this method uses numeric differentiation,
+     * i.e. forward or backward differences (2nd order), or central differences (3rd
+     * order). The method used is determined by the "NumericDifferenceMethod" property:
+     *
+     * - If the value of this property is smaller than 0, backward differences are used,
+     *   i.e.:
+     *   \f[
+     *     \frac{\partial f(x)}{\partial x} \approx \frac{f(x) - f(x - \epsilon)}{\epsilon}
+     *   \f]
+     *
+     * - If the value of this property is 0, central differences are used, i.e.:
+     *   \f[
+     *     \frac{\partial f(x)}{\partial x} \approx
+     *          \frac{f(x + \epsilon) - f(x - \epsilon)}{2 \epsilon}
+     *   \f]
+     *
+     * - if the value of this property is larger than 0, forward
+     *   differences are used, i.e.:
+     *   \f[
+           \frac{\partial f(x)}{\partial x} \approx \frac{f(x + \epsilon) - f(x)}{\epsilon}
+     *   \f]
+     *
+     * Here, \f$ f \f$ is the residual function for all equations, \f$x\f$ is the value
+     * of a sub-control volume's primary variable at the evaluation point and
+     * \f$\epsilon\f$ is a small value larger than 0.
+     *
+     * \param elemCtx The element context for which the local partial
+     *                derivative ought to be calculated
+     * \param dofIdx The sub-control volume index of the current
+     *               finite element for which the partial derivative
+     *               ought to be calculated
+     * \param pvIdx The index of the primary variable at the dofIdx'
+     *              sub-control volume of the current finite element
+     *              for which the partial derivative ought to be
+     *              calculated
+     */
+    void evalPartialDerivative_(ElementContext& elemCtx,
+                                unsigned dofIdx,
+                                unsigned pvIdx)
+    {
+        // save all quantities which depend on the specified primary
+        // variable at the given sub control volume
+        elemCtx.stashIntensiveQuantities(dofIdx);
+
+        PrimaryVariables priVars(elemCtx.primaryVars(dofIdx, /*timeIdx=*/0));
+        Scalar eps = asImp_().numericEpsilon(elemCtx, dofIdx, pvIdx);
+        Scalar delta = 0.0;
+
+        if (numericDifferenceMethod_() >= 0) {
+            // we are not using backward differences, i.e. we need to
+            // calculate f(x + \epsilon)
+
+            // deflect primary variables
+            priVars[pvIdx] += eps;
+            delta += eps;
+
+            // calculate the deflected residual
+            elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
+            elemCtx.updateAllExtensiveQuantities();
+            localResidual_.eval(derivResidual_, elemCtx);
+        }
+        else {
+            // we are using backward differences, i.e. we don't need
+            // to calculate f(x + \epsilon) and we can recycle the
+            // (already calculated) residual f(x)
+            derivResidual_ = residual_;
+        }
+
+        if (numericDifferenceMethod_() <= 0) {
+            // we are not using forward differences, i.e. we don't
+            // need to calculate f(x - \epsilon)
+
+            // deflect the primary variables
+            priVars[pvIdx] -= delta + eps;
+            delta += eps;
+
+            // calculate the deflected residual again, this time we use the local
+            // residual's internal storage.
+            elemCtx.updateIntensiveQuantities(priVars, dofIdx, /*timeIdx=*/0);
+            elemCtx.updateAllExtensiveQuantities();
+            localResidual_.eval(elemCtx);
+
+            derivResidual_ -= localResidual_.residual();
+        }
+        else {
+            // we are using forward differences, i.e. we don't need to
+            // calculate f(x - \epsilon) and we can recycle the
+            // (already calculated) residual f(x)
+            derivResidual_ -= residual_;
+        }
+
+        assert(delta > 0);
+
+        // divide difference in residuals by the magnitude of the
+        // deflections between the two function evaluation
+        derivResidual_ /= delta;
+
+        // restore the original state of the element's volume
+        // variables
+        elemCtx.restoreIntensiveQuantities(dofIdx);
+
+#ifndef NDEBUG
+        for (unsigned i = 0; i < derivResidual_.size(); ++i)
+            Opm::Valgrind::CheckDefined(derivResidual_[i]);
+#endif
+    }
+
+    /*!
+     * \brief Updates the current local Jacobian matrix with the partial derivatives of
+     *        all equations for primary variable 'pvIdx' at the degree of freedom
+     *        associated with 'focusDofIdx'.
+     */
+    void updateLocalJacobian_(const ElementContext& elemCtx,
+                              unsigned focusDofIdx,
+                              unsigned pvIdx)
+    {
+        size_t numDof = elemCtx.numDof(/*timeIdx=*/0);
+        for (unsigned dofIdx = 0; dofIdx < numDof; dofIdx++) {
+            for (unsigned eqIdx = 0; eqIdx < numEq; eqIdx++) {
+                // A[dofIdx][focusDofIdx][eqIdx][pvIdx] is the partial derivative of the
+                // residual function 'eqIdx' for the degree of freedom 'dofIdx' with
+                // regard to the primary variable 'pvIdx' of the degree of freedom
+                // 'focusDofIdx'
+                jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx] = derivResidual_[dofIdx][eqIdx];
+                Opm::Valgrind::CheckDefined(jacobian_[dofIdx][focusDofIdx][eqIdx][pvIdx]);
+            }
+        }
+    }
+
+    Simulator *simulatorPtr_;
+    Model *modelPtr_;
+
+    ElementContext *internalElemContext_;
+
+    LocalEvalBlockVector residual_;
+    LocalEvalBlockVector derivResidual_;
+    ScalarLocalBlockMatrix jacobian_;
+
+    LocalResidual localResidual_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbasegradientcalculator.hh

@@ -0,0 +1,358 @@
+// -*- 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::FvBaseGradientCalculator
+ */
+#ifndef EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH
+#define EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/material/common/Unused.hpp>
+
+#include <dune/common/fvector.hh>
+
+namespace Opm {
+template<class TypeTag>
+class EcfvDiscretization;
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief This class calculates gradients of arbitrary quantities at
+ *        flux integration points using the two-point approximation scheme
+ */
+template<class TypeTag>
+class FvBaseGradientCalculator
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    enum { dim = GridView::dimension };
+    enum { dimWorld = GridView::dimensionworld };
+
+    // maximum number of flux approximation points. to calculate this,
+    // we assume that the geometry with the most pointsq is a cube.
+    enum { maxFap = 2 << dim };
+
+    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
+    typedef Dune::FieldVector<Evaluation, dimWorld> EvalDimVector;
+
+public:
+    /*!
+     * \brief Register all run-time parameters for the gradient calculator
+     *        of the base class of the discretization.
+     */
+    static void registerParameters()
+    { }
+
+    /*!
+     * \brief Precomputes the common values to calculate gradients and values of
+     *        quantities at every interior flux approximation point.
+     *
+     * \param elemCtx The current execution context
+     * \param timeIdx The index used by the time discretization.
+     */
+    template <bool prepareValues = true, bool prepareGradients = true>
+    void prepare(const ElementContext& elemCtx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
+    { /* noting to do */ }
+
+    /*!
+     * \brief Calculates the value of an arbitrary scalar quantity at any interior flux
+     *        approximation point.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point in the current
+     *               element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    auto calculateScalarValue(const ElementContext& elemCtx,
+                              unsigned fapIdx,
+                              const QuantityCallback& quantityCallback) const
+        -> typename std::remove_reference<decltype(quantityCallback.operator()(0))>::type
+    {
+        typedef decltype(quantityCallback.operator()(0)) RawReturnType;
+        typedef typename std::remove_const<typename std::remove_reference<RawReturnType>::type>::type ReturnType;
+
+        Scalar interiorDistance;
+        Scalar exteriorDistance;
+        computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);
+
+        const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(fapIdx);
+        auto i = face.interiorIndex();
+        auto j = face.exteriorIndex();
+        auto focusDofIdx = elemCtx.focusDofIndex();
+
+        // use the average weighted by distance...
+        ReturnType value;
+        if (i == focusDofIdx)
+            value = quantityCallback(i)*interiorDistance;
+        else
+            value = Opm::getValue(quantityCallback(i))*interiorDistance;
+
+        if (j == focusDofIdx)
+            value += quantityCallback(j)*exteriorDistance;
+        else
+            value += Opm::getValue(quantityCallback(j))*exteriorDistance;
+
+        value /= interiorDistance + exteriorDistance;
+
+        return value;
+    }
+
+    /*!
+     * \brief Calculates the value of an arbitrary vectorial quantity at any interior flux
+     *        approximation point.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point in the current
+     *               element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    auto calculateVectorValue(const ElementContext& elemCtx,
+                              unsigned fapIdx,
+                              const QuantityCallback& quantityCallback) const
+        -> typename std::remove_reference<decltype(quantityCallback.operator()(0))>::type
+    {
+        typedef decltype(quantityCallback.operator()(0)) RawReturnType;
+        typedef typename std::remove_const<typename std::remove_reference<RawReturnType>::type>::type ReturnType;
+
+        Scalar interiorDistance;
+        Scalar exteriorDistance;
+        computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);
+
+        const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(fapIdx);
+        auto i = face.interiorIndex();
+        auto j = face.exteriorIndex();
+        auto focusDofIdx = elemCtx.focusDofIndex();
+
+        // use the average weighted by distance...
+        ReturnType value;
+        if (i == focusDofIdx) {
+            value = quantityCallback(i);
+            for (int k = 0; k < value.size(); ++k)
+                value[k] *= interiorDistance;
+        }
+        else {
+            const auto& dofVal = Opm::getValue(quantityCallback(i));
+            for (int k = 0; k < dofVal.size(); ++k)
+                value[k] = Opm::getValue(dofVal[k])*interiorDistance;
+        }
+
+        if (j == focusDofIdx) {
+            const auto& dofVal = quantityCallback(j);
+            for (int k = 0; k < dofVal.size(); ++k)
+                value[k] += dofVal[k]*exteriorDistance;
+        }
+        else {
+            const auto& dofVal = quantityCallback(j);
+            for (int k = 0; k < dofVal.size(); ++k)
+                value[k] += Opm::getValue(dofVal[k])*exteriorDistance;
+        }
+
+        Scalar totDistance = interiorDistance + exteriorDistance;
+        for (int k = 0; k < value.size(); ++k)
+            value[k] /= totDistance;
+
+        return value;
+    }
+
+    /*!
+     * \brief Calculates the gradient of an arbitrary quantity at any
+     *        flux approximation point.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point
+     *               in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity given the index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    void calculateGradient(EvalDimVector& quantityGrad,
+                           const ElementContext& elemCtx,
+                           unsigned fapIdx,
+                           const QuantityCallback& quantityCallback) const
+    {
+        const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
+        const auto& face = stencil.interiorFace(fapIdx);
+
+        auto i = face.interiorIndex();
+        auto j = face.exteriorIndex();
+        auto focusIdx = elemCtx.focusDofIndex();
+
+        const auto& interiorPos = stencil.subControlVolume(i).globalPos();
+        const auto& exteriorPos = stencil.subControlVolume(j).globalPos();
+
+        Evaluation deltay;
+        if (i == focusIdx) {
+            deltay =
+                Opm::getValue(quantityCallback(j))
+                - quantityCallback(i);
+        }
+        else if (j == focusIdx) {
+            deltay =
+                quantityCallback(j)
+                - Opm::getValue(quantityCallback(i));
+        }
+        else
+            deltay =
+                Opm::getValue(quantityCallback(j))
+                - Opm::getValue(quantityCallback(i));
+
+        Scalar distSquared = 0.0;
+        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
+            Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
+            distSquared += tmp*tmp;
+        }
+
+        // divide the gradient by the squared distance between the centers of the
+        // sub-control volumes: the gradient is the normalized directional vector between
+        // the two centers times the ratio of the difference of the values and their
+        // distance, i.e., d/abs(d) * delta y / abs(d) = d*delta y / abs(d)^2.
+        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
+            Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
+            quantityGrad[dimIdx] = deltay*(tmp/distSquared);
+        }
+    }
+
+    /*!
+     * \brief Calculates the value of an arbitrary quantity at any
+     *        flux approximation point on the grid boundary.
+     *
+     * Boundary values are always calculated using the two-point
+     * approximation.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point
+     *               in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity given the index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    auto calculateBoundaryValue(const ElementContext& elemCtx OPM_UNUSED,
+                                unsigned fapIdx OPM_UNUSED,
+                                const QuantityCallback& quantityCallback)
+        -> decltype(quantityCallback.boundaryValue())
+    { return quantityCallback.boundaryValue(); }
+
+    /*!
+     * \brief Calculates the gradient of an arbitrary quantity at any
+     *        flux approximation point on the boundary.
+     *
+     * Boundary gradients are always calculated using the two-point
+     * approximation.
+     *
+     * \param elemCtx The current execution context
+     * \param faceIdx The local index of the flux approximation point
+     *                in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    void calculateBoundaryGradient(EvalDimVector& quantityGrad,
+                                   const ElementContext& elemCtx,
+                                   unsigned faceIdx,
+                                   const QuantityCallback& quantityCallback) const
+    {
+        const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
+        const auto& face = stencil.boundaryFace(faceIdx);
+
+        Evaluation deltay;
+        if (face.interiorIndex() == elemCtx.focusDofIndex())
+            deltay = quantityCallback.boundaryValue() - quantityCallback(face.interiorIndex());
+        else
+            deltay =
+                Opm::getValue(quantityCallback.boundaryValue())
+                - Opm::getValue(quantityCallback(face.interiorIndex()));
+
+        const auto& boundaryFacePos = face.integrationPos();
+        const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).center();
+
+        Scalar distSquared = 0;
+        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
+            Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
+            distSquared += tmp*tmp;
+        }
+
+        // divide the gradient by the squared distance between the center of the
+        // sub-control and the center of the boundary face: the gradient is the
+        // normalized directional vector between the two centers times the ratio of the
+        // difference of the values and their distance, i.e., d/abs(d) * deltay / abs(d)
+        // = d*deltay / abs(d)^2.
+        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
+            Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
+            quantityGrad[dimIdx] = deltay*(tmp/distSquared);
+        }
+    }
+
+private:
+    void computeDistances_(Scalar& interiorDistance,
+                           Scalar& exteriorDistance,
+                           const ElementContext& elemCtx,
+                           unsigned fapIdx) const
+    {
+        const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
+        const auto& face = stencil.interiorFace(fapIdx);
+
+        // calculate the distances of the position of the interior and of the exterior
+        // finite volume to the position of the integration point.
+        const auto& normal = face.normal();
+        auto i = face.interiorIndex();
+        auto j = face.exteriorIndex();
+        const auto& interiorPos = stencil.subControlVolume(i).globalPos();
+        const auto& exteriorPos = stencil.subControlVolume(j).globalPos();
+        const auto& integrationPos = face.integrationPos();
+
+        interiorDistance = 0.0;
+        exteriorDistance = 0.0;
+        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
+            interiorDistance +=
+                (interiorPos[dimIdx] - integrationPos[dimIdx])
+                * normal[dimIdx];
+
+            exteriorDistance +=
+                (exteriorPos[dimIdx] - integrationPos[dimIdx])
+                * normal[dimIdx];
+        }
+
+        interiorDistance = std::sqrt(std::abs(interiorDistance));
+        exteriorDistance = std::sqrt(std::abs(exteriorDistance));
+    }
+};
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseintensivequantities.hh

@@ -0,0 +1,103 @@
+// -*- 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::FvBaseIntensiveQuantities
+ */
+#ifndef EWOMS_FV_BASE_INTENSIVE_QUANTITIES_HH
+#define EWOMS_FV_BASE_INTENSIVE_QUANTITIES_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Base class for the model specific class which provides access to all intensive
+ *        (i.e., volume averaged) quantities.
+ */
+template <class TypeTag>
+class FvBaseIntensiveQuantities
+{
+    typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+public:
+    // default constructor
+    FvBaseIntensiveQuantities()
+    { }
+
+    // copy constructor
+    FvBaseIntensiveQuantities(const FvBaseIntensiveQuantities& v) = default;
+
+    /*!
+     * \brief Register all run-time parameters for the intensive quantities.
+     */
+    static void registerParameters()
+    { }
+
+    /*!
+     * \brief Update all quantities for a given control volume.
+     */
+    void update(const ElementContext& elemCtx,
+                unsigned dofIdx,
+                unsigned timeIdx)
+    { extrusionFactor_ = elemCtx.problem().extrusionFactor(elemCtx, dofIdx, timeIdx); }
+
+    /*!
+     * \brief Return how much a given sub-control volume is extruded.
+     *
+     * This means the factor by which a lower-dimensional (1D or 2D)
+     * entity needs to be expanded to get a full dimensional cell. The
+     * default is 1.0 which means that 1D problems are actually
+     * thought as pipes with a cross section of 1 m^2 and 2D problems
+     * are assumed to extend 1 m to the back.
+     */
+    Scalar extrusionFactor() const
+    { return extrusionFactor_; }
+
+    /*!
+     * \brief If running in valgrind this makes sure that all
+     *        quantities in the intensive quantities are defined.
+     */
+    void checkDefined() const
+    { }
+
+private:
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    Scalar extrusionFactor_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaselinearizer.hh

@@ -0,0 +1,601 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \copydoc Opm::FvBaseLinearizer
+ */
+#ifndef EWOMS_FV_BASE_LINEARIZER_HH
+#define EWOMS_FV_BASE_LINEARIZER_HH
+
+#include "fvbaseproperties.hh"
+
+#include <opm/models/parallel/gridcommhandles.hh>
+#include <opm/models/parallel/threadmanager.hh>
+#include <opm/models/parallel/threadedentityiterator.hh>
+#include <opm/models/discretization/common/baseauxiliarymodule.hh>
+
+#include <opm/material/common/Exceptions.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 FvBaseLinearizer
+{
+//! \cond SKIP_THIS
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, Discretization) Discretization;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, SparseMatrixAdapter) SparseMatrixAdapter;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
+    typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
+    typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
+
+    typedef typename GET_PROP_TYPE(TypeTag, GridCommHandleFactory) GridCommHandleFactory;
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+    typedef GlobalEqVector Vector;
+
+    typedef typename SparseMatrixAdapter::IstlMatrix IstlMatrix;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { historySize = GET_PROP_VALUE(TypeTag, TimeDiscHistorySize) };
+
+    typedef typename SparseMatrixAdapter::MatrixBlock MatrixBlock;
+    typedef Dune::FieldVector<Scalar, numEq> VectorBlock;
+
+    static const bool linearizeNonLocalElements = GET_PROP_VALUE(TypeTag, LinearizeNonLocalElements);
+
+    // copying the linearizer is not a good idea
+    FvBaseLinearizer(const FvBaseLinearizer&);
+//! \endcond
+
+public:
+    FvBaseLinearizer()
+        : jacobian_()
+    {
+        simulatorPtr_ = 0;
+    }
+
+    ~FvBaseLinearizer()
+    {
+        auto it = elementCtx_.begin();
+        const auto& endIt = elementCtx_.end();
+        for (; it != endIt; ++it)
+            delete *it;
+    }
+
+    /*!
+     * \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.
+     *
+     * This means the spatial domain plus 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;
+        }
+#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,5)
+        catch (const Dune::Exception& e)
+        {
+            std::cout << "rank " << simulator_().gridView().comm().rank()
+                      << " caught an exception while linearizing:" << e.what()
+                      << "\n"  << std::flush;
+            succeeded = 0;
+        }
+#endif
+        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 Opm::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;
+            }
+#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,5)
+            catch (const Dune::Exception& e)
+            {
+                succeeded = false;
+
+                std::cout << "rank " << simulator_().gridView().comm().rank()
+                          << " caught an exception while linearizing:" << e.what()
+                          << "\n"  << std::flush;
+            }
+#endif
+
+            succeeded = comm.min(succeeded);
+
+            if (!succeeded)
+                throw Opm::NumericalIssue("linearization of an auxilary 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_; }
+
+    /*!
+     * \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 constraintsMap_; }
+
+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(); }
+
+    const ElementMapper& elementMapper_() const
+    { return model_().elementMapper(); }
+
+    const DofMapper& dofMapper_() const
+    { return model_().dofMapper(); }
+
+    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_();
+
+        // create the per-thread context objects
+        elementCtx_.resize(ThreadManager::maxThreads());
+        for (unsigned threadId = 0; threadId != ThreadManager::maxThreads(); ++ threadId)
+            elementCtx_[threadId] = new ElementContext(simulator_());
+    }
+
+    // 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
+        typedef std::set< unsigned > NeighborSet;
+        std::vector<NeighborSet> sparsityPattern(model.numTotalDof());
+
+        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);
+
+                for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
+                    unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
+                    sparsityPattern[myIdx].insert(neighborIdx);
+                }
+            }
+        }
+
+        // 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();
+    }
+
+    // query the problem for all constraint degrees of freedom. note that this method is
+    // quite involved and is thus relatively slow.
+    void updateConstraintsMap_()
+    {
+        if (!enableConstraints_())
+            // constraints are not explictly enabled, so we don't need to consider them!
+            return;
+
+        constraintsMap_.clear();
+
+        // loop over all elements...
+        ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_());
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+        {
+            unsigned threadId = ThreadManager::threadId();
+            ElementIterator elemIt = threadedElemIt.beginParallel();
+            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
+                // create an element context (the solution-based quantities are not
+                // available here!)
+                const Element& elem = *elemIt;
+                ElementContext& elemCtx = *elementCtx_[threadId];
+                elemCtx.updateStencil(elem);
+
+                // check if the problem wants to constrain any degree of the current
+                // element's freedom. if yes, add the constraint to the map.
+                for (unsigned primaryDofIdx = 0;
+                     primaryDofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0);
+                     ++ primaryDofIdx)
+                {
+                    Constraints constraints;
+                    elemCtx.problem().constraints(constraints,
+                                                  elemCtx,
+                                                  primaryDofIdx,
+                                                  /*timeIdx=*/0);
+                    if (constraints.isActive()) {
+                        unsigned globI = elemCtx.globalSpaceIndex(primaryDofIdx, /*timeIdx=*/0);
+                        constraintsMap_[globI] = constraints;
+                        continue;
+                    }
+                }
+            }
+        }
+    }
+
+    // linearize the whole system
+    void linearize_()
+    {
+        resetSystem_();
+
+        // before the first iteration of each time step, we need to update the
+        // constraints. (i.e., we assume that constraints can be time dependent, but they
+        // can't depend on the solution.)
+        if (model_().newtonMethod().numIterations() == 0)
+            updateConstraintsMap_();
+
+        applyConstraintsToSolution_();
+
+        // to avoid a race condition if two threads handle an exception at the same time,
+        // we use an explicit lock to control access to the exception storage object
+        // amongst thread-local handlers
+        std::mutex exceptionLock;
+
+        // storage to any exception that needs to be bridged out of the
+        // parallel block below. initialized to null to indicate no exception
+        std::exception_ptr exceptionPtr = nullptr;
+
+        // relinearize the elements...
+        ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_());
+#ifdef _OPENMP
+#pragma omp parallel
+#endif
+        {
+            ElementIterator elemIt = threadedElemIt.beginParallel();
+            ElementIterator nextElemIt = elemIt;
+            try {
+                for (; !threadedElemIt.isFinished(elemIt); elemIt = nextElemIt) {
+                    // give the model and the problem a chance to prefetch the data required
+                    // to linearize the next element, but only if we need to consider it
+                    nextElemIt = threadedElemIt.increment();
+                    if (!threadedElemIt.isFinished(nextElemIt)) {
+                        const auto& nextElem = *nextElemIt;
+                        if (linearizeNonLocalElements
+                            || nextElem.partitionType() == Dune::InteriorEntity)
+                        {
+                            model_().prefetch(nextElem);
+                            problem_().prefetch(nextElem);
+                        }
+                    }
+
+                    const Element& elem = *elemIt;
+                    if (!linearizeNonLocalElements && elem.partitionType() != Dune::InteriorEntity)
+                        continue;
+
+                    linearizeElement_(elem);
+                }
+            }
+            // If an exception occurs in the parallel block, it won't escape the
+            // block; terminate() is called instead of a handler outside!  hence, we
+            // tuck any exceptions that occur away in the pointer. If an exception
+            // occurs in more than one thread at the same time, we must pick one of
+            // them to be rethrown as we cannot have two active exceptions at the
+            // same time. This solution essentially picks one at random. This will
+            // only be a problem if two different kinds of exceptions are thrown, for
+            // instance if one thread experiences a (recoverable) numerical issue
+            // while another is out of memory.
+            catch(...) {
+                std::lock_guard<std::mutex> take(exceptionLock);
+                exceptionPtr = std::current_exception();
+                threadedElemIt.setFinished();
+            }
+        }  // parallel block
+
+        // after reduction from the parallel block, exceptionPtr will point to
+        // a valid exception if one occurred in one of the threads; rethrow
+        // it here to let the outer handler take care of it properly
+        if(exceptionPtr) {
+            std::rethrow_exception(exceptionPtr);
+        }
+
+        applyConstraintsToLinearization_();
+    }
+
+    // linearize an element in the interior of the process' grid partition
+    void linearizeElement_(const Element& elem)
+    {
+        unsigned threadId = ThreadManager::threadId();
+
+        ElementContext *elementCtx = elementCtx_[threadId];
+        auto& localLinearizer = model_().localLinearizer(threadId);
+
+        // the actual work of linearization is done by the local linearizer class
+        localLinearizer.linearize(*elementCtx, elem);
+
+        // update the right hand side and the Jacobian matrix
+        if (GET_PROP_VALUE(TypeTag, UseLinearizationLock))
+            globalMatrixMutex_.lock();
+
+        size_t numPrimaryDof = elementCtx->numPrimaryDof(/*timeIdx=*/0);
+        for (unsigned primaryDofIdx = 0; primaryDofIdx < numPrimaryDof; ++ primaryDofIdx) {
+            unsigned globI = elementCtx->globalSpaceIndex(/*spaceIdx=*/primaryDofIdx, /*timeIdx=*/0);
+
+            // update the right hand side
+            residual_[globI] += localLinearizer.residual(primaryDofIdx);
+
+            // update the global Jacobian matrix
+            for (unsigned dofIdx = 0; dofIdx < elementCtx->numDof(/*timeIdx=*/0); ++ dofIdx) {
+                unsigned globJ = elementCtx->globalSpaceIndex(/*spaceIdx=*/dofIdx, /*timeIdx=*/0);
+
+                jacobian_->addToBlock(globJ, globI, localLinearizer.jacobian(dofIdx, primaryDofIdx));
+            }
+        }
+
+        if (GET_PROP_VALUE(TypeTag, UseLinearizationLock))
+            globalMatrixMutex_.unlock();
+    }
+
+    // apply the constraints to the solution. (i.e., the solution of constraint degrees
+    // of freedom is set to the value of the constraint.)
+    void applyConstraintsToSolution_()
+    {
+        if (!enableConstraints_())
+            return;
+
+        // TODO: assuming a history size of 2 only works for Euler time discretizations!
+        auto& sol = model_().solution(/*timeIdx=*/0);
+        auto& oldSol = model_().solution(/*timeIdx=*/1);
+
+        auto it = constraintsMap_.begin();
+        const auto& endIt = constraintsMap_.end();
+        for (; it != endIt; ++it) {
+            sol[it->first] = it->second;
+            oldSol[it->first] = it->second;
+        }
+    }
+
+    // apply the constraints to the linearization. (i.e., for constrain degrees of
+    // freedom the Jacobian matrix maps to identity and the residual is zero)
+    void applyConstraintsToLinearization_()
+    {
+        if (!enableConstraints_())
+            return;
+
+        auto it = constraintsMap_.begin();
+        const auto& endIt = constraintsMap_.end();
+        for (; it != endIt; ++it) {
+            unsigned constraintDofIdx = it->first;
+
+            // reset the column of the Jacobian matrix
+            // put an identity matrix on the main diagonal of the Jacobian
+            jacobian_->clearRow(constraintDofIdx, Scalar(1.0));
+
+            // make the right-hand side of constraint DOFs zero
+            residual_[constraintDofIdx] = 0.0;
+        }
+    }
+
+    static bool enableConstraints_()
+    { return GET_PROP_VALUE(TypeTag, EnableConstraints); }
+
+    Simulator *simulatorPtr_;
+    std::vector<ElementContext*> elementCtx_;
+
+    // The constraint equations (only non-empty if the
+    // EnableConstraints property is true)
+    std::map<unsigned, Constraints> constraintsMap_;
+
+    // the jacobian matrix
+    std::unique_ptr<SparseMatrixAdapter> jacobian_;
+
+    // the right-hand side
+    GlobalEqVector residual_;
+
+
+    std::mutex globalMatrixMutex_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaselocalresidual.hh

@@ -0,0 +1,638 @@
+// -*- 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_RESIDUAL_HH
+#define EWOMS_FV_BASE_LOCAL_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 <opm/material/common/Exceptions.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 FvBaseLocalResidual
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryContext) BoundaryContext;
+
+    static constexpr bool useVolumetricResidual = GET_PROP_VALUE(TypeTag, UseVolumetricResidual);
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { extensiveStorageTerm = GET_PROP_VALUE(TypeTag, ExtensiveStorageTerm) };
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+    typedef Dune::FieldVector<Evaluation, numEq> EvalVector;
+
+    // copying the local residual class is not a good idea
+    FvBaseLocalResidual(const FvBaseLocalResidual& )
+    {}
+
+public:
+    typedef Dune::BlockVector<EvalVector, Opm::aligned_allocator<EvalVector, alignof(EvalVector)> > LocalEvalBlockVector;
+
+    FvBaseLocalResidual()
+    { }
+
+    ~FvBaseLocalResidual()
+    { }
+
+    /*!
+     * \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));
+                    Opm::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) {
+                Opm::Valgrind::CheckDefined(storage[dofIdx][eqIdx]);
+                assert(Opm::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();
+
+            Opm::Valgrind::SetUndefined(flux);
+            asImp_().computeFlux(flux, /*context=*/elemCtx, scvfIdx, timeIdx);
+            Opm::Valgrind::CheckDefined(flux);
+#ifndef NDEBUG
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                assert(Opm::isfinite(flux[eqIdx]));
+#endif
+
+            Scalar alpha = elemCtx.extensiveQuantities(scvfIdx, timeIdx).extrusionFactor();
+            alpha *= face.area();
+            Opm::Valgrind::CheckDefined(alpha);
+            assert(alpha > 0.0);
+            assert(Opm::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(Opm::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(Opm::isfinite(residual[i][j]));
+                Opm::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& storage OPM_UNUSED,
+                        const ElementContext& elemCtx OPM_UNUSED,
+                        unsigned dofIdx OPM_UNUSED,
+                        unsigned timeIdx OPM_UNUSED) 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& flux OPM_UNUSED,
+                     const ElementContext& elemCtx OPM_UNUSED,
+                     unsigned scvfIdx OPM_UNUSED,
+                     unsigned timeIdx OPM_UNUSED) 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& source OPM_UNUSED,
+                       const ElementContext& elemCtx OPM_UNUSED,
+                       unsigned dofIdx OPM_UNUSED,
+                       unsigned timeIdx OPM_UNUSED) 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;
+
+        BoundaryContext boundaryCtx(elemCtx);
+
+        // 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(Opm::isfinite(residual[i][j]));
+                Opm::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
+    {
+        BoundaryRateVector values;
+
+        Opm::Valgrind::SetUndefined(values);
+        boundaryCtx.problem().boundary(values, boundaryCtx, boundaryFaceIdx, timeIdx);
+        Opm::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();
+
+            Opm::Valgrind::CheckDefined(values[eqIdx]);
+            assert(Opm::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();
+            Opm::Valgrind::CheckDefined(extrusionFactor);
+            assert(Opm::isfinite(extrusionFactor));
+            assert(extrusionFactor > 0.0);
+            Scalar scvVolume =
+               elemCtx.stencil(/*timeIdx=*/0).subControlVolume(dofIdx).volume() * extrusionFactor;
+            Opm::Valgrind::CheckDefined(scvVolume);
+            assert(Opm::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
+            Opm::Valgrind::CheckDefined(tmp);
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
+                assert(Opm::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]);
+                    }
+
+                    Opm::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);
+                    Opm::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);
+                Opm::Valgrind::CheckDefined(tmp2);
+            }
+
+            // Use the implicit Euler time discretization
+            for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                tmp[eqIdx] -= tmp2[eqIdx];
+                tmp[eqIdx] *= scvVolume / elemCtx.simulator().timeStepSize();
+
+                residual[dofIdx][eqIdx] += tmp[eqIdx];
+            }
+
+            Opm::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] -= Opm::scalarValue(sourceRate[eqIdx])*scvVolume;
+            }
+            else {
+                for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
+                    sourceRate[eqIdx] *= scvVolume;
+                    residual[dofIdx][eqIdx] -= sourceRate[eqIdx];
+                }
+            }
+
+            Opm::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(Opm::isfinite(residual[i][j]));
+                Opm::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

opm/models/discretization/common/fvbasenewtonconvergencewriter.hh

@@ -0,0 +1,156 @@
+// -*- 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::FvBaseNewtonConvergenceWriter
+ */
+#ifndef EWOMS_FV_BASE_NEWTON_CONVERGENCE_WRITER_HH
+#define EWOMS_FV_BASE_NEWTON_CONVERGENCE_WRITER_HH
+
+#include <ewoms/io/vtkmultiwriter.hh>
+#include <opm/models/utils/propertysystem.hh>
+
+#include <iostream>
+
+//! \cond SKIP_THIS
+BEGIN_PROPERTIES
+
+// forward declaration of the required property tags
+NEW_PROP_TAG(GridView);
+NEW_PROP_TAG(NewtonMethod);
+NEW_PROP_TAG(SolutionVector);
+NEW_PROP_TAG(GlobalEqVector);
+NEW_PROP_TAG(VtkOutputFormat);
+
+END_PROPERTIES
+//! \endcond
+
+namespace Opm {
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Writes the intermediate solutions during the Newton scheme
+ *        for models using a finite volume discretization
+ */
+template <class TypeTag>
+class FvBaseNewtonConvergenceWriter
+{
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
+
+    static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat);
+    typedef Opm::VtkMultiWriter<GridView, vtkFormat> VtkMultiWriter;
+
+public:
+    FvBaseNewtonConvergenceWriter(NewtonMethod& nm)
+        : newtonMethod_(nm)
+    {
+        timeStepIdx_ = 0;
+        iteration_ = 0;
+        vtkMultiWriter_ = 0;
+    }
+
+    ~FvBaseNewtonConvergenceWriter()
+    { delete vtkMultiWriter_; }
+
+    /*!
+     * \brief Called by the Newton method before the actual algorithm
+     *        is started for any given timestep.
+     */
+    void beginTimeStep()
+    {
+        ++timeStepIdx_;
+        iteration_ = 0;
+    }
+
+    /*!
+     * \brief Called by the Newton method before an iteration of the
+     *        Newton algorithm is started.
+     */
+    void beginIteration()
+    {
+        ++ iteration_;
+        if (!vtkMultiWriter_)
+            vtkMultiWriter_ =
+                new VtkMultiWriter(/*async=*/false,
+                                   newtonMethod_.problem().gridView(),
+                                   newtonMethod_.problem().outputDir(),
+                                   "convergence");
+        vtkMultiWriter_->beginWrite(timeStepIdx_ + iteration_ / 100.0);
+    }
+
+    /*!
+     * \brief Write the Newton update to disk.
+     *
+     * Called after the linear solution is found for an iteration.
+     *
+     * \param uLastIter The solution vector of the previous iteration.
+     * \param deltaU The negative difference between the solution
+     *        vectors of the previous and the current iteration.
+     */
+    void writeFields(const SolutionVector& uLastIter,
+                     const GlobalEqVector& deltaU)
+    {
+        try {
+            newtonMethod_.problem().model().addConvergenceVtkFields(*vtkMultiWriter_,
+                                                                    uLastIter,
+                                                                    deltaU);
+        }
+        catch (...) {
+            std::cout << "Oops: exception thrown on rank "
+                      << newtonMethod_.problem().gridView().comm().rank()
+                      << " while writing the convergence\n"  << std::flush;
+        }
+    }
+
+    /*!
+     * \brief Called by the Newton method after an iteration of the
+     *        Newton algorithm has been completed.
+     */
+    void endIteration()
+    { vtkMultiWriter_->endWrite(); }
+
+    /*!
+     * \brief Called by the Newton method after Newton algorithm
+     *        has been completed for any given timestep.
+     *
+     * This method is called regardless of whether the Newton method
+     * converged or not.
+     */
+    void endTimeStep()
+    { iteration_ = 0; }
+
+private:
+    int timeStepIdx_;
+    int iteration_;
+    VtkMultiWriter *vtkMultiWriter_;
+    NewtonMethod& newtonMethod_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbasenewtonmethod.hh

@@ -0,0 +1,179 @@
+// -*- 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::FvBaseNewtonMethod
+ */
+#ifndef EWOMS_FV_BASE_NEWTON_METHOD_HH
+#define EWOMS_FV_BASE_NEWTON_METHOD_HH
+
+#include "fvbasenewtonconvergencewriter.hh"
+
+#include <ewoms/nonlinear/newtonmethod.hh>
+#include <opm/models/utils/propertysystem.hh>
+
+namespace Opm {
+
+template <class TypeTag>
+class FvBaseNewtonMethod;
+
+template <class TypeTag>
+class FvBaseNewtonConvergenceWriter;
+} // namespace Opm
+
+BEGIN_PROPERTIES
+
+//! create a type tag for the Newton method of the finite-volume discretization
+NEW_TYPE_TAG(FvBaseNewtonMethod, INHERITS_FROM(NewtonMethod));
+
+//! The class dealing with the balance equations
+NEW_PROP_TAG(Model);
+
+//! The class storing primary variables plus pseudo primary variables
+NEW_PROP_TAG(PrimaryVariables);
+
+//! The class storing values of conservation equations (e.g., a "naked" primary varible
+//! vector)
+NEW_PROP_TAG(EqVector);
+
+//! The number of balance equations.
+NEW_PROP_TAG(NumEq);
+
+//! The discretization specific part of he implementing the Newton algorithm
+NEW_PROP_TAG(DiscNewtonMethod);
+
+//! The class implementing the Newton algorithm
+NEW_PROP_TAG(NewtonMethod);
+
+// set default values
+SET_TYPE_PROP(FvBaseNewtonMethod, DiscNewtonMethod,
+              Opm::FvBaseNewtonMethod<TypeTag>);
+SET_TYPE_PROP(FvBaseNewtonMethod, NewtonMethod,
+              typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod));
+SET_TYPE_PROP(FvBaseNewtonMethod, NewtonConvergenceWriter,
+              Opm::FvBaseNewtonConvergenceWriter<TypeTag>);
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief A Newton method for models using a finite volume discretization.
+ *
+ * This class is sufficient for most models which use an Element or a
+ * Vertex Centered Finite Volume discretization.
+ */
+template <class TypeTag>
+class FvBaseNewtonMethod : public NewtonMethod<TypeTag>
+{
+    typedef Opm::NewtonMethod<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
+    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
+    typedef typename GET_PROP_TYPE(TypeTag, GlobalEqVector) GlobalEqVector;
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
+
+
+public:
+    FvBaseNewtonMethod(Simulator& simulator)
+        : ParentType(simulator)
+    { }
+
+protected:
+    friend class Opm::NewtonMethod<TypeTag>;
+
+    /*!
+     * \brief Update the current solution with a delta vector.
+     *
+     * The error estimates required for the converged() and
+     * proceed() methods should be updated inside this method.
+     *
+     * Different update strategies, such as line search and chopped
+     * updates can be implemented. The default behavior is just to
+     * subtract deltaU from uLastIter, i.e.
+     * \f[ u^{k+1} = u^k - \Delta u^k \f]
+     *
+     * \param nextSolution The solution vector at the end of the current iteration
+     * \param currentSolution The solution vector at the beginning of the current iteration
+     * \param solutionUpdate The delta as calculated by solving the linear system of
+     *                       equations. This parameter also stores the updated solution.
+     * \param currentResidual The residual (i.e., right-hand-side) of the current solution.
+     */
+    void update_(SolutionVector& nextSolution,
+                 const SolutionVector& currentSolution,
+                 const GlobalEqVector& solutionUpdate,
+                 const GlobalEqVector& currentResidual)
+    {
+        ParentType::update_(nextSolution, currentSolution, solutionUpdate, currentResidual);
+
+        // make sure that the intensive quantities get recalculated at the next
+        // linearization
+        if (model_().storeIntensiveQuantities()) {
+            for (unsigned dofIdx = 0; dofIdx < model_().numGridDof(); ++dofIdx)
+                model_().setIntensiveQuantitiesCacheEntryValidity(dofIdx,
+                                                                  /*timeIdx=*/0,
+                                                                  /*valid=*/false);
+        }
+    }
+
+    /*!
+     * \brief Indicates the beginning of a Newton iteration.
+     */
+    void beginIteration_()
+    {
+        model_().syncOverlap();
+
+        ParentType::beginIteration_();
+    }
+
+    /*!
+     * \brief Returns a reference to the model.
+     */
+    Model& model_()
+    { return ParentType::model(); }
+
+    /*!
+     * \brief Returns a reference to the model.
+     */
+    const Model& model_() const
+    { return ParentType::model(); }
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseprimaryvariables.hh

@@ -0,0 +1,166 @@
+// -*- 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::FvBasePrimaryVariables
+ */
+#ifndef EWOMS_FV_BASE_PRIMARY_VARIABLES_HH
+#define EWOMS_FV_BASE_PRIMARY_VARIABLES_HH
+
+#include <type_traits>
+
+#include "fvbaseproperties.hh"
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <dune/common/fvector.hh>
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Represents the primary variables used by the a model.
+ */
+template <class TypeTag>
+class FvBasePrimaryVariables
+    : public Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
+                               GET_PROP_VALUE(TypeTag, NumEq)>
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+
+    typedef Opm::MathToolbox<Evaluation> Toolbox;
+    typedef Dune::FieldVector<Scalar, numEq> ParentType;
+
+public:
+    FvBasePrimaryVariables()
+        : ParentType()
+    { Opm::Valgrind::SetUndefined(*this); }
+
+    /*!
+     * \brief Construction from a scalar value
+     */
+    FvBasePrimaryVariables(Scalar value)
+        : ParentType(value)
+    { }
+
+    /*!
+     * \brief Assignment from another primary variables object
+     */
+    FvBasePrimaryVariables(const FvBasePrimaryVariables& value) = default;
+
+    /*!
+     * \brief Assignment from another primary variables object
+     */
+    FvBasePrimaryVariables& operator=(const FvBasePrimaryVariables& value) = default;
+
+    /*!
+     * \brief Return a primary variable intensive evaluation.
+     *
+     * i.e., the result represents the function f = x_i if the time index is zero, else
+     * it represents the a constant f = x_i. (the difference is that in the first case,
+     * the derivative w.r.t. x_i is 1, while it is 0 in the second case.
+     */
+    Evaluation makeEvaluation(unsigned varIdx, unsigned timeIdx) const
+    {
+        if (std::is_same<Evaluation, Scalar>::value)
+            return (*this)[varIdx]; // finite differences
+        else {
+            // automatic differentiation
+            if (timeIdx == 0)
+                return Toolbox::createVariable((*this)[varIdx], varIdx);
+            else
+                return Toolbox::createConstant((*this)[varIdx]);
+        }
+    }
+
+    /*!
+     * \brief Assign the primary variables "somehow" from a fluid state
+     *
+     * That is without considering any consistency issues which the
+     * fluid state might have. This method is guaranteed to produce
+     * consistent results if the fluid state is consistent to the
+     * properties at a given spatial location. (Where "consistent
+     * results" means that the same fluid state can be reconstructed
+     * from the primary variables.)
+     */
+    template <class FluidState>
+    void assignNaive(const FluidState& fluidState OPM_UNUSED)
+    {
+        throw std::runtime_error("The PrimaryVariables class does not define "
+                                 "an assignNaive() method");
+    }
+
+    /*!
+     * \brief Instruct valgrind to check the definedness of all attributes of this class.
+     */
+    void checkDefined() const
+    {
+        Opm::Valgrind::CheckDefined(*static_cast<const ParentType*>(this));
+    }
+};
+
+} // namespace Opm
+
+namespace Dune {
+
+  /** Compatibility traits class for DenseVector and DenseMatrix.
+   */
+  template<class TypeTag, bool>
+  struct FieldTraitsImpl;
+
+  /** FieldTraitsImpl for classes derived from
+   * Opm::FvBasePrimaryVariables: use FieldVector's FieldTraits implementation) */
+  template<class TypeTag>
+  struct FieldTraitsImpl< TypeTag, true >
+      : public FieldTraits<FieldVector<typename GET_PROP_TYPE(TypeTag, Scalar),
+                                       GET_PROP_VALUE(TypeTag, NumEq)> >
+  {
+  };
+
+  /** FieldTraitsImpl for classes not derived from
+   * Opm::FvBasePrimaryVariables, fall bakc to existing implementation */
+  template<class T>
+  struct FieldTraitsImpl< T, false >
+    : public FieldTraits< T >
+  {
+  };
+
+
+  /** Specialization of FieldTraits for all PrimaryVariables derived from Opm::FvBasePrimaryVariables */
+  template<class TypeTag, template <class> class EwomsPrimaryVariable>
+  struct FieldTraits< EwomsPrimaryVariable< TypeTag > >
+    : public FieldTraitsImpl< TypeTag,
+                              std::is_base_of< Opm::FvBasePrimaryVariables< TypeTag >,
+                                               EwomsPrimaryVariable< TypeTag > > :: value >
+  {
+  };
+}
+
+#endif

opm/models/discretization/common/fvbaseproblem.hh

@@ -0,0 +1,841 @@
+// -*- 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::FvBaseProblem
+ */
+#ifndef EWOMS_FV_BASE_PROBLEM_HH
+#define EWOMS_FV_BASE_PROBLEM_HH
+
+#include "fvbaseproperties.hh"
+
+#include <ewoms/io/vtkmultiwriter.hh>
+#include <ewoms/io/restart.hh>
+#include <opm/models/discretization/common/restrictprolong.hh>
+
+#include <opm/material/common/Unused.hpp>
+#include <opm/material/common/Exceptions.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <iostream>
+#include <limits>
+#include <string>
+
+#include <sys/stat.h>
+
+namespace Opm {
+
+/*!
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Base class for all problems which use a finite volume spatial discretization.
+ *
+ * \note All quantities are specified assuming a threedimensional world. Problems
+ *       discretized using 2D grids are assumed to be extruded by \f$1 m\f$ and 1D grids
+ *       are assumed to have a cross section of \f$1m \times 1m\f$.
+ */
+template<class TypeTag>
+class FvBaseProblem
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    static const int vtkOutputFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat);
+    typedef Opm::VtkMultiWriter<GridView, vtkOutputFormat> VtkMultiWriter;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager;
+    typedef typename GET_PROP_TYPE(TypeTag, NewtonMethod) NewtonMethod;
+
+    typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
+
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
+
+    enum {
+        dim = GridView::dimension,
+        dimWorld = GridView::dimensionworld
+    };
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<dim>::Entity Vertex;
+    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+
+    typedef typename GridView::Grid::ctype CoordScalar;
+    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+
+public:
+    // the default restriction and prolongation for adaptation is simply an empty one
+    typedef EmptyRestrictProlong  RestrictProlongOperator;
+
+private:
+    // copying a problem is not a good idea
+    FvBaseProblem(const FvBaseProblem& ) = delete;
+
+public:
+    /*!
+     * \copydoc Doxygen::defaultProblemConstructor
+     *
+     * \param simulator The time manager of the simulation
+     * \param gridView The view on the DUNE grid which ought to be
+     *                 used (normally the leaf grid view)
+     */
+    FvBaseProblem(Simulator& simulator)
+        : nextTimeStepSize_(0.0)
+        , gridView_(simulator.gridView())
+#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
+        , elementMapper_(gridView_, Dune::mcmgElementLayout())
+        , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
+#else
+        , elementMapper_(gridView_)
+        , vertexMapper_(gridView_)
+#endif
+        , boundingBoxMin_(std::numeric_limits<double>::max())
+        , boundingBoxMax_(-std::numeric_limits<double>::max())
+        , simulator_(simulator)
+        , defaultVtkWriter_(0)
+    {
+        // calculate the bounding box of the local partition of the grid view
+        VertexIterator vIt = gridView_.template begin<dim>();
+        const VertexIterator vEndIt = gridView_.template end<dim>();
+        for (; vIt!=vEndIt; ++vIt) {
+            for (unsigned i=0; i<dim; i++) {
+                boundingBoxMin_[i] = std::min(boundingBoxMin_[i], vIt->geometry().corner(0)[i]);
+                boundingBoxMax_[i] = std::max(boundingBoxMax_[i], vIt->geometry().corner(0)[i]);
+            }
+        }
+
+        // communicate to get the bounding box of the whole domain
+        for (unsigned i = 0; i < dim; ++i) {
+            boundingBoxMin_[i] = gridView_.comm().min(boundingBoxMin_[i]);
+            boundingBoxMax_[i] = gridView_.comm().max(boundingBoxMax_[i]);
+        }
+
+        if (enableVtkOutput_()) {
+            bool asyncVtkOutput =
+                simulator_.gridView().comm().size() == 1 &&
+                EWOMS_GET_PARAM(TypeTag, bool, EnableAsyncVtkOutput);
+
+            // asynchonous VTK output currently does not work in conjunction with grid
+            // adaptivity because the async-IO code assumes that the grid stays
+            // constant. complain about that case.
+            bool enableGridAdaptation = EWOMS_GET_PARAM(TypeTag, bool, EnableGridAdaptation);
+            if (asyncVtkOutput && enableGridAdaptation)
+                throw std::runtime_error("Asynchronous VTK output currently cannot be used "
+                                         "at the same time as grid adaptivity");
+
+            std::string outputDir = asImp_().outputDir();
+
+            defaultVtkWriter_ =
+                new VtkMultiWriter(asyncVtkOutput, gridView_, outputDir, asImp_().name());
+        }
+    }
+
+    ~FvBaseProblem()
+    { delete defaultVtkWriter_; }
+
+    /*!
+     * \brief Registers all available parameters for the problem and
+     *        the model.
+     */
+    static void registerParameters()
+    {
+        Model::registerParameters();
+        EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxTimeStepSize,
+                             "The maximum size to which all time steps are limited to [s]");
+        EWOMS_REGISTER_PARAM(TypeTag, Scalar, MinTimeStepSize,
+                             "The minimum size to which all time steps are limited to [s]");
+        EWOMS_REGISTER_PARAM(TypeTag, unsigned, MaxTimeStepDivisions,
+                             "The maximum number of divisions by two of the timestep size "
+                             "before the simulation bails out");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, EnableAsyncVtkOutput,
+                             "Dispatch a separate thread to write the VTK output");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, ContinueOnConvergenceError,
+                             "Continue with a non-converged solution instead of giving up "
+                             "if we encounter a time step size smaller than the minimum time "
+                             "step size.");
+    }
+
+    /*!
+     * \brief Return if the storage term of the first iteration is identical to the storage
+     *        term for the solution of the previous time step.
+     *
+     * This is only relevant if the storage cache is enabled and is usually the case,
+     * i.e., this method only needs to be overwritten in rare corner cases.
+     */
+    bool recycleFirstIterationStorage() const
+    { return true; }
+
+    /*!
+     * \brief Determine the directory for simulation output.
+     *
+     * The actual problem may chose to transform the value of the OutputDir parameter and
+     * it can e.g. choose to create the directory on demand if it does not exist. The
+     * default behaviour is to just return the OutputDir parameter and to throw an
+     * exception if no directory with this name exists.
+     */
+    std::string outputDir() const
+    {
+        std::string outputDir = EWOMS_GET_PARAM(TypeTag, std::string, OutputDir);
+
+        if (outputDir == "")
+            throw std::runtime_error("No directory for output specified");
+
+        // TODO: replace this by std::filesystem once we require c++-2017
+        struct stat st;
+        if (::stat(outputDir.c_str(), &st) != 0)
+            throw std::runtime_error("Could not access output directory '"+outputDir+"':"
+                                     +strerror(errno));
+        if (!S_ISDIR(st.st_mode))
+            throw std::runtime_error("Path to output directory '"+outputDir+"' exists but is not a directory");
+
+        if (access(outputDir.c_str(), W_OK) != 0)
+            throw std::runtime_error("Output directory '"+outputDir+"' exists but is not writeable");
+
+        return outputDir;
+    }
+
+    /*!
+     * \brief Returns the string that is printed before the list of command line
+     *        parameters in the help message.
+     *
+     * If the returned string is empty, no help message will be generated.
+     */
+    static std::string helpPreamble(int argc OPM_UNUSED,
+                                    const char **argv)
+    {
+        std::string desc = Implementation::briefDescription();
+        if (!desc.empty())
+            desc = desc + "\n";
+
+        return
+            "Usage: "+std::string(argv[0]) + " [OPTIONS]\n"
+            + desc;
+    }
+
+    /*!
+     * \brief Returns a human readable description of the problem for the help message
+     *
+     * The problem description is printed as part of the --help message. It is optional
+     * and should not exceed one or two lines of text.
+     */
+    static std::string briefDescription()
+    { return ""; }
+
+    // TODO (?): detailedDescription()
+
+    /*!
+     * \brief Handles positional command line parameters.
+     *
+     * Positional parameters are parameters that are not prefixed by any parameter name.
+     *
+     * \param seenParams The parameters which have already been seen in the current context
+     * \param errorMsg If the positional argument cannot be handled, this is the reason why
+     * \param argc The total number of command line parameters
+     * \param argv The string value of the command line parameters
+     * \param paramIdx The index of the positional parameter in the array of all parameters
+     * \param posParamIdx The number of the positional parameter encountered so far
+     *
+     * \return The number of array entries which ought to be skipped before processing
+     *         the next regular parameter. If this is less than 1, it indicated that the
+     *         positional parameter was invalid.
+     */
+    static int handlePositionalParameter(std::set<std::string>& seenParams OPM_UNUSED,
+                                         std::string& errorMsg,
+                                         int argc OPM_UNUSED,
+                                         const char** argv,
+                                         int paramIdx,
+                                         int posParamIdx OPM_UNUSED)
+    {
+        errorMsg = std::string("Illegal parameter \"")+argv[paramIdx]+"\".";
+        return 0;
+    }
+
+    /*!
+     * \brief Called by the Opm::Simulator in order to initialize the problem.
+     *
+     * If you overload this method don't forget to call ParentType::finishInit()
+     */
+    void finishInit()
+    { }
+
+    /*!
+     * \brief Allows to improve the performance by prefetching all data which is
+     *        associated with a given element.
+     */
+    void prefetch(const Element& elem OPM_UNUSED) const
+    {
+        // do nothing by default
+    }
+
+    /*!
+     * \brief Handle changes of the grid
+     */
+    void gridChanged()
+    {
+        elementMapper_.update();
+        vertexMapper_.update();
+
+        if (enableVtkOutput_())
+            defaultVtkWriter_->gridChanged();
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a boundary segment.
+     *
+     * \param values Stores the fluxes over the boundary segment.
+     * \param context The object representing the execution context from
+     *                which this method is called.
+     * \param spaceIdx The local index of the spatial entity which represents the boundary segment.
+     * \param timeIdx The index used for the time discretization
+     */
+    template <class Context>
+    void boundary(BoundaryRateVector& values OPM_UNUSED,
+                  const Context& context OPM_UNUSED,
+                  unsigned spaceIdx OPM_UNUSED,
+                  unsigned timeIdx OPM_UNUSED) const
+    { throw std::logic_error("Problem does not provide a boundary() method"); }
+
+    /*!
+     * \brief Evaluate the constraints for a control volume.
+     *
+     * \param constraints Stores the values of the primary variables at a
+     *                    given spatial and temporal location.
+     * \param context The object representing the execution context from
+     *                which this method is called.
+     * \param spaceIdx The local index of the spatial entity which represents the boundary segment.
+     * \param timeIdx The index used for the time discretization
+     */
+    template <class Context>
+    void constraints(Constraints& constrs OPM_UNUSED,
+                     const Context& context OPM_UNUSED,
+                     unsigned spaceIdx OPM_UNUSED,
+                     unsigned timeIdx OPM_UNUSED) const
+    { throw std::logic_error("Problem does not provide a constraints() method"); }
+
+    /*!
+     * \brief Evaluate the source term for all phases within a given
+     *        sub-control-volume.
+     *
+     * \param rate Stores the values of the volumetric creation/anihilition
+     *             rates of the conserved quantities.
+     * \param context The object representing the execution context from which
+     *                this method is called.
+     * \param spaceIdx The local index of the spatial entity which represents
+     *                 the boundary segment.
+     * \param timeIdx The index used for the time discretization
+     */
+    template <class Context>
+    void source(RateVector& rate OPM_UNUSED,
+                const Context& context OPM_UNUSED,
+                unsigned spaceIdx OPM_UNUSED,
+                unsigned timeIdx OPM_UNUSED) const
+    { throw std::logic_error("Problem does not provide a source() method"); }
+
+    /*!
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * \param values Stores the primary variables.
+     * \param context The object representing the execution context from which
+     *                this method is called.
+     * \param spaceIdx The local index of the spatial entity which represents
+     *                 the boundary segment.
+     * \param timeIdx The index used for the time discretization
+     */
+    template <class Context>
+    void initial(PrimaryVariables& values OPM_UNUSED,
+                 const Context& context OPM_UNUSED,
+                 unsigned spaceIdx OPM_UNUSED,
+                 unsigned timeIdx OPM_UNUSED) const
+    { throw std::logic_error("Problem does not provide a initial() method"); }
+
+    /*!
+     * \brief Return how much the domain is extruded at a given sub-control volume.
+     *
+     * This means the factor by which a lower-dimensional (1D or 2D)
+     * entity needs to be expanded to get a full dimensional cell. The
+     * default is 1.0 which means that 1D problems are actually
+     * thought as pipes with a cross section of 1 m^2 and 2D problems
+     * are assumed to extend 1 m to the back.
+     *
+     * \param context The object representing the execution context from which
+     *                this method is called.
+     * \param spaceIdx The local index of the spatial entity which represents
+     *                 the boundary segment.
+     * \param timeIdx The index used for the time discretization
+     */
+    template <class Context>
+    Scalar extrusionFactor(const Context& context OPM_UNUSED,
+                           unsigned spaceIdx OPM_UNUSED,
+                           unsigned timeIdx OPM_UNUSED) const
+    { return asImp_().extrusionFactor(); }
+
+    Scalar extrusionFactor() const
+    { return 1.0; }
+
+    /*!
+     * \brief Callback used by the model to indicate that the initial solution has been
+     *        determined for all degrees of freedom.
+     */
+    void initialSolutionApplied()
+    {}
+
+    /*!
+     * \brief Called at the beginning of an simulation episode.
+     */
+    void beginEpisode()
+    { }
+
+    /*!
+     * \brief Called by the simulator before each time integration.
+     */
+    void beginTimeStep()
+    { }
+
+    /*!
+     * \brief Called by the simulator before each Newton-Raphson iteration.
+     */
+    void beginIteration()
+    { }
+
+    /*!
+     * \brief Called by the simulator after each Newton-Raphson update.
+     */
+    void endIteration()
+    { }
+
+    /*!
+     * \brief Called by the simulator after each time integration.
+     *
+     * This method is intended to do some post processing of the
+     * solution. (e.g., some additional output)
+     */
+    void endTimeStep()
+    { }
+
+    /*!
+     * \brief Called when the end of an simulation episode is reached.
+     *
+     * Typically, a new episode is started in this method.
+     */
+    void endEpisode()
+    {
+        std::cerr << "The end of episode " << simulator().episodeIndex() + 1 << " has been "
+                  << "reached, but the problem does not override the endEpisode() method. "
+                  << "Doing nothing!\n";
+    }
+
+    /*!
+     * \brief Called after the simulation has been run sucessfully.
+     */
+    void finalize()
+    {
+        const auto& executionTimer = simulator().executionTimer();
+
+        Scalar executionTime = executionTimer.realTimeElapsed();
+        Scalar setupTime = simulator().setupTimer().realTimeElapsed();
+        Scalar prePostProcessTime = simulator().prePostProcessTimer().realTimeElapsed();
+        Scalar localCpuTime = executionTimer.cpuTimeElapsed();
+        Scalar globalCpuTime = executionTimer.globalCpuTimeElapsed();
+        Scalar writeTime = simulator().writeTimer().realTimeElapsed();
+        Scalar linearizeTime = simulator().linearizeTimer().realTimeElapsed();
+        Scalar solveTime = simulator().solveTimer().realTimeElapsed();
+        Scalar updateTime = simulator().updateTimer().realTimeElapsed();
+        unsigned numProcesses = static_cast<unsigned>(this->gridView().comm().size());
+        unsigned threadsPerProcess = ThreadManager::maxThreads();
+        if (gridView().comm().rank() == 0) {
+            std::cout << std::setprecision(3)
+                      << "Simulation of problem '" << asImp_().name() << "' finished.\n"
+                      << "\n"
+                      << "------------------------ Timing receipt ------------------------\n"
+                      << "Setup time: " << setupTime << " seconds" << Simulator::humanReadableTime(setupTime)
+                      << ", " << setupTime/(executionTime + setupTime)*100 << "%\n"
+                      << "Simulation time: " << executionTime << " seconds" << Simulator::humanReadableTime(executionTime)
+                      << ", " << executionTime/(executionTime + setupTime)*100 << "%\n"
+                      << "    Linearization time: " << linearizeTime << " seconds" << Simulator::humanReadableTime(linearizeTime)
+                      << ", " << linearizeTime/executionTime*100 << "%\n"
+                      << "    Linear solve time: "  << solveTime << " seconds" << Simulator::humanReadableTime(solveTime)
+                      << ", " << solveTime/executionTime*100 << "%\n"
+                      << "    Newton update time: "  << updateTime << " seconds" << Simulator::humanReadableTime(updateTime)
+                      << ", " << updateTime/executionTime*100 << "%\n"
+                      << "    Pre/postprocess time: "  << prePostProcessTime << " seconds" << Simulator::humanReadableTime(prePostProcessTime)
+                      << ", " << prePostProcessTime/executionTime*100 << "%\n"
+                      << "    Output write time: "  << writeTime << " seconds" << Simulator::humanReadableTime(writeTime)
+                      << ", " << writeTime/executionTime*100 << "%\n"
+                      << "First process' simulation CPU time: "  << localCpuTime << " seconds" <<  Simulator::humanReadableTime(localCpuTime) << "\n"
+                      << "Number of processes: " << numProcesses << "\n"
+                      << "Threads per processes: " << threadsPerProcess << "\n"
+                      << "Total CPU time: " << globalCpuTime << " seconds" << Simulator::humanReadableTime(globalCpuTime) << "\n"
+                      << "\n"
+                      << "Note 1: If not stated otherwise, all times are wall clock times\n"
+                      << "Note 2: Taxes and administrative overhead are "
+                      << (executionTime - (linearizeTime+solveTime+updateTime+prePostProcessTime+writeTime))/executionTime*100
+                      << "%\n"
+                      << "\n"
+                      << "Our simulation hours are 24/7. Thank you for choosing us.\n"
+                      << "----------------------------------------------------------------\n"
+                      << "\n"
+                      << std::flush;
+        }
+    }
+
+    /*!
+     * \brief Called by Opm::Simulator in order to do a time
+     *        integration on the model.
+     */
+    void timeIntegration()
+    {
+        unsigned maxFails = asImp_().maxTimeIntegrationFailures();
+        Scalar minTimeStepSize = asImp_().minTimeStepSize();
+
+        std::string errorMessage;
+        for (unsigned i = 0; i < maxFails; ++i) {
+            bool converged = model().update();
+            if (converged)
+                return;
+
+            Scalar dt = simulator().timeStepSize();
+            Scalar nextDt = dt / 2.0;
+            if (dt < minTimeStepSize*(1 + 1e-9)) {
+                if (asImp_().continueOnConvergenceError()) {
+                    if (gridView().comm().rank() == 0)
+                        std::cout << "Newton solver did not converge with minimum time step of "
+                                  << dt << " seconds. Continuing with unconverged solution!\n"
+                                  << std::flush;
+                    return;
+                }
+                else {
+                    errorMessage =
+                        "Time integration did not succeed with the minumum time step size of "
+                        + std::to_string(double(minTimeStepSize)) + " seconds. Giving up!";
+                    break; // give up: we can't make the time step smaller anymore!
+                }
+            }
+            else if (nextDt < minTimeStepSize)
+                nextDt = minTimeStepSize;
+            simulator().setTimeStepSize(nextDt);
+
+            // update failed
+            if (gridView().comm().rank() == 0)
+                std::cout << "Newton solver did not converge with "
+                          << "dt=" << dt << " seconds. Retrying with time step of "
+                          << nextDt << " seconds\n" << std::flush;
+        }
+
+        if (errorMessage.empty())
+            errorMessage =
+                "Newton solver didn't converge after "
+                +std::to_string(maxFails)+" time-step divisions. dt="
+                +std::to_string(double(simulator().timeStepSize()));
+        throw std::runtime_error(errorMessage);
+    }
+
+    /*!
+     * \brief Returns the minimum allowable size of a time step.
+     */
+    Scalar minTimeStepSize() const
+    { return EWOMS_GET_PARAM(TypeTag, Scalar, MinTimeStepSize); }
+
+    /*!
+     * \brief Returns the maximum number of subsequent failures for the time integration
+     *        before giving up.
+     */
+    unsigned maxTimeIntegrationFailures() const
+    { return EWOMS_GET_PARAM(TypeTag, unsigned, MaxTimeStepDivisions); }
+
+    /*!
+     * \brief Returns if we should continue with a non-converged solution instead of
+     *        giving up if we encounter a time step size smaller than the minimum time
+     *        step size.
+     */
+    bool continueOnConvergenceError() const
+    { return EWOMS_GET_PARAM(TypeTag, unsigned, ContinueOnConvergenceError); }
+
+    /*!
+     * \brief Impose the next time step size to be used externally.
+     */
+    void setNextTimeStepSize(Scalar dt)
+    { nextTimeStepSize_ = dt; }
+
+    /*!
+     * \brief Called by Opm::Simulator whenever a solution for a
+     *        time step has been computed and the simulation time has
+     *        been updated.
+     */
+    Scalar nextTimeStepSize() const
+    {
+        if (nextTimeStepSize_ > 0.0)
+            return nextTimeStepSize_;
+
+        Scalar dtNext = std::min(EWOMS_GET_PARAM(TypeTag, Scalar, MaxTimeStepSize),
+                                 newtonMethod().suggestTimeStepSize(simulator().timeStepSize()));
+
+        if (dtNext < simulator().maxTimeStepSize()
+            && simulator().maxTimeStepSize() < dtNext*2)
+        {
+            dtNext = simulator().maxTimeStepSize()/2 * 1.01;
+        }
+
+        return dtNext;
+    }
+
+    /*!
+     * \brief Returns true if a restart file should be written to
+     *        disk.
+     *
+     * The default behavior is to write one restart file every 10 time
+     * steps. This method should be overwritten by the
+     * implementation if the default behavior is deemed insufficient.
+     */
+    bool shouldWriteRestartFile() const
+    {
+        return simulator().timeStepIndex() > 0 &&
+            (simulator().timeStepIndex() % 10 == 0);
+    }
+
+    /*!
+     * \brief Returns true if the current solution should be written to
+     *        disk (i.e. as a VTK file)
+     *
+     * The default behavior is to write out the solution for every
+     * time step. This method is should be overwritten by the
+     * implementation if the default behavior is deemed insufficient.
+     */
+    bool shouldWriteOutput() const
+    { return true; }
+
+    /*!
+     * \brief Called by the simulator after everything which can be
+     *        done about the current time step is finished and the
+     *        model should be prepared to do the next time integration.
+     */
+    void advanceTimeLevel()
+    { model().advanceTimeLevel(); }
+
+    /*!
+     * \brief The problem name.
+     *
+     * This is used as a prefix for files generated by the simulation.
+     * It is highly recommend to overwrite this method in the concrete
+     * problem which is simulated.
+     */
+    std::string name() const
+    { return "sim"; }
+
+    /*!
+     * \brief The GridView which used by the problem.
+     */
+    const GridView& gridView() const
+    { return gridView_; }
+
+    /*!
+     * \brief The coordinate of the corner of the GridView's bounding
+     *        box with the smallest values.
+     */
+    const GlobalPosition& boundingBoxMin() const
+    { return boundingBoxMin_; }
+
+    /*!
+     * \brief The coordinate of the corner of the GridView's bounding
+     *        box with the largest values.
+     */
+    const GlobalPosition& boundingBoxMax() const
+    { return boundingBoxMax_; }
+
+    /*!
+     * \brief Returns the mapper for vertices to indices.
+     */
+    const VertexMapper& vertexMapper() const
+    { return vertexMapper_; }
+
+    /*!
+     * \brief Returns the mapper for elements to indices.
+     */
+    const ElementMapper& elementMapper() const
+    { return elementMapper_; }
+
+    /*!
+     * \brief Returns Simulator object used by the simulation
+     */
+    Simulator& simulator()
+    { return simulator_; }
+
+    /*!
+     * \copydoc simulator()
+     */
+    const Simulator& simulator() const
+    { return simulator_; }
+
+    /*!
+     * \brief Returns numerical model used for the problem.
+     */
+    Model& model()
+    { return simulator_.model(); }
+
+    /*!
+     * \copydoc model()
+     */
+    const Model& model() const
+    { return simulator_.model(); }
+
+    /*!
+     * \brief Returns object which implements the Newton method.
+     */
+    NewtonMethod& newtonMethod()
+    { return model().newtonMethod(); }
+
+    /*!
+     * \brief Returns object which implements the Newton method.
+     */
+    const NewtonMethod& newtonMethod() const
+    { return model().newtonMethod(); }
+    // \}
+
+    /*!
+     * \brief return restriction and prolongation operator
+     * \note This method has to be overloaded by the implementation.
+     */
+    RestrictProlongOperator restrictProlongOperator()
+    {
+        return RestrictProlongOperator();
+    }
+
+    /*!
+     * \brief Mark grid cells for refinement or coarsening
+     * \note This method has to be overloaded in derived classes to proper implement
+     *       marking of grid entities.
+     *
+     * \return number of marked cells (default is 0)
+     */
+    unsigned markForGridAdaptation()
+    {
+        return 0;
+    }
+
+    /*!
+     * \brief This method writes the complete state of the problem
+     *        to the harddisk.
+     *
+     * The file will start with the prefix returned by the name()
+     * method, has the current time of the simulation clock in it's
+     * name and uses the extension <tt>.ers</tt>. (Ewoms ReStart
+     * file.)  See Opm::Restart for details.
+     *
+     * \tparam Restarter The serializer type
+     *
+     * \param res The serializer object
+     */
+    template <class Restarter>
+    void serialize(Restarter& res)
+    {
+        if (enableVtkOutput_())
+            defaultVtkWriter_->serialize(res);
+    }
+
+    /*!
+     * \brief This method restores the complete state of the problem
+     *        from disk.
+     *
+     * It is the inverse of the serialize() method.
+     *
+     * \tparam Restarter The deserializer type
+     *
+     * \param res The deserializer object
+     */
+    template <class Restarter>
+    void deserialize(Restarter& res)
+    {
+        if (enableVtkOutput_())
+            defaultVtkWriter_->deserialize(res);
+    }
+
+    /*!
+     * \brief Write the relevant secondary variables of the current
+     *        solution into an VTK output file.
+     *
+     * \param verbose Specify if a message should be printed whenever a file is written
+     */
+    void writeOutput(bool verbose = true)
+    {
+        if (!enableVtkOutput_())
+            return;
+
+        if (verbose && gridView().comm().rank() == 0)
+            std::cout << "Writing visualization results for the current time step.\n"
+                      << std::flush;
+
+        // calculate the time _after_ the time was updated
+        Scalar t = simulator().time() + simulator().timeStepSize();
+
+        defaultVtkWriter_->beginWrite(t);
+        model().prepareOutputFields();
+        model().appendOutputFields(*defaultVtkWriter_);
+        defaultVtkWriter_->endWrite();
+
+    }
+
+    /*!
+     * \brief Method to retrieve the VTK writer which should be used
+     *        to write the default ouput after each time step to disk.
+     */
+    VtkMultiWriter& defaultVtkWriter() const
+    { return defaultVtkWriter_; }
+
+protected:
+    Scalar nextTimeStepSize_;
+
+private:
+    bool enableVtkOutput_() const
+    { return EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput); }
+
+    //! Returns the implementation of the problem (i.e. static polymorphism)
+    Implementation& asImp_()
+    { return *static_cast<Implementation *>(this); }
+
+    //! \copydoc asImp_()
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation *>(this); }
+
+    // Grid management stuff
+    const GridView gridView_;
+    ElementMapper elementMapper_;
+    VertexMapper vertexMapper_;
+    GlobalPosition boundingBoxMin_;
+    GlobalPosition boundingBoxMax_;
+
+    // Attributes required for the actual simulation
+    Simulator& simulator_;
+    mutable VtkMultiWriter *defaultVtkWriter_;
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/common/fvbaseproperties.hh

@@ -0,0 +1,324 @@
+// -*- 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
+ * \ingroup FiniteVolumeDiscretizations
+ *
+ * \brief Declare the properties used by the infrastructure code of
+ *        the finite volume discretizations.
+ */
+#ifndef EWOMS_FV_BASE_PROPERTIES_HH
+#define EWOMS_FV_BASE_PROPERTIES_HH
+
+#include "fvbasenewtonmethod.hh"
+#include "fvbaseproperties.hh"
+#include "fvbasefdlocallinearizer.hh"
+
+#include <opm/models/utils/basicproperties.hh>
+#include <ewoms/io/vtkprimaryvarsmodule.hh>
+#include <opm/simulators/linalg/parallelbicgstabbackend.hh>
+
+BEGIN_PROPERTIES
+
+//! The type tag for models based on the finite volume schemes
+NEW_TYPE_TAG(FvBaseDiscretization,
+             INHERITS_FROM(ImplicitModel,
+                           FvBaseNewtonMethod,
+                           VtkPrimaryVars));
+
+
+//! set the splices for the finite volume discretizations
+NEW_PROP_TAG(LinearSolverSplice);
+NEW_PROP_TAG(ParallelBiCGStabLinearSolver);
+
+NEW_PROP_TAG(LocalLinearizerSplice);
+NEW_PROP_TAG(FiniteDifferenceLocalLinearizer);
+
+SET_SPLICES(FvBaseDiscretization, LinearSolverSplice, LocalLinearizerSplice);
+
+//! use a parallel BiCGStab linear solver by default
+SET_TAG_PROP(FvBaseDiscretization, LinearSolverSplice, ParallelBiCGStabLinearSolver);
+
+//! by default, use finite differences to linearize the system of PDEs
+SET_TAG_PROP(FvBaseDiscretization, LocalLinearizerSplice, FiniteDifferenceLocalLinearizer);
+
+/*!
+ * \brief Representation of a function evaluation and all necessary derivatives with
+ *        regard to the intensive quantities of the primary variables.
+ *
+ * Depending on the chosen linearization method, this property may be the same as the
+ * "Scalar" property (if the finite difference linearizer is used), or it may be more
+ * complex (for the linearizer which uses automatic differentiation).
+ */
+NEW_PROP_TAG(Evaluation);
+
+//! The type of the DUNE grid
+NEW_PROP_TAG(Grid);
+//! The type of the grid view
+NEW_PROP_TAG(GridView);
+
+//! The class describing the stencil of the spatial discretization
+NEW_PROP_TAG(Stencil);
+
+//! The class describing the discrete function space when dune-fem is used, otherwise it points to the stencil class
+NEW_PROP_TAG(DiscreteFunctionSpace);
+
+//! The type of the problem
+NEW_PROP_TAG(Problem);
+//! The type of the base class for all problems which use this model
+NEW_PROP_TAG(BaseProblem);
+//! The type of the model
+NEW_PROP_TAG(Model);
+//! Number of equations in the system of PDEs
+NEW_PROP_TAG(NumEq);
+
+//! The type of the spatial discretization used by the model
+NEW_PROP_TAG(Discretization);
+//! The discretization specific part of the local residual
+NEW_PROP_TAG(DiscLocalResidual);
+//! The type of the local residual function
+NEW_PROP_TAG(LocalResidual);
+//! The type of the local linearizer
+NEW_PROP_TAG(LocalLinearizer);
+//! Specify if elements that do not belong to the local process' grid partition should be
+//! skipped
+NEW_PROP_TAG(LinearizeNonLocalElements);
+
+//! Linearizes the global non-linear system of equations
+NEW_PROP_TAG(BaseLinearizer);
+//! The class that allows to manipulate sparse matrices
+NEW_PROP_TAG(SparseMatrixAdapter);
+
+//! A vector of holding a quantity for each equation (usually at a given spatial location)
+NEW_PROP_TAG(EqVector);
+//! A vector of holding a quantity for each equation for each DOF of an element
+NEW_PROP_TAG(ElementEqVector);
+//! Vector containing a quantity of for equation for each DOF of the whole grid
+NEW_PROP_TAG(GlobalEqVector);
+
+//! Vector containing volumetric or areal rates of quantities
+NEW_PROP_TAG(RateVector);
+//! Type of object for specifying boundary conditions
+NEW_PROP_TAG(BoundaryRateVector);
+//! The class which represents a constraint degree of freedom
+NEW_PROP_TAG(Constraints);
+
+//! Vector containing all primary variables of the grid
+NEW_PROP_TAG(SolutionVector);
+
+//! A vector of primary variables within a sub-control volume
+NEW_PROP_TAG(PrimaryVariables);
+//! The secondary variables within a sub-control volume
+NEW_PROP_TAG(IntensiveQuantities);
+//! The discretization specific part of the intensive quantities
+NEW_PROP_TAG(DiscIntensiveQuantities);
+
+//! The secondary variables of all degrees of freedom in an element's stencil
+NEW_PROP_TAG(ElementContext);
+//! The secondary variables of a boundary segment
+NEW_PROP_TAG(BoundaryContext);
+//! The secondary variables of a constraint degree of freedom
+NEW_PROP_TAG(ConstraintsContext);
+//! Data required to calculate a flux over a face
+NEW_PROP_TAG(ExtensiveQuantities);
+//! Calculates gradients of arbitrary quantities at flux integration points
+NEW_PROP_TAG(GradientCalculator);
+
+//! The part of the intensive quantities which is specific to the spatial discretization
+NEW_PROP_TAG(DiscBaseIntensiveQuantities);
+
+//! The part of the extensive quantities which is specific to the spatial discretization
+NEW_PROP_TAG(DiscExtensiveQuantities);
+
+//! The part of the VTK ouput modules which is specific to the spatial discretization
+NEW_PROP_TAG(DiscBaseOutputModule);
+
+//! The class to create grid communication handles
+NEW_PROP_TAG(GridCommHandleFactory);
+
+/*!
+ * \brief The OpenMP threads manager
+ */
+NEW_PROP_TAG(ThreadManager);
+NEW_PROP_TAG(ThreadsPerProcess);
+
+//! use locking to prevent race conditions when linearizing the global system of
+//! equations in multi-threaded mode. (setting this property to true is always save, but
+//! it may slightly deter performance in multi-threaded simlations and some
+//! discretizations do not need this.)
+NEW_PROP_TAG(UseLinearizationLock);
+
+// high-level simulation control
+
+//! Manages the simulation time
+NEW_PROP_TAG(Simulator);
+
+/*!
+ * \brief Switch to enable or disable grid adaptation
+ *
+ * Currently grid adaptation requires the presence of the dune-FEM module. If it is not
+ * available and grid adaptation is enabled, an exception is thrown.
+ */
+NEW_PROP_TAG(EnableGridAdaptation);
+
+/*!
+ * \brief The directory to which simulation output ought to be written to.
+ */
+NEW_PROP_TAG(OutputDir);
+
+/*!
+ * \brief Global switch to enable or disable the writing of VTK output files
+ *
+ * If writing VTK files is disabled, then the WriteVtk$FOO options do
+ * not have any effect...
+ */
+NEW_PROP_TAG(EnableVtkOutput);
+
+/*!
+ * \brief Determines if the VTK output is written to disk asynchronously
+ *
+ * I.e. written to disk using a separate thread. This has only an effect if
+ * EnableVtkOutput is true and if the simulation is run sequentially. The reasons for
+ * this not being used for MPI-parallel simulations are that Dune's VTK output code does
+ * not support multi-threaded multi-process VTK output and even if it would, the result
+ * would be slower than when using synchronous output.
+ */
+NEW_PROP_TAG(EnableAsyncVtkOutput);
+
+/*!
+ * \brief Specify the format the VTK output is written to disk
+ *
+ * Possible values are:
+ *   - Dune::VTK::ascii (default)
+ *   - Dune::VTK::base64
+ *   - Dune::VTK::appendedraw
+ *   - Dune::VTK::appendedbase64
+ */
+NEW_PROP_TAG(VtkOutputFormat);
+
+//! Specify whether the some degrees of fredom can be constraint
+NEW_PROP_TAG(EnableConstraints);
+
+/*!
+ * \brief Specify the maximum size of a time integration [s].
+ *
+ * The default is to not limit the step size.
+ */
+NEW_PROP_TAG(MaxTimeStepSize);
+
+/*!
+ * \brief Specify the minimal size of a time integration [s].
+ *
+ * The default is to not limit the step size.
+ */
+NEW_PROP_TAG(MinTimeStepSize);
+
+/*!
+ * \brief The maximum allowed number of timestep divisions for the
+ *        Newton solver.
+ */
+NEW_PROP_TAG(MaxTimeStepDivisions);
+
+/*!
+ * \brief Continue with a non-converged solution instead of giving up
+ *        if we encounter a time step size smaller than the minimum time
+ *        step size.
+ */
+NEW_PROP_TAG(ContinueOnConvergenceError);
+
+/*!
+ * \brief Specify whether all intensive quantities for the grid should be
+ *        cached in the discretization.
+ *
+ * This potentially reduces the CPU time, but comes at the cost of
+ * higher memory consumption. In turn, the higher memory requirements
+ * may cause the simulation to exhibit worse cache coherence behavior
+ * which eats some of the computational benefits again.
+ */
+NEW_PROP_TAG(EnableIntensiveQuantityCache);
+
+/*!
+ * \brief Specify whether the storage terms for previous solutions should be cached.
+ *
+ * This potentially reduces the CPU time, but comes at the cost of higher memory
+ * consumption.
+ */
+NEW_PROP_TAG(EnableStorageCache);
+
+/*!
+ * \brief Specify whether to use the already calculated solutions as
+ *        starting values of the intensive quantities.
+ *
+ * This only makes sense if the calculation of the intensive quantities is
+ * very expensive (e.g. for non-linear fugacity functions where the
+ * solver converges faster).
+ */
+NEW_PROP_TAG(EnableThermodynamicHints);
+
+// mappers from local to global DOF indices
+
+/*!
+ * \brief The mapper to find the global index of a vertex.
+ */
+NEW_PROP_TAG(VertexMapper);
+
+/*!
+ * \brief The mapper to find the global index of an element.
+ */
+NEW_PROP_TAG(ElementMapper);
+
+/*!
+ * \brief The mapper to find the global index of a degree of freedom.
+ */
+NEW_PROP_TAG(DofMapper);
+
+/*!
+ * \brief The class which marks the border indices associated with the
+ *        degrees of freedom on a process boundary.
+ *
+ * This is required for the algebraic overlap stuff.
+ */
+NEW_PROP_TAG(BorderListCreator);
+
+/*!
+ * \brief The history size required by the time discretization
+ */
+NEW_PROP_TAG(TimeDiscHistorySize);
+
+/*!
+ * \brief Specify whether the storage terms use extensive quantities or not.
+ *
+ * Most models don't need this, but the (Navier-)Stokes ones do...
+ */
+NEW_PROP_TAG(ExtensiveStorageTerm);
+
+//! \brief Specify whether to use volumetric residuals or not
+NEW_PROP_TAG(UseVolumetricResidual);
+
+
+//! Specify if experimental features should be enabled or not.
+NEW_PROP_TAG(EnableExperiments);
+
+END_PROPERTIES
+
+#endif

opm/models/discretization/common/restrictprolong.hh

@@ -0,0 +1,205 @@
+// -*- 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.
+*/
+#ifndef EWOMS_COPYRESTRICTPROLONG_HH
+#define EWOMS_COPYRESTRICTPROLONG_HH
+
+#include <opm/material/common/Unused.hpp>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/space/common/restrictprolonginterface.hh>
+#endif
+
+namespace Opm
+{
+    template < class Grid, class Container >
+    class CopyRestrictProlong;
+
+    template < class Grid, class Container >
+    struct CopyRestrictProlongTraits
+    {
+      typedef typename Grid::ctype DomainFieldType;
+      typedef CopyRestrictProlong< Grid, Container >  RestProlImp;
+    };
+
+    template< class Grid, class Container >
+    class CopyRestrictProlong
+#if HAVE_DUNE_FEM
+    : public Dune::Fem::RestrictProlongInterfaceDefault< CopyRestrictProlongTraits< Grid, Container > >
+#endif
+    {
+      typedef CopyRestrictProlong< Grid, Container > ThisType;
+
+      Container& container_;
+    public:
+      typedef typename Grid::ctype DomainFieldType;
+
+      explicit CopyRestrictProlong( Container& container )
+        : container_( container )
+      {}
+
+      /** \brief explicit set volume ratio of son and father
+       *
+       *  \param[in]  weight  volume of son / volume of father
+       *
+       *  \note If this ratio is set, it is assume to be constant.
+       */
+      template <class Field>
+      void setFatherChildWeight(const Field& weight OPM_UNUSED) const
+      {}
+
+      //! restrict data to father
+      template< class Entity >
+      void restrictLocal ( const Entity& father, const Entity& son, bool initialize ) const
+      {
+        container_.resize();
+        assert( container_.codimension() == 0 );
+        if( initialize )
+        {
+          // copy values from son to father (only once)
+          container_[ father ] = container_[ son ];
+        }
+      }
+
+      //! restrict data to father
+      template< class Entity, class LocalGeometry >
+      void restrictLocal(const Entity& father,
+                         const Entity& son,
+                         const LocalGeometry& geometryInFather OPM_UNUSED,
+                         bool initialize) const
+      { restrictLocal(father, son, initialize); }
+
+      //! prolong data to children
+      template< class Entity >
+      void prolongLocal(const Entity& father,
+                        const Entity& son,
+                        bool initialize OPM_UNUSED) const
+      {
+        container_.resize();
+        assert( container_.codimension() == 0 );
+        // copy values from father to son all sons
+        container_[ son ] = container_[ father ];
+      }
+
+      //! prolong data to children
+      template< class Entity, class LocalGeometry >
+      void prolongLocal(const Entity& father,
+                        const Entity& son,
+                        const LocalGeometry& geometryInFather OPM_UNUSED,
+                        bool initialize) const
+      { prolongLocal(father, son, initialize); }
+
+      /** \brief add discrete function to communicator
+       *  \param[in]  comm  Communicator to add the discrete functions to
+       */
+      template< class Communicator >
+      void addToList(Communicator& comm OPM_UNUSED)
+      {
+        // TODO
+      }
+
+      /** \brief add discrete function to load balancer
+       *  \param[in]  lb LoadBalancer to add the discrete functions to
+       */
+      template< class LoadBalancer >
+      void addToLoadBalancer(LoadBalancer& lb OPM_UNUSED)
+      {
+        // TODO
+      }
+
+    };
+
+
+    class EmptyRestrictProlong;
+
+    struct EmptyRestrictProlongTraits
+    {
+      typedef double                DomainFieldType;
+      typedef EmptyRestrictProlong  RestProlImp;
+    };
+
+    class EmptyRestrictProlong
+#if HAVE_DUNE_FEM
+    : public Dune::Fem::RestrictProlongInterfaceDefault< EmptyRestrictProlongTraits >
+#endif
+    {
+      typedef EmptyRestrictProlong ThisType;
+
+    public:
+      /** \brief explicit set volume ratio of son and father
+       *
+       *  \param[in]  weight  volume of son / volume of father
+       *
+       *  \note If this ratio is set, it is assume to be constant.
+       */
+      template <class Field>
+      void setFatherChildWeight(const Field& weight OPM_UNUSED) const
+      { }
+
+      //! restrict data to father
+      template< class Entity >
+      void restrictLocal(const Entity& father OPM_UNUSED,
+                         const Entity& son OPM_UNUSED,
+                         bool initialize OPM_UNUSED) const
+      { }
+
+      //! restrict data to father
+      template< class Entity, class LocalGeometry >
+      void restrictLocal(const Entity& father OPM_UNUSED,
+                         const Entity& son OPM_UNUSED,
+                         const LocalGeometry& geometryInFather OPM_UNUSED,
+                         bool initialize OPM_UNUSED) const
+      { }
+
+      //! prolong data to children
+      template< class Entity >
+      void prolongLocal(const Entity& father OPM_UNUSED,
+                        const Entity& son OPM_UNUSED,
+                        bool initialize OPM_UNUSED) const
+      { }
+
+      //! prolong data to children
+      template< class Entity, class LocalGeometry >
+      void prolongLocal(const Entity& father OPM_UNUSED,
+                        const Entity& son OPM_UNUSED,
+                        const LocalGeometry& geometryInFather OPM_UNUSED,
+                        bool initialize OPM_UNUSED) const
+      { }
+
+      /** \brief add discrete function to communicator
+       *  \param[in]  comm  Communicator to add the discrete functions to
+       */
+      template< class Communicator >
+      void addToList(Communicator& comm OPM_UNUSED)
+      { }
+
+      /** \brief add discrete function to load balancer
+       *  \param[in]  lb LoadBalancer to add the discrete functions to
+       */
+      template< class LoadBalancer >
+      void addToLoadBalancer(LoadBalancer& lb OPM_UNUSED)
+      { }
+    };
+
+} // namespace Opm
+
+#endif

opm/models/discretization/ecfv/ecfvbaseoutputmodule.hh

@@ -0,0 +1,81 @@
+// -*- 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::EcfvBaseOutputModule
+ */
+#ifndef EWOMS_ECFV_VTK_BASE_OUTPUT_MODULE_HH
+#define EWOMS_ECFV_VTK_BASE_OUTPUT_MODULE_HH
+
+#include "ecfvproperties.hh"
+
+#include <ewoms/io/vtkmultiwriter.hh>
+
+namespace Opm {
+/*!
+ * \ingroup EcfvDiscretization
+ *
+ * \brief Implements the discretization specific parts of writing files.
+ */
+template<class TypeTag>
+class EcfvBaseOutputModule
+{
+public:
+    typedef BaseOutputWriter::Scalar Scalar;
+    typedef BaseOutputWriter::Vector Vector;
+    typedef BaseOutputWriter::ScalarBuffer ScalarBuffer;
+    typedef BaseOutputWriter::VectorBuffer VectorBuffer;
+    typedef BaseOutputWriter::TensorBuffer TensorBuffer;
+
+    /*!
+     * \brief Add a buffer where the data is associated with the
+     *        degrees of freedom to the current VTK output file.
+     */
+    static void attachScalarDofData_(BaseOutputWriter& baseWriter,
+                                     ScalarBuffer& buffer,
+                                     const std::string& name)
+    { baseWriter.attachScalarElementData(buffer, name.c_str()); }
+
+    /*!
+     * \brief Add a buffer where the data is associated with the
+     *        degrees of freedom to the current VTK output file.
+     */
+    static void attachVectorDofData_(BaseOutputWriter& baseWriter,
+                                     VectorBuffer& buffer,
+                                     const std::string& name)
+    { baseWriter.attachVectorElementData(buffer, name.c_str()); }
+
+    /*!
+     * \brief Add a buffer where the data is associated with the
+     *        degrees of freedom to the current VTK output file.
+     */
+    static void attachTensorDofData_(BaseOutputWriter& baseWriter,
+                                     TensorBuffer& buffer,
+                                     const std::string& name)
+    { baseWriter.attachTensorElementData(buffer, name.c_str()); }
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/ecfv/ecfvdiscretization.hh

@@ -0,0 +1,216 @@
+// -*- 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::EcfvDiscretization
+ */
+#ifndef EWOMS_ECFV_DISCRETIZATION_HH
+#define EWOMS_ECFV_DISCRETIZATION_HH
+
+#include <opm/material/densead/Math.hpp>
+
+#include "ecfvproperties.hh"
+#include "ecfvstencil.hh"
+#include "ecfvgridcommhandlefactory.hh"
+#include "ecfvbaseoutputmodule.hh"
+
+#include <opm/simulators/linalg/elementborderlistfromgrid.hh>
+#include <opm/models/discretization/common/fvbasediscretization.hh>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/space/common/functionspace.hh>
+#include <dune/fem/space/finitevolume.hh>
+#endif
+
+namespace Opm {
+template <class TypeTag>
+class EcfvDiscretization;
+}
+
+BEGIN_PROPERTIES
+
+//! Set the stencil
+SET_PROP(EcfvDiscretization, Stencil)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+public:
+    typedef Opm::EcfvStencil<Scalar, GridView> type;
+};
+
+//! Mapper for the degrees of freedoms.
+SET_TYPE_PROP(EcfvDiscretization, DofMapper, typename GET_PROP_TYPE(TypeTag, ElementMapper));
+
+//! The concrete class which manages the spatial discretization
+SET_TYPE_PROP(EcfvDiscretization, Discretization, Opm::EcfvDiscretization<TypeTag>);
+
+//! The base class for the output modules (decides whether to write
+//! element or vertex based fields)
+SET_TYPE_PROP(EcfvDiscretization, DiscBaseOutputModule,
+              Opm::EcfvBaseOutputModule<TypeTag>);
+
+//! The class to create grid communication handles
+SET_TYPE_PROP(EcfvDiscretization, GridCommHandleFactory,
+              Opm::EcfvGridCommHandleFactory<TypeTag>);
+
+#if HAVE_DUNE_FEM
+//! Set the DiscreteFunctionSpace
+SET_PROP(EcfvDiscretization, DiscreteFunctionSpace)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar)   Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridPart) GridPart;
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    typedef Dune::Fem::FunctionSpace<typename GridPart::GridType::ctype,
+                                     Scalar,
+                                     GridPart::GridType::dimensionworld,
+                                     numEq> FunctionSpace;
+public:
+    typedef Dune::Fem::FiniteVolumeSpace< FunctionSpace, GridPart, 0 > type;
+};
+#endif
+
+//! Set the border list creator for to the one of an element based
+//! method
+SET_PROP(EcfvDiscretization, BorderListCreator)
+{ private:
+    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+public:
+    typedef Opm::Linear::ElementBorderListFromGrid<GridView, ElementMapper> type;
+};
+
+//! For the element centered finite volume method, ghost and overlap elements must be
+//! assembled to calculate the fluxes over the process boundary faces of the local
+//! process' grid partition
+SET_BOOL_PROP(EcfvDiscretization, LinearizeNonLocalElements, true);
+
+//! locking is not required for the element centered finite volume method because race
+//! conditions cannot occur since each matrix/vector entry is written exactly once
+SET_BOOL_PROP(EcfvDiscretization, UseLinearizationLock, false);
+
+END_PROPERTIES
+
+namespace Opm {
+/*!
+ * \ingroup EcfvDiscretization
+ *
+ * \brief The base class for the element-centered finite-volume discretization scheme.
+ */
+template<class TypeTag>
+class EcfvDiscretization : public FvBaseDiscretization<TypeTag>
+{
+    typedef FvBaseDiscretization<TypeTag> ParentType;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+
+public:
+    EcfvDiscretization(Simulator& simulator)
+        : ParentType(simulator)
+    { }
+
+    /*!
+     * \brief Returns a string of discretization's human-readable name
+     */
+    static std::string discretizationName()
+    { return "ecfv"; }
+
+    /*!
+     * \brief Returns the number of global degrees of freedom (DOFs) due to the grid
+     */
+    size_t numGridDof() const
+    { return static_cast<size_t>(this->gridView_.size(/*codim=*/0)); }
+
+    /*!
+     * \brief Mapper to convert the Dune entities of the
+     *        discretization's degrees of freedoms are to indices.
+     */
+    const DofMapper& dofMapper() const
+    { return this->elementMapper(); }
+
+    /*!
+     * \brief Syncronize the values of the primary variables on the
+     *        degrees of freedom that overlap with the neighboring
+     *        processes.
+     *
+     * For the Element Centered Finite Volume discretization, this
+     * method retrieves the primary variables corresponding to
+     * overlap/ghost elements from their respective master process.
+     */
+    void syncOverlap()
+    {
+        // syncronize the solution on the ghost and overlap elements
+        typedef GridCommHandleGhostSync<PrimaryVariables,
+                                        SolutionVector,
+                                        DofMapper,
+                                        /*commCodim=*/0> GhostSyncHandle;
+
+        auto ghostSync = GhostSyncHandle(this->solution(/*timeIdx=*/0),
+                                         asImp_().dofMapper());
+        this->gridView().communicate(ghostSync,
+                                     Dune::InteriorBorder_All_Interface,
+                                     Dune::ForwardCommunication);
+    }
+
+    /*!
+     * \brief Serializes the current state of the model.
+     *
+     * \tparam Restarter The type of the serializer class
+     *
+     * \param res The serializer object
+     */
+    template <class Restarter>
+    void serialize(Restarter& res)
+    { res.template serializeEntities</*codim=*/0>(asImp_(), this->gridView_); }
+
+    /*!
+     * \brief Deserializes the state of the model.
+     *
+     * \tparam Restarter The type of the serializer class
+     *
+     * \param res The serializer object
+     */
+    template <class Restarter>
+    void deserialize(Restarter& res)
+    {
+        res.template deserializeEntities</*codim=*/0>(asImp_(), this->gridView_);
+        this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0);
+    }
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+} // namespace Opm
+
+#endif

opm/models/discretization/ecfv/ecfvgridcommhandlefactory.hh

@@ -0,0 +1,88 @@
+// -*- 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::EcfvGridCommHandleFactory
+ */
+#ifndef EWOMS_ECFV_GRID_COMM_HANDLE_FACTORY_HH
+#define EWOMS_ECFV_GRID_COMM_HANDLE_FACTORY_HH
+
+#include "ecfvproperties.hh"
+
+#include <opm/models/parallel/gridcommhandles.hh>
+
+namespace Opm {
+/*!
+ * \ingroup EcfvDiscretization
+ *
+ * \brief A class which provides types for DUNE grid handles for
+ *        communication.
+ *
+ * This is required for parallel computations
+ */
+template<class TypeTag>
+class EcfvGridCommHandleFactory
+{
+    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+
+public:
+    /*!
+     * \brief Return a handle which computes the minimum of a value
+     *        for each overlapping degree of freedom across all processes.
+     */
+    template <class ValueType, class ArrayType>
+    static std::shared_ptr<GridCommHandleMin<ValueType, ArrayType,  DofMapper, /*commCodim=*/0> >
+    minHandle(ArrayType& array, const DofMapper& dofMapper)
+    {
+        typedef GridCommHandleMin<ValueType, ArrayType,  DofMapper, /*commCodim=*/0> Handle;
+        return  std::shared_ptr<Handle>(new Handle(array, dofMapper));
+    }
+
+    /*!
+     * \brief Return a handle which computes the maximum of a value
+     *        for each overlapping degree of freedom across all processes.
+     */
+    template <class ValueType, class ArrayType>
+    static std::shared_ptr<GridCommHandleMax<ValueType, ArrayType,  DofMapper, /*commCodim=*/0> >
+    maxHandle(ArrayType& array, const DofMapper& dofMapper)
+    {
+        typedef GridCommHandleMax<ValueType, ArrayType,  DofMapper, /*commCodim=*/0> Handle;
+        return  std::shared_ptr<Handle>(new Handle(array, dofMapper));
+    }
+
+    /*!
+     * \brief Return a handle which computes the sum of all values
+     *        all overlapping degrees of freedom across all processes.
+     */
+    template <class ValueType, class ArrayType>
+    static std::shared_ptr<GridCommHandleSum<ValueType, ArrayType,  DofMapper, /*commCodim=*/0> >
+    sumHandle(ArrayType& array, const DofMapper& dofMapper)
+    {
+        typedef GridCommHandleSum<ValueType, ArrayType,  DofMapper, /*commCodim=*/0> Handle;
+        return  std::shared_ptr<Handle>(new Handle(array, dofMapper));
+    }
+};
+} // namespace Opm
+
+#endif

opm/models/discretization/ecfv/ecfvproperties.hh

@@ -0,0 +1,43 @@
+// -*- 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
+ *
+ * \brief Declare the basic properties used by the common infrastructure of
+ *        the element-centered finite volume discretization.
+ */
+#ifndef EWOMS_ECFV_PROPERTIES_HH
+#define EWOMS_ECFV_PROPERTIES_HH
+
+#include <opm/models/discretization/common/fvbaseproperties.hh>
+
+#include <opm/models/utils/propertysystem.hh>
+
+BEGIN_PROPERTIES
+
+//! The type tag for models based on the ECFV-scheme
+NEW_TYPE_TAG(EcfvDiscretization, INHERITS_FROM(FvBaseDiscretization));
+
+END_PROPERTIES
+
+#endif

opm/models/discretization/ecfv/ecfvstencil.hh

@@ -0,0 +1,399 @@
+// -*- 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::EcfvStencil
+ */
+#ifndef EWOMS_ECFV_STENCIL_HH
+#define EWOMS_ECFV_STENCIL_HH
+
+#include <opm/models/utils/quadraturegeometries.hh>
+
+#include <opm/material/common/ConditionalStorage.hpp>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/grid/common/intersectioniterator.hh>
+#include <dune/geometry/type.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/version.hh>
+
+#include <vector>
+
+namespace Opm {
+/*!
+ * \ingroup EcfvDiscretization
+ *
+ * \brief Represents the stencil (finite volume geometry) of a single
+ *        element in the ECFV discretization.
+ *
+ * The ECFV discretization is a element centered finite volume
+ * approach. This means that each element corresponds to a control
+ * volume.
+ */
+template <class Scalar,
+          class GridView,
+          bool needFaceIntegrationPos = true,
+          bool needFaceNormal = true>
+class EcfvStencil
+{
+    enum { dimWorld = GridView::dimensionworld };
+
+    typedef typename GridView::ctype CoordScalar;
+    typedef typename GridView::Intersection Intersection;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+#if DUNE_VERSION_NEWER(DUNE_GRID, 2,6)
+    typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView> ElementMapper;
+#else
+    typedef Dune::MultipleCodimMultipleGeomTypeMapper<GridView, Dune::MCMGElementLayout> ElementMapper;
+#endif
+
+    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
+
+    typedef Dune::FieldVector<Scalar, dimWorld> WorldVector;
+
+public:
+    typedef Element        Entity;
+    typedef ElementMapper  Mapper;
+
+    typedef typename Element::Geometry LocalGeometry;
+
+    /*!
+     * \brief Represents a sub-control volume.
+     *
+     * For element centered finite volumes, this is equivalent to the
+     * element, in the vertex centered finite volume approach, this
+     * corresponds to the intersection of a finite volume and the
+     * grid element.
+     */
+    class SubControlVolume
+    {
+    public:
+        // default construct an uninitialized object.
+        // this is only here because std::vector needs it...
+        SubControlVolume()
+        {}
+
+        SubControlVolume(const Element& element)
+            : element_(element)
+        { update(); }
+
+        void update(const Element& element)
+        { element_ = element; }
+
+        void update()
+        {
+            const auto& geometry = element_.geometry();
+            centerPos_ = geometry.center();
+            volume_ = geometry.volume();
+        }
+
+        /*!
+         * \brief The global position associated with the sub-control volume
+         */
+        const GlobalPosition& globalPos() const
+        { return centerPos_; }
+
+        /*!
+         * \brief The center of the sub-control volume
+         */
+        const GlobalPosition& center() const
+        { return centerPos_; }
+
+        /*!
+         * \brief The volume [m^3] occupied by the sub-control volume
+         */
+        Scalar volume() const
+        { return volume_; }
+
+        /*!
+         * \brief The geometry of the sub-control volume.
+         */
+        const LocalGeometry geometry() const
+        { return element_.geometry(); }
+
+        /*!
+         * \brief Geometry of the sub-control volume relative to parent.
+         */
+        const LocalGeometry localGeometry() const
+        { return element_.geometryInFather(); }
+
+    private:
+        GlobalPosition centerPos_;
+        Scalar volume_;
+        Element element_;
+    };
+
+    /*!
+     * \brief Represents a face of a sub-control volume.
+     */
+    template <bool needIntegrationPos, bool needNormal>
+    class EcfvSubControlVolumeFace
+    {
+    public:
+        EcfvSubControlVolumeFace()
+        {}
+
+        EcfvSubControlVolumeFace(const Intersection& intersection, unsigned localNeighborIdx)
+        {
+            exteriorIdx_ = static_cast<unsigned short>(localNeighborIdx);
+
+            if (needNormal)
+                (*normal_) = intersection.centerUnitOuterNormal();
+
+            const auto& geometry = intersection.geometry();
+            if (needIntegrationPos)
+                (*integrationPos_) = geometry.center();
+            area_ = geometry.volume();
+        }
+
+        /*!
+         * \brief Returns the local index of the degree of freedom to
+         *        the face's interior.
+         */
+        unsigned short interiorIndex() const
+        {
+            // The local index of the control volume in the interior
+            // of a face of the stencil in the element centered finite
+            // volume discretization is always the "central"
+            // element. In this implementation this element always has
+            // index 0....
+            return 0;
+        }
+
+        /*!
+         * \brief Returns the local index of the degree of freedom to
+         *        the face's outside.
+         */
+        unsigned short exteriorIndex() const
+        { return exteriorIdx_; }
+
+        /*!
+         * \brief Returns the global position of the face's
+         *        integration point.
+         */
+        const GlobalPosition& integrationPos() const
+        { return *integrationPos_; }
+
+        /*!
+         * \brief Returns the outer unit normal at the face's
+         *        integration point.
+         */
+        const WorldVector& normal() const
+        { return *normal_; }
+
+        /*!
+         * \brief Returns the area [m^2] of the face
+         */
+        Scalar area() const
+        { return area_; }
+
+    private:
+        Opm::ConditionalStorage<needIntegrationPos, GlobalPosition> integrationPos_;
+        Opm::ConditionalStorage<needNormal, WorldVector> normal_;
+        Scalar area_;
+
+        unsigned short exteriorIdx_;
+    };
+
+    typedef EcfvSubControlVolumeFace<needFaceIntegrationPos, needFaceNormal> SubControlVolumeFace;
+    typedef EcfvSubControlVolumeFace</*needFaceIntegrationPos=*/true, needFaceNormal> BoundaryFace;
+
+    EcfvStencil(const GridView& gridView, const Mapper& mapper)
+        : gridView_(gridView)
+        , elementMapper_(mapper)
+    {
+        // try to ensure that the mapper passed indeed maps elements
+        assert(int(gridView.size(/*codim=*/0)) == int(elementMapper_.size()));
+    }
+
+    void updateTopology(const Element& element)
+    {
+        auto isIt = gridView_.ibegin(element);
+        const auto& endIsIt = gridView_.iend(element);
+
+        // add the "center" element of the stencil
+        subControlVolumes_.clear();
+        subControlVolumes_.emplace_back(/*SubControlVolume(*/element/*)*/);
+        elements_.clear();
+        elements_.emplace_back(element);
+
+        interiorFaces_.clear();
+        boundaryFaces_.clear();
+
+        for (; isIt != endIsIt; ++isIt) {
+            const auto& intersection = *isIt;
+            // if the current intersection has a neighbor, add a
+            // degree of freedom and an internal face, else add a
+            // boundary face
+            if (intersection.neighbor()) {
+                elements_.emplace_back( intersection.outside() );
+                subControlVolumes_.emplace_back(/*SubControlVolume(*/elements_.back()/*)*/);
+                interiorFaces_.emplace_back(/*SubControlVolumeFace(*/intersection, subControlVolumes_.size() - 1/*)*/);
+            }
+            else {
+                boundaryFaces_.emplace_back(/*SubControlVolumeFace(*/intersection, - 10000/*)*/);
+            }
+        }
+    }
+
+    void updatePrimaryTopology(const Element& element)
+    {
+        // add the "center" element of the stencil
+        subControlVolumes_.clear();
+        subControlVolumes_.emplace_back(/*SubControlVolume(*/element/*)*/);
+        elements_.clear();
+        elements_.emplace_back(element);
+    }
+
+    void update(const Element& element)
+    {
+        updateTopology(element);
+    }
+
+    void updateCenterGradients()
+    {
+        assert(false); // not yet implemented
+    }
+
+    /*!
+     * \brief Return the element to which the stencil refers.
+     */
+    const Element& element() const
+    { return element( 0 ); }
+
+    /*!
+     * \brief Return the grid view of the element to which the stencil
+     *        refers.
+     */
+    const GridView& gridView() const
+    { return *gridView_; }
+
+    /*!
+     * \brief Returns the number of degrees of freedom which the
+     *        current element interacts with.
+     */
+    size_t numDof() const
+    { return subControlVolumes_.size(); }
+
+    /*!
+     * \brief Returns the number of degrees of freedom which are contained
+     *        by within the current element.
+     *
+     * Primary DOFs are always expected to have a lower index than
+     * "secondary" DOFs.
+     *
+     * For element centered finite elements, this is only the central DOF.
+     */
+    size_t numPrimaryDof() const
+    { return 1; }
+
+    /*!
+     * \brief Return the global space index given the index of a degree of
+     *        freedom.
+     */
+    unsigned globalSpaceIndex(unsigned dofIdx) const
+    {
+        assert(0 <= dofIdx && dofIdx < numDof());
+
+        return static_cast<unsigned>(elementMapper_.index(element(dofIdx)));
+    }
+
+    /*!
+     * \brief Return partition type of a given degree of freedom
+     */
+    Dune::PartitionType partitionType(unsigned dofIdx) const
+    { return elements_[dofIdx]->partitionType(); }
+
+    /*!
+     * \brief Return the element given the index of a degree of
+     *        freedom.
+     *
+     * If no degree of freedom index is passed, the element which was
+     * passed to the update() method is returned...
+     */
+    const Element& element(unsigned dofIdx) const
+    {
+        assert(0 <= dofIdx && dofIdx < numDof());
+
+        return elements_[dofIdx];
+    }
+
+    /*!
+     * \brief Return the entity given the index of a degree of
+     *        freedom.
+     */
+    const Entity& entity(unsigned dofIdx) const
+    {
+        return element( dofIdx );
+    }
+
+    /*!
+     * \brief Returns the sub-control volume object belonging to a
+     *        given degree of freedom.
+     */
+    const SubControlVolume& subControlVolume(unsigned dofIdx) const
+    { return subControlVolumes_[dofIdx]; }
+
+    /*!
+     * \brief Returns the number of interior faces of the stencil.
+     */
+    size_t numInteriorFaces() const
+    { return interiorFaces_.size(); }
+
+    /*!
+     * \brief Returns the face object belonging to a given face index
+     *        in the interior of the domain.
+     */
+    const SubControlVolumeFace& interiorFace(unsigned faceIdx) const
+    { return interiorFaces_[faceIdx]; }
+
+    /*!
+     * \brief Returns the number of boundary faces of the stencil.
+     */
+    size_t numBoundaryFaces() const
+    { return boundaryFaces_.size(); }
+
+    /*!
+     * \brief Returns the boundary face object belonging to a given
+     *        boundary face index.
+     */
+    const BoundaryFace& boundaryFace(unsigned bfIdx) const
+    { return boundaryFaces_[bfIdx]; }
+
+protected:
+    const GridView&       gridView_;
+    const ElementMapper&  elementMapper_;
+
+    std::vector<Element> elements_;
+    std::vector<SubControlVolume>      subControlVolumes_;
+    std::vector<SubControlVolumeFace>  interiorFaces_;
+    std::vector<BoundaryFace>  boundaryFaces_;
+};
+
+} // namespace Opm
+
+
+#endif
+

opm/models/discretization/vcfv/p1fegradientcalculator.hh

@@ -0,0 +1,365 @@
+// -*- 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::P1FeGradientCalculator
+ */
+#ifndef EWOMS_P1FE_GRADIENT_CALCULATOR_HH
+#define EWOMS_P1FE_GRADIENT_CALCULATOR_HH
+
+#include "vcfvproperties.hh"
+
+#include <opm/models/discretization/common/fvbasegradientcalculator.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/version.hh>
+
+#include <dune/geometry/type.hh>
+
+#if HAVE_DUNE_LOCALFUNCTIONS
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+#endif // HAVE_DUNE_LOCALFUNCTIONS
+
+#include <dune/common/fvector.hh>
+
+#include <vector>
+
+#ifdef HAVE_DUNE_LOCALFUNCTIONS
+#define EWOMS_NO_LOCALFUNCTIONS_UNUSED
+#else
+#define EWOMS_NO_LOCALFUNCTIONS_UNUSED  OPM_UNUSED
+#endif
+
+namespace Opm {
+/*!
+ * \ingroup FiniteElementDiscretizations
+ *
+ * \brief This class calculates gradients of arbitrary quantities at flux integration
+ *        points using first order finite elemens ansatz functions.
+ *
+ * This approach can also be used for the vertex-centered finite volume (VCFV)
+ * discretization.
+ */
+template<class TypeTag>
+class P1FeGradientCalculator : public FvBaseGradientCalculator<TypeTag>
+{
+    typedef FvBaseGradientCalculator<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+    enum { dim = GridView::dimension };
+
+    // set the maximum number of degrees of freedom and the maximum
+    // number of flux approximation points per elements. For this, we
+    // assume cubes as the type of element with the most vertices...
+    enum { maxDof = (2 << dim) };
+    enum { maxFap = maxDof };
+
+    typedef typename GridView::ctype CoordScalar;
+    typedef Dune::FieldVector<Scalar, dim> DimVector;
+
+#if HAVE_DUNE_LOCALFUNCTIONS
+    typedef Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1> LocalFiniteElementCache;
+    typedef typename LocalFiniteElementCache::FiniteElementType LocalFiniteElement;
+    typedef typename LocalFiniteElement::Traits::LocalBasisType::Traits LocalBasisTraits;
+    typedef typename LocalBasisTraits::JacobianType ShapeJacobian;
+#endif // HAVE_DUNE_LOCALFUNCTIONS
+
+public:
+    /*!
+     * \brief Precomputes the common values to calculate gradients and
+     *        values of quantities at any flux approximation point.
+     *
+     * \param elemCtx The current execution context
+     */
+    template <bool prepareValues = true, bool prepareGradients = true>
+    void prepare(const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                 unsigned timeIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED)
+    {
+        if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
+#if !HAVE_DUNE_LOCALFUNCTIONS
+            // The dune-localfunctions module is required for P1 finite element gradients
+            throw std::logic_error("The dune-localfunctions module is required in oder to use"
+                                   " finite element gradients");
+#else
+            const auto& stencil = elemCtx.stencil(timeIdx);
+
+            const LocalFiniteElement& localFE = feCache_.get(elemCtx.element().type());
+            localFiniteElement_ = &localFE;
+
+            // loop over all face centeres
+            for (unsigned faceIdx = 0; faceIdx < stencil.numInteriorFaces(); ++faceIdx) {
+                const auto& localFacePos = stencil.interiorFace(faceIdx).localPos();
+
+                // Evaluate the P1 shape functions and their gradients at all
+                // flux approximation points.
+                if (prepareValues)
+                    localFE.localBasis().evaluateFunction(localFacePos, p1Value_[faceIdx]);
+
+                if (prepareGradients) {
+                    // first, get the shape function's gradient in local coordinates
+                    std::vector<ShapeJacobian> localGradient;
+                    localFE.localBasis().evaluateJacobian(localFacePos, localGradient);
+
+                    // convert to a gradient in global space by
+                    // multiplying with the inverse transposed jacobian of
+                    // the position
+                    const auto& geom = elemCtx.element().geometry();
+                    const auto& jacInvT =
+                        geom.jacobianInverseTransposed(localFacePos);
+
+                    size_t numVertices = elemCtx.numDof(timeIdx);
+                    for (unsigned vertIdx = 0; vertIdx < numVertices; vertIdx++) {
+                        jacInvT.mv(/*xVector=*/localGradient[vertIdx][0],
+                                   /*destVector=*/p1Gradient_[faceIdx][vertIdx]);
+                    }
+                }
+            }
+#endif
+        }
+        else
+            ParentType::template prepare<prepareValues, prepareGradients>(elemCtx, timeIdx);
+    }
+
+    /*!
+     * \brief Calculates the value of an arbitrary quantity at any
+     *        interior flux approximation point.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point
+     *               in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    auto calculateScalarValue(const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                              unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                              const QuantityCallback& quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const
+        ->  typename std::remove_reference<typename QuantityCallback::ResultType>::type
+    {
+        if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
+#if !HAVE_DUNE_LOCALFUNCTIONS
+            // The dune-localfunctions module is required for P1 finite element gradients
+            throw std::logic_error("The dune-localfunctions module is required in oder to use"
+                                   " finite element gradients");
+#else
+            typedef typename std::remove_reference<typename QuantityCallback::ResultType>::type QuantityConstType;
+            typedef typename std::remove_const<QuantityConstType>::type QuantityType;
+            typedef Opm::MathToolbox<QuantityType> Toolbox;
+
+            // If the user does not want to use two-point gradients, we
+            // use P1 finite element gradients..
+            QuantityType value(0.0);
+            for (unsigned vertIdx = 0; vertIdx < elemCtx.numDof(/*timeIdx=*/0); ++vertIdx) {
+                if (std::is_same<QuantityType, Scalar>::value ||
+                    elemCtx.focusDofIndex() == vertIdx)
+                    value += quantityCallback(vertIdx)*p1Value_[fapIdx][vertIdx];
+                else
+                    value += Toolbox::value(quantityCallback(vertIdx))*p1Value_[fapIdx][vertIdx];
+            }
+
+            return value;
+#endif
+        }
+        else
+            return ParentType::calculateScalarValue(elemCtx, fapIdx, quantityCallback);
+    }
+
+    /*!
+     * \brief Calculates the value of an arbitrary quantity at any
+     *        interior flux approximation point.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point
+     *               in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    auto calculateVectorValue(const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                              unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                              const QuantityCallback& quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const
+        ->  typename std::remove_reference<typename QuantityCallback::ResultType>::type
+    {
+        if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
+#if !HAVE_DUNE_LOCALFUNCTIONS
+            // The dune-localfunctions module is required for P1 finite element gradients
+            throw std::logic_error("The dune-localfunctions module is required in oder to use"
+                                   " finite element gradients");
+#else
+            typedef typename std::remove_reference<typename QuantityCallback::ResultType>::type QuantityConstType;
+            typedef typename std::remove_const<QuantityConstType>::type QuantityType;
+
+            typedef decltype(std::declval<QuantityType>()[0]) RawFieldType;
+            typedef typename std::remove_const<typename std::remove_reference<RawFieldType>::type>::type FieldType;
+
+            typedef Opm::MathToolbox<FieldType> Toolbox;
+
+            // If the user does not want to use two-point gradients, we
+            // use P1 finite element gradients..
+            QuantityType value(0.0);
+            for (unsigned vertIdx = 0; vertIdx < elemCtx.numDof(/*timeIdx=*/0); ++vertIdx) {
+                if (std::is_same<QuantityType, Scalar>::value ||
+                    elemCtx.focusDofIndex() == vertIdx)
+                {
+                    const auto& tmp = quantityCallback(vertIdx);
+                    for (unsigned k = 0; k < tmp.size(); ++k)
+                        value[k] += tmp[k]*p1Value_[fapIdx][vertIdx];
+                }
+                else {
+                    const auto& tmp = quantityCallback(vertIdx);
+                    for (unsigned k = 0; k < tmp.size(); ++k)
+                        value[k] += Toolbox::value(tmp[k])*p1Value_[fapIdx][vertIdx];
+                }
+            }
+
+            return value;
+#endif
+        }
+        else
+            return ParentType::calculateVectorValue(elemCtx, fapIdx, quantityCallback);
+    }
+
+    /*!
+     * \brief Calculates the gradient of an arbitrary quantity at any
+     *        flux approximation point.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point
+     *               in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback, class EvalDimVector>
+    void calculateGradient(EvalDimVector& quantityGrad EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                           const ElementContext& elemCtx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                           unsigned fapIdx EWOMS_NO_LOCALFUNCTIONS_UNUSED,
+                           const QuantityCallback& quantityCallback EWOMS_NO_LOCALFUNCTIONS_UNUSED) const
+    {
+        if (GET_PROP_VALUE(TypeTag, UseP1FiniteElementGradients)) {
+#if !HAVE_DUNE_LOCALFUNCTIONS
+            // The dune-localfunctions module is required for P1 finite element gradients
+            throw std::logic_error("The dune-localfunctions module is required in oder to use"
+                                   " finite element gradients");
+#else
+            typedef typename std::remove_reference<typename QuantityCallback::ResultType>::type QuantityConstType;
+            typedef typename std::remove_const<QuantityConstType>::type QuantityType;
+
+            // If the user does not want two-point gradients, we use P1 finite element
+            // gradients...
+            quantityGrad = 0.0;
+            for (unsigned vertIdx = 0; vertIdx < elemCtx.numDof(/*timeIdx=*/0); ++vertIdx) {
+                if (std::is_same<QuantityType, Scalar>::value ||
+                    elemCtx.focusDofIndex() == vertIdx)
+                {
+                    const auto& dofVal = quantityCallback(vertIdx);
+                    const auto& tmp = p1Gradient_[fapIdx][vertIdx];
+                    for (int dimIdx = 0; dimIdx < dim; ++ dimIdx)
+                        quantityGrad[dimIdx] += dofVal*tmp[dimIdx];
+                }
+                else {
+                    const auto& dofVal = quantityCallback(vertIdx);
+                    const auto& tmp = p1Gradient_[fapIdx][vertIdx];
+                    for (int dimIdx = 0; dimIdx < dim; ++ dimIdx)
+                        quantityGrad[dimIdx] += Opm::scalarValue(dofVal)*tmp[dimIdx];
+                }
+            }
+#endif
+        }
+        else
+            ParentType::calculateGradient(quantityGrad, elemCtx, fapIdx, quantityCallback);
+    }
+
+    /*!
+     * \brief Calculates the value of an arbitrary quantity at any
+     *        flux approximation point on the grid boundary.
+     *
+     * Boundary values are always calculated using the two-point
+     * approximation.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point
+     *               in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback>
+    auto calculateBoundaryValue(const ElementContext& elemCtx,
+                                unsigned fapIdx,
+                                const QuantityCallback& quantityCallback)
+        ->  decltype(ParentType::calculateBoundaryValue(elemCtx, fapIdx, quantityCallback))
+    { return ParentType::calculateBoundaryValue(elemCtx, fapIdx, quantityCallback); }
+
+    /*!
+     * \brief Calculates the gradient of an arbitrary quantity at any
+     *        flux approximation point on the boundary.
+     *
+     * Boundary gradients are always calculated using the two-point
+     * approximation.
+     *
+     * \param elemCtx The current execution context
+     * \param fapIdx The local index of the flux approximation point
+     *               in the current element's stencil.
+     * \param quantityCallback A callable object returning the value
+     *               of the quantity at an index of a degree of
+     *               freedom
+     */
+    template <class QuantityCallback, class EvalDimVector>
+    void calculateBoundaryGradient(EvalDimVector& quantityGrad,
+                                   const ElementContext& elemCtx,
+                                   unsigned fapIdx,
+                                   const QuantityCallback& quantityCallback) const
+    { ParentType::calculateBoundaryGradient(quantityGrad, elemCtx, fapIdx, quantityCallback); }
+
+#if HAVE_DUNE_LOCALFUNCTIONS
+    static LocalFiniteElementCache& localFiniteElementCache()
+    { return feCache_; }
+#endif
+
+private:
+#if HAVE_DUNE_LOCALFUNCTIONS
+    static LocalFiniteElementCache feCache_;
+
+    const LocalFiniteElement* localFiniteElement_;
+    std::vector<Dune::FieldVector<Scalar, 1>> p1Value_[maxFap];
+    DimVector p1Gradient_[maxFap][maxDof];
+#endif // HAVE_DUNE_LOCALFUNCTIONS
+};
+
+#if HAVE_DUNE_LOCALFUNCTIONS
+template<class TypeTag>
+typename P1FeGradientCalculator<TypeTag>::LocalFiniteElementCache
+P1FeGradientCalculator<TypeTag>::feCache_;
+#endif
+} // namespace Opm
+
+#undef EWOMS_NO_LOCALFUNCTIONS_UNUSED
+
+#endif

opm/models/discretization/vcfv/vcfvbaseoutputmodule.hh

@@ -0,0 +1,85 @@
+// -*- 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::VcfvBaseOutputModule
+ */
+#ifndef EWOMS_VCFV_VTK_BASE_OUTPUT_MODULE_HH
+#define EWOMS_VCFV_VTK_BASE_OUTPUT_MODULE_HH
+
+#include "vcfvproperties.hh"
+
+#include <ewoms/io/baseoutputwriter.hh>
+
+#include <string>
+#include <vector>
+
+namespace Opm {
+/*!
+ * \ingroup VcfvDiscretization
+ *
+ * \brief Implements the discretization specific parts of writing files.
+ */
+template<class TypeTag>
+class VcfvBaseOutputModule
+{
+public:
+    typedef BaseOutputWriter::Scalar Scalar;
+    typedef BaseOutputWriter::Vector Vector;
+    typedef BaseOutputWriter::ScalarBuffer ScalarBuffer;
+    typedef BaseOutputWriter::VectorBuffer VectorBuffer;
+    typedef BaseOutputWriter::TensorBuffer TensorBuffer;
+
+    /*!
+     * \brief Add a buffer where the data is associated with the
+     *        degrees of freedom to the current VTK output file.
+     */
+    static void attachScalarDofData_(BaseOutputWriter& baseWriter,
+                                     ScalarBuffer& buffer,
+                                     const std::string& name)
+    { baseWriter.attachScalarVertexData(buffer, name.c_str()); }
+
+    /*!
+     * \brief Add a buffer where the data is associated with the
+     *        degrees of freedom to the current VTK output file.
+     */
+    static void attachVectorDofData_(BaseOutputWriter& baseWriter,
+                                     VectorBuffer& buffer,
+                                     const std::string& name)
+    { baseWriter.attachVectorVertexData(buffer, name.c_str()); }
+
+
+    /*!
+     * \brief Add a buffer where the data is associated with the
+     *        degrees of freedom to the current VTK output file.
+     */
+    static void attachTensorDofData_(BaseOutputWriter& baseWriter,
+                                     TensorBuffer& buffer,
+                                     const std::string& name)
+    { baseWriter.attachTensorVertexData(buffer, name.c_str()); }
+};
+
+} // namespace Opm
+
+#endif

opm/models/discretization/vcfv/vcfvdiscretization.hh

@@ -0,0 +1,198 @@
+// -*- 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::VcfvDiscretization
+ */
+#ifndef EWOMS_VCFV_DISCRETIZATION_HH
+#define EWOMS_VCFV_DISCRETIZATION_HH
+
+#include <opm/material/densead/Math.hpp>
+
+#include "vcfvproperties.hh"
+#include "vcfvstencil.hh"
+#include "p1fegradientcalculator.hh"
+#include "vcfvgridcommhandlefactory.hh"
+#include "vcfvbaseoutputmodule.hh"
+
+#include <opm/simulators/linalg/vertexborderlistfromgrid.hh>
+#include <opm/models/discretization/common/fvbasediscretization.hh>
+
+#if HAVE_DUNE_FEM
+#include <dune/fem/space/common/functionspace.hh>
+#include <dune/fem/space/lagrange.hh>
+#endif
+
+namespace Opm {
+template <class TypeTag>
+class VcfvDiscretization;
+
+} // namespace Opm
+
+BEGIN_PROPERTIES
+
+//! Set the stencil
+SET_PROP(VcfvDiscretization, Stencil)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::ctype CoordScalar;
+
+public:
+    typedef Opm::VcfvStencil<CoordScalar, GridView> type;
+};
+
+//! Mapper for the degrees of freedoms.
+SET_TYPE_PROP(VcfvDiscretization, DofMapper, typename GET_PROP_TYPE(TypeTag, VertexMapper));
+
+//! The concrete class which manages the spatial discretization
+SET_TYPE_PROP(VcfvDiscretization, Discretization, Opm::VcfvDiscretization<TypeTag>);
+
+//! The base class for the output modules (decides whether to write
+//! element or vertex based fields)
+SET_TYPE_PROP(VcfvDiscretization, DiscBaseOutputModule,
+              Opm::VcfvBaseOutputModule<TypeTag>);
+
+//! Calculates the gradient of any quantity given the index of a flux approximation point
+SET_TYPE_PROP(VcfvDiscretization, GradientCalculator,
+              Opm::P1FeGradientCalculator<TypeTag>);
+
+//! The class to create grid communication handles
+SET_TYPE_PROP(VcfvDiscretization, GridCommHandleFactory,
+              Opm::VcfvGridCommHandleFactory<TypeTag>);
+
+//! Use two-point gradients by default for the vertex centered finite volume scheme.
+SET_BOOL_PROP(VcfvDiscretization, UseP1FiniteElementGradients, false);
+
+#if HAVE_DUNE_FEM
+//! Set the DiscreteFunctionSpace
+SET_PROP(VcfvDiscretization, DiscreteFunctionSpace)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar)   Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, GridPart) GridPart;
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    typedef Dune::Fem::FunctionSpace<typename GridPart::GridType::ctype,
+                                     Scalar,
+                                     GridPart::GridType::dimensionworld,
+                                     numEq> FunctionSpace;
+public:
+    // Lagrange discrete function space with unknowns at the cell vertices
+    typedef Dune::Fem::LagrangeDiscreteFunctionSpace< FunctionSpace, GridPart, 1 > type;
+};
+#endif
+
+//! Set the border list creator for vertices
+SET_PROP(VcfvDiscretization, BorderListCreator)
+{ private:
+    typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+public:
+    typedef Opm::Linear::VertexBorderListFromGrid<GridView, VertexMapper> type;
+};
+
+//! For the vertex centered finite volume method, ghost and overlap elements must _not_
+//! be assembled to avoid accounting twice for the fluxes over the process boundary faces
+//! of the local process' grid partition
+SET_BOOL_PROP(VcfvDiscretization, LinearizeNonLocalElements, false);
+
+
+END_PROPERTIES
+
+namespace Opm {
+
+/*!
+ * \ingroup VcfvDiscretization
+ *
+ * \brief The base class for the vertex centered finite volume discretization scheme.
+ */
+template<class TypeTag>
+class VcfvDiscretization : public FvBaseDiscretization<TypeTag>
+{
+    typedef FvBaseDiscretization<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
+    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+
+    enum { dim = GridView::dimension };
+
+public:
+    VcfvDiscretization(Simulator& simulator)
+        : ParentType(simulator)
+    { }
+
+    /*!
+     * \brief Returns a string of discretization's human-readable name
+     */
+    static std::string discretizationName()
+    { return "vcfv"; }
+
+    /*!
+     * \brief Returns the number of global degrees of freedom (DOFs) due to the grid
+     */
+    size_t numGridDof() const
+    { return static_cast<size_t>(this->gridView_.size(/*codim=*/dim)); }
+
+    /*!
+     * \brief Mapper to convert the Dune entities of the
+     *        discretization's degrees of freedoms are to indices.
+     */
+    const DofMapper& dofMapper() const
+    { return this->vertexMapper(); }
+
+    /*!
+     * \brief Serializes the current state of the model.
+     *
+     * \tparam Restarter The type of the serializer class
+     *
+     * \param res The serializer object
+     */
+    template <class Restarter>
+    void serialize(Restarter& res)
+    { res.template serializeEntities</*codim=*/dim>(asImp_(), this->gridView_); }
+
+    /*!
+     * \brief Deserializes the state of the model.
+     *
+     * \tparam Restarter The type of the serializer class
+     *
+     * \param res The serializer object
+     */
+    template <class Restarter>
+    void deserialize(Restarter& res)
+    {
+        res.template deserializeEntities</*codim=*/dim>(asImp_(), this->gridView_);
+        this->solution(/*timeIdx=*/1) = this->solution(/*timeIdx=*/0);
+    }
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+    const Implementation& asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+} // namespace Opm
+
+#endif

opm/models/discretization/vcfv/vcfvgridcommhandlefactory.hh

@@ -0,0 +1,91 @@
+// -*- 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::VcfvGridCommHandleFactory
+ */
+#ifndef EWOMS_VCFV_GRID_COMM_HANDLE_FACTORY_HH
+#define EWOMS_VCFV_GRID_COMM_HANDLE_FACTORY_HH
+
+#include "vcfvproperties.hh"
+
+#include <opm/models/parallel/gridcommhandles.hh>
+
+namespace Opm {
+/*!
+ * \ingroup VcfvDiscretization
+ *
+ * \brief A class which provides types for DUNE grid handles for
+ *        communication.
+ *
+ * This is required for parallel computations
+ */
+template<class TypeTag>
+class VcfvGridCommHandleFactory
+{
+    typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    static const int dim = GridView::dimension;
+
+public:
+    /*!
+     * \brief Return a handle which computes the minimum of a value
+     *        for each overlapping degree of freedom across all processes.
+     */
+    template <class ValueType, class ArrayType>
+    static std::shared_ptr<GridCommHandleMin<ValueType, ArrayType,  DofMapper, /*commCodim=*/dim> >
+    minHandle(ArrayType& array, const DofMapper& dofMapper)
+    {
+        typedef GridCommHandleMin<ValueType, ArrayType,  DofMapper, /*commCodim=*/dim> Handle;
+        return  std::shared_ptr<Handle>(new Handle(array, dofMapper));
+    }
+
+    /*!
+     * \brief Return a handle which computes the maximum of a value
+     *        for each overlapping degree of freedom across all processes.
+     */
+    template <class ValueType, class ArrayType>
+    static std::shared_ptr<GridCommHandleMax<ValueType, ArrayType,  DofMapper, /*commCodim=*/dim> >
+    maxHandle(ArrayType& array, const DofMapper& dofMapper)
+    {
+        typedef GridCommHandleMax<ValueType, ArrayType,  DofMapper, /*commCodim=*/dim> Handle;
+        return  std::shared_ptr<Handle>(new Handle(array, dofMapper));
+    }
+
+    /*!
+     * \brief Return a handle which computes the sum of all values
+     *        all overlapping degrees of freedom across all processes.
+     */
+    template <class ValueType, class ArrayType>
+    static std::shared_ptr<GridCommHandleSum<ValueType, ArrayType,  DofMapper, /*commCodim=*/dim> >
+    sumHandle(ArrayType& array, const DofMapper& dofMapper)
+    {
+        typedef GridCommHandleSum<ValueType, ArrayType,  DofMapper, /*commCodim=*/dim> Handle;
+        return  std::shared_ptr<Handle>(new Handle(array, dofMapper));
+    }
+};
+} // namespace Opm
+
+#endif

opm/models/discretization/vcfv/vcfvproperties.hh

@@ -0,0 +1,47 @@
+// -*- 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
+ *
+ * \brief Declares the basic properties used by the common infrastructure of
+ *        the vertex-centered finite volume discretization.
+ */
+#ifndef EWOMS_VCFV_PROPERTIES_HH
+#define EWOMS_VCFV_PROPERTIES_HH
+
+#include <opm/models/discretization/common/fvbaseproperties.hh>
+
+#include <opm/models/utils/propertysystem.hh>
+
+BEGIN_PROPERTIES
+
+//! The type tag for models based on the VCFV-scheme
+NEW_TYPE_TAG(VcfvDiscretization, INHERITS_FROM(FvBaseDiscretization));
+
+//! Use P1 finite-elements gradients instead of two-point gradients. Note that setting
+//! this property to true requires the dune-localfunctions module to be available.
+NEW_PROP_TAG(UseP1FiniteElementGradients);
+
+END_PROPERTIES
+
+#endif

tests/models/test_quadrature.cpp

@@ -39,7 +39,7 @@
 #include <dune/common/version.hh>
 #include <dune/geometry/quadraturerules.hh>
 
-#include <ewoms/disc/vcfv/vcfvstencil.hh>
+#include <opm/models/discretization/vcfv/vcfvstencil.hh>
 
 #include <opm/material/common/Unused.hpp>