fvbasediscretization.hh

Go to the documentation of this file.

// -*- 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_FV_BASE_DISCRETIZATION_HH
#define EWOMS_FV_BASE_DISCRETIZATION_HH
 
#include <dune/common/version.hh>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/bvector.hh>
 
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/densead/Math.hpp>
 
#include <opm/models/discretization/common/baseauxiliarymodule.hh>
#include <opm/models/discretization/common/fvbaseadlocallinearizer.hh>
#include <opm/models/discretization/common/fvbaseboundarycontext.hh>
#include <opm/models/discretization/common/fvbaseconstraints.hh>
#include <opm/models/discretization/common/fvbaseconstraintscontext.hh>
#include <opm/models/discretization/common/fvbaseelementcontext.hh>
#include <opm/models/discretization/common/fvbaseextensivequantities.hh>
#include <opm/models/discretization/common/fvbasefdlocallinearizer.hh>
#include <opm/models/discretization/common/fvbasegradientcalculator.hh>
#include <opm/models/discretization/common/fvbaseintensivequantities.hh>
#include <opm/models/discretization/common/fvbaselinearizer.hh>
#include <opm/models/discretization/common/fvbaselocalresidual.hh>
#include <opm/models/discretization/common/fvbasenewtonmethod.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
 
#include <opm/models/io/vtkprimaryvarsmodule.hpp>
 
#include <opm/models/parallel/gridcommhandles.hh>
#include <opm/models/parallel/threadmanager.hpp>
 
#include <opm/models/utils/alignedallocator.hh>
#include <opm/models/utils/simulator.hh>
#include <opm/models/utils/timer.hpp>
#include <opm/models/utils/timerguard.hh>
 
#include <opm/simulators/linalg/linalgparameters.hh>
#include <opm/simulators/linalg/nullborderlistmanager.hh>
 
#include <algorithm>
#include <cstddef>
#include <list>
#include <stdexcept>
#include <sstream>
#include <string>
#include <type_traits>
#include <vector>
 
namespace Opm {
template<class TypeTag>
class FvBaseDiscretizationNoAdapt;
template<class TypeTag>
class FvBaseDiscretization;
 
} // namespace Opm
 
namespace Opm::Properties {
 
template<class TypeTag>
struct Simulator<TypeTag, TTag::FvBaseDiscretization>
{ using type = ::Opm::Simulator<TypeTag>; };
 
template<class TypeTag>
struct VertexMapper<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
 
template<class TypeTag>
struct ElementMapper<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::MultipleCodimMultipleGeomTypeMapper<GetPropType<TypeTag, Properties::GridView>>; };
 
template<class TypeTag>
struct BorderListCreator<TypeTag, TTag::FvBaseDiscretization>
{
    using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
 
public:
    using type = Linear::NullBorderListCreator<GridView, DofMapper>;
};
 
template<class TypeTag>
struct DiscLocalResidual<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseLocalResidual<TypeTag>; };
 
template<class TypeTag>
struct DiscIntensiveQuantities<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseIntensiveQuantities<TypeTag>; };
 
template<class TypeTag>
struct DiscExtensiveQuantities<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseExtensiveQuantities<TypeTag>; };
 
template<class TypeTag>
struct GradientCalculator<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseGradientCalculator<TypeTag>; };
 
template<class TypeTag>
struct EqVector<TypeTag, TTag::FvBaseDiscretization>
{
    using type = Dune::FieldVector<GetPropType<TypeTag, Properties::Scalar>,
                                   getPropValue<TypeTag, Properties::NumEq>()>;
};
 
template<class TypeTag>
struct RateVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = GetPropType<TypeTag, Properties::EqVector>; };
 
template<class TypeTag>
struct BoundaryRateVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = GetPropType<TypeTag, Properties::RateVector>; };
 
template<class TypeTag>
struct Constraints<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseConstraints<TypeTag>; };
 
template<class TypeTag>
struct ElementEqVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
 
template<class TypeTag>
struct GlobalEqVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::EqVector>>; };
 
template<class TypeTag>
struct PrimaryVariables<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBasePrimaryVariables<TypeTag>; };
 
template<class TypeTag>
struct SolutionVector<TypeTag, TTag::FvBaseDiscretization>
{ using type = Dune::BlockVector<GetPropType<TypeTag, Properties::PrimaryVariables>>; };
 
template<class TypeTag>
struct IntensiveQuantities<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseIntensiveQuantities<TypeTag>; };
 
template<class TypeTag>
struct ElementContext<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseElementContext<TypeTag>; };
 
template<class TypeTag>
struct BoundaryContext<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseBoundaryContext<TypeTag>; };
 
template<class TypeTag>
struct ConstraintsContext<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseConstraintsContext<TypeTag>; };
 
template<class TypeTag>
struct ThreadManager<TypeTag, TTag::FvBaseDiscretization>
{ using type = ::Opm::ThreadManager; };
 
template<class TypeTag>
struct UseLinearizationLock<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = true; };
 
template<class TypeTag>
struct Linearizer<TypeTag, TTag::FvBaseDiscretization>
{ using type = FvBaseLinearizer<TypeTag>; };
 
template<class TypeTag>
struct VtkOutputFormat<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr int value = Dune::VTK::ascii; };
 
// disable constraints by default
template<class TypeTag>
struct EnableConstraints<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = false; };
 
template<class TypeTag>
struct TimeDiscHistorySize<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr int value = 2; };
 
template<class TypeTag>
struct ExtensiveStorageTerm<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = false; };
 
// use volumetric residuals is default
template<class TypeTag>
struct UseVolumetricResidual<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = true; };
 
template<class TypeTag>
struct EnableExperiments<TypeTag, TTag::FvBaseDiscretization>
{ static constexpr bool value = true; };
 
template <class TypeTag, class MyTypeTag>
struct BaseDiscretizationType { using type = UndefinedProperty; };
 
#if !HAVE_DUNE_FEM
template<class TypeTag>
struct BaseDiscretizationType<TypeTag,TTag::FvBaseDiscretization>
{ using type = FvBaseDiscretizationNoAdapt<TypeTag>; };
 
template<class TypeTag>
struct DiscreteFunction<TypeTag, TTag::FvBaseDiscretization>
{
    using BaseDiscretization = FvBaseDiscretization<TypeTag>;
    using type = typename BaseDiscretization::BlockVectorWrapper;
};
#endif
 
} // namespace Opm::Properties
 
namespace Opm {
 
template<class TypeTag>
class FvBaseDiscretization
{
    using Implementation = GetPropType<TypeTag, Properties::Model>;
    using Discretization = GetPropType<TypeTag, Properties::Discretization>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Grid = GetPropType<TypeTag, Properties::Grid>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using ElementMapper = GetPropType<TypeTag, Properties::ElementMapper>;
    using VertexMapper = GetPropType<TypeTag, Properties::VertexMapper>;
    using DofMapper = GetPropType<TypeTag, Properties::DofMapper>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
    using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using Linearizer = GetPropType<TypeTag, Properties::Linearizer>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using BoundaryContext = GetPropType<TypeTag, Properties::BoundaryContext>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
    using GradientCalculator = GetPropType<TypeTag, Properties::GradientCalculator>;
    using Stencil = GetPropType<TypeTag, Properties::Stencil>;
    using DiscBaseOutputModule = GetPropType<TypeTag, Properties::DiscBaseOutputModule>;
    using GridCommHandleFactory = GetPropType<TypeTag, Properties::GridCommHandleFactory>;
    using NewtonMethod = GetPropType<TypeTag, Properties::NewtonMethod>;
    using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
 
    using LocalLinearizer = GetPropType<TypeTag, Properties::LocalLinearizer>;
    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
 
    enum {
        numEq = getPropValue<TypeTag, Properties::NumEq>(),
        historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>(),
    };
 
    using IntensiveQuantitiesVector = std::vector<IntensiveQuantities, aligned_allocator<IntensiveQuantities, alignof(IntensiveQuantities)> >;
 
    using Element = typename GridView::template Codim<0>::Entity;
    using ElementIterator = typename GridView::template Codim<0>::Iterator;
 
    using Toolbox = MathToolbox<Evaluation>;
    using VectorBlock = Dune::FieldVector<Evaluation, numEq>;
    using EvalEqVector = Dune::FieldVector<Evaluation, numEq>;
 
    using LocalEvalBlockVector = typename LocalResidual::LocalEvalBlockVector;
 
public:
    class BlockVectorWrapper
    {
    protected:
        SolutionVector blockVector_;
    public:
        BlockVectorWrapper(const std::string&, const size_t size)
            : blockVector_(size)
        {}
 
        BlockVectorWrapper() = default;
 
        static BlockVectorWrapper serializationTestObject()
        {
            BlockVectorWrapper result("dummy", 3);
            result.blockVector_[0] = 1.0;
            result.blockVector_[1] = 2.0;
            result.blockVector_[2] = 3.0;
 
            return result;
        }
 
        SolutionVector& blockVector()
        { return blockVector_; }
        const SolutionVector& blockVector() const
        { return blockVector_; }
 
        bool operator==(const BlockVectorWrapper& wrapper) const
        {
            return std::equal(this->blockVector_.begin(), this->blockVector_.end(),
                              wrapper.blockVector_.begin(), wrapper.blockVector_.end());
        }
 
        template<class Serializer>
        void serializeOp(Serializer& serializer)
        {
            serializer(blockVector_);
        }
    };
 
private:
    using DiscreteFunctionSpace = GetPropType<TypeTag, Properties::DiscreteFunctionSpace>;
    using DiscreteFunction = GetPropType<TypeTag, Properties::DiscreteFunction>;
 
    // copying a discretization object is not a good idea
    FvBaseDiscretization(const FvBaseDiscretization& );
 
public:
    // this constructor required to be explicitly specified because
    // we've defined a constructor above which deletes all implicitly
    // generated constructors in C++.
    FvBaseDiscretization(Simulator& simulator)
        : simulator_(simulator)
        , gridView_(simulator.gridView())
        , elementMapper_(gridView_, Dune::mcmgElementLayout())
        , vertexMapper_(gridView_, Dune::mcmgVertexLayout())
        , newtonMethod_(simulator)
        , localLinearizer_(ThreadManager::maxThreads())
        , linearizer_(new Linearizer())
        , enableGridAdaptation_(Parameters::Get<Parameters::EnableGridAdaptation>() )
        , enableIntensiveQuantityCache_(Parameters::Get<Parameters::EnableIntensiveQuantityCache>())
        , enableStorageCache_(Parameters::Get<Parameters::EnableStorageCache>())
        , enableThermodynamicHints_(Parameters::Get<Parameters::EnableThermodynamicHints>())
    {
        bool isEcfv = std::is_same<Discretization, EcfvDiscretization<TypeTag> >::value;
        if (enableGridAdaptation_ && !isEcfv)
            throw std::invalid_argument("Grid adaptation currently only works for the "
                                        "element-centered finite volume discretization (is: "
                                        +Dune::className<Discretization>()+")");
 
        PrimaryVariables::init();
        size_t numDof = asImp_().numGridDof();
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            if (storeIntensiveQuantities()) {
                intensiveQuantityCache_[timeIdx].resize(numDof);
                intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof, /*value=*/false);
            }
 
            if (enableStorageCache_)
                storageCache_[timeIdx].resize(numDof);
        }
 
        resizeAndResetIntensiveQuantitiesCache_();
        asImp_().registerOutputModules_();
    }
 
    ~FvBaseDiscretization()
    {
        // delete all output modules
        auto modIt = outputModules_.begin();
        const auto& modEndIt = outputModules_.end();
        for (; modIt != modEndIt; ++modIt)
            delete *modIt;
 
        delete linearizer_;
    }
 
    static void registerParameters()
    {
        Linearizer::registerParameters();
        LocalLinearizer::registerParameters();
        LocalResidual::registerParameters();
        GradientCalculator::registerParameters();
        IntensiveQuantities::registerParameters();
        ExtensiveQuantities::registerParameters();
        NewtonMethod::registerParameters();
        Linearizer::registerParameters();
        PrimaryVariables::registerParameters();
        // register runtime parameters of the output modules
        VtkPrimaryVarsModule<TypeTag>::registerParameters();
 
        Parameters::Register<Parameters::EnableGridAdaptation>
            ("Enable adaptive grid refinement/coarsening");
        Parameters::Register<Parameters::EnableVtkOutput>
            ("Global switch for turning on writing VTK files");
        Parameters::Register<Parameters::EnableThermodynamicHints>
            ("Enable thermodynamic hints");
        Parameters::Register<Parameters::EnableIntensiveQuantityCache>
            ("Turn on caching of intensive quantities");
        Parameters::Register<Parameters::EnableStorageCache>
            ("Store previous storage terms and avoid re-calculating them.");
        Parameters::Register<Parameters::OutputDir>
            ("The directory to which result files are written");
    }
 
    void finishInit()
    {
        // initialize the volume of the finite volumes to zero
        size_t numDof = asImp_().numGridDof();
        dofTotalVolume_.resize(numDof);
        std::fill(dofTotalVolume_.begin(), dofTotalVolume_.end(), 0.0);
 
        ElementContext elemCtx(simulator_);
        gridTotalVolume_ = 0.0;
 
        // iterate through the grid and evaluate the initial condition
        for (const auto& elem : elements(gridView_)) {
            const bool isInteriorElement = elem.partitionType() == Dune::InteriorEntity;
            // ignore everything which is not in the interior if the
            // current process' piece of the grid
            if (!isInteriorElement)
                continue;
 
            // deal with the current element
            elemCtx.updateStencil(elem);
            const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
 
            // loop over all element vertices, i.e. sub control volumes
            for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++) {
                // map the local degree of freedom index to the global one
                unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
 
                Scalar dofVolume = stencil.subControlVolume(dofIdx).volume();
                dofTotalVolume_[globalIdx] += dofVolume;
                if (isInteriorElement)
                    gridTotalVolume_ += dofVolume;
            }
        }
 
        // determine which DOFs should be considered to lie fully in the interior of the
        // local process grid partition: those which do not have a non-zero volume
        // before taking the peer processes into account...
        isLocalDof_.resize(numDof);
        for (unsigned dofIdx = 0; dofIdx < numDof; ++dofIdx)
            isLocalDof_[dofIdx] = (dofTotalVolume_[dofIdx] != 0.0);
 
        // add the volumes of the DOFs on the process boundaries
        const auto sumHandle =
            GridCommHandleFactory::template sumHandle<Scalar>(dofTotalVolume_,
                                                              asImp_().dofMapper());
        gridView_.communicate(*sumHandle,
                              Dune::InteriorBorder_All_Interface,
                              Dune::ForwardCommunication);
 
        // sum up the volumes of the grid partitions
        gridTotalVolume_ = gridView_.comm().sum(gridTotalVolume_);
 
        linearizer_->init(simulator_);
        for (unsigned threadId = 0; threadId < ThreadManager::maxThreads(); ++threadId)
            localLinearizer_[threadId].init(simulator_);
 
        resizeAndResetIntensiveQuantitiesCache_();
        if (storeIntensiveQuantities()) {
            // invalidate all cached intensive quantities
            for (unsigned timeIdx = 0; timeIdx < historySize; ++ timeIdx)
                invalidateIntensiveQuantitiesCache(timeIdx);
        }
 
        newtonMethod_.finishInit();
    }
 
    bool enableGridAdaptation() const
    { return enableGridAdaptation_; }
 
    void applyInitialSolution()
    {
        // first set the whole domain to zero
        SolutionVector& uCur = asImp_().solution(/*timeIdx=*/0);
        uCur = Scalar(0.0);
 
        ElementContext elemCtx(simulator_);
 
        // iterate through the grid and evaluate the initial condition
        for (const auto& elem : elements(gridView_)) {
            // ignore everything which is not in the interior if the
            // current process' piece of the grid
            if (elem.partitionType() != Dune::InteriorEntity)
                continue;
 
            // deal with the current element
            elemCtx.updateStencil(elem);
 
            // loop over all element vertices, i.e. sub control volumes
            for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); dofIdx++)
            {
                // map the local degree of freedom index to the global one
                unsigned globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
 
                // let the problem do the dirty work of nailing down
                // the initial solution.
                simulator_.problem().initial(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
                asImp_().supplementInitialSolution_(uCur[globalIdx], elemCtx, dofIdx, /*timeIdx=*/0);
                uCur[globalIdx].checkDefined();
            }
        }
 
        // synchronize the ghost DOFs (if necessary)
        asImp_().syncOverlap();
 
        // also set the solutions of the "previous" time steps to the initial solution.
        for (unsigned timeIdx = 1; timeIdx < historySize; ++timeIdx)
            solution(timeIdx) = solution(/*timeIdx=*/0);
 
        simulator_.problem().initialSolutionApplied();
 
#ifndef NDEBUG
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx)  {
            const auto& sol = solution(timeIdx);
            for (unsigned dofIdx = 0; dofIdx < sol.size(); ++dofIdx)
                sol[dofIdx].checkDefined();
        }
#endif // NDEBUG
    }
 
    void prefetch(const Element&) const
    {
        // do nothing by default
    }
 
    NewtonMethod& newtonMethod()
    { return newtonMethod_; }
 
    const NewtonMethod& newtonMethod() const
    { return newtonMethod_; }
 
    const IntensiveQuantities* thermodynamicHint(unsigned globalIdx, unsigned timeIdx) const
    {
        if (!enableThermodynamicHints_)
            return 0;
 
        // the intensive quantities cache doubles as thermodynamic hint
        return cachedIntensiveQuantities(globalIdx, timeIdx);
    }
 
    const IntensiveQuantities* cachedIntensiveQuantities(unsigned globalIdx, unsigned timeIdx) const
    {
        if (!enableIntensiveQuantityCache_ || !intensiveQuantityCacheUpToDate_[timeIdx][globalIdx]) {
            return nullptr;
        }
 
        // With the storage cache enabled, usually only the
        // intensive quantities for the most recent time step are
        // cached. However, this may be false for some Problem
        // variants, so we should check if the cache exists for
        // the timeIdx in question.
        if (timeIdx > 0 && enableStorageCache_ && intensiveQuantityCache_[timeIdx].empty()) {
            return nullptr;
        }
 
        return &intensiveQuantityCache_[timeIdx][globalIdx];
    }
 
    void updateCachedIntensiveQuantities(const IntensiveQuantities& intQuants,
                                         unsigned globalIdx,
                                         unsigned timeIdx) const
    {
        if (!storeIntensiveQuantities())
            return;
 
        intensiveQuantityCache_[timeIdx][globalIdx] = intQuants;
        intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = 1;
    }
 
    void setIntensiveQuantitiesCacheEntryValidity(unsigned globalIdx,
                                                  unsigned timeIdx,
                                                  bool newValue) const
    {
        if (!storeIntensiveQuantities())
            return;
 
        intensiveQuantityCacheUpToDate_[timeIdx][globalIdx] = newValue ? 1 : 0;
    }
 
    void invalidateIntensiveQuantitiesCache(unsigned timeIdx) const
    {
        if (storeIntensiveQuantities()) {
            std::fill(intensiveQuantityCacheUpToDate_[timeIdx].begin(),
                      intensiveQuantityCacheUpToDate_[timeIdx].end(),
                      /*value=*/0);
        }
    }
 
    void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx) const
    {
        invalidateIntensiveQuantitiesCache(timeIdx);
 
        // loop over all elements...
        ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
#ifdef _OPENMP
#pragma omp parallel
#endif
        {
            ElementContext elemCtx(simulator_);
            ElementIterator elemIt = threadedElemIt.beginParallel();
            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
                const Element& elem = *elemIt;
                elemCtx.updatePrimaryStencil(elem);
                elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
            }
        }
    }
 
    template <class GridViewType>
    void invalidateAndUpdateIntensiveQuantities(unsigned timeIdx, const GridViewType& gridView) const
    {
        // loop over all elements...
        ThreadedEntityIterator<GridViewType, /*codim=*/0> threadedElemIt(gridView);
#ifdef _OPENMP
#pragma omp parallel
#endif
        {
 
            ElementContext elemCtx(simulator_);
            auto elemIt = threadedElemIt.beginParallel();
            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
                if (elemIt->partitionType() != Dune::InteriorEntity) {
                    continue;
                }
                const Element& elem = *elemIt;
                elemCtx.updatePrimaryStencil(elem);
                // Mark cache for this element as invalid.
                const std::size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
                for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
                    const unsigned globalIndex = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
                    setIntensiveQuantitiesCacheEntryValidity(globalIndex, timeIdx, false);
                }
                // Update for this element.
                elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
            }
        }
    }
 
    void shiftIntensiveQuantityCache(unsigned numSlots = 1)
    {
        if (!storeIntensiveQuantities())
            return;
 
        if (enableStorageCache() && simulator_.problem().recycleFirstIterationStorage()) {
            // If the storage term is cached, the intensive quantities of the previous
            // time steps do not need to be accessed, and we can thus spare ourselves to
            // copy the objects for the intensive quantities.
            // However, if the storage term at the start of the timestep cannot be deduced
            // from the primary variables, we must calculate it from the old intensive
            // quantities, and need to shift them.
            return;
        }
 
        assert(numSlots > 0);
 
        for (unsigned timeIdx = 0; timeIdx < historySize - numSlots; ++ timeIdx) {
            intensiveQuantityCache_[timeIdx + numSlots] = intensiveQuantityCache_[timeIdx];
            intensiveQuantityCacheUpToDate_[timeIdx + numSlots] = intensiveQuantityCacheUpToDate_[timeIdx];
        }
 
        // the cache for the most recent time indices do not need to be invalidated
        // because the solution for them did not change (TODO: that assumes that there is
        // no post-processing of the solution after a time step! fix it?)
    }
 
    bool enableStorageCache() const
    { return enableStorageCache_; }
 
    void setEnableStorageCache(bool enableStorageCache)
    { enableStorageCache_= enableStorageCache; }
 
    const EqVector& cachedStorage(unsigned globalIdx, unsigned timeIdx) const
    {
        assert(enableStorageCache_);
        return storageCache_[timeIdx][globalIdx];
    }
 
    void updateCachedStorage(unsigned globalIdx, unsigned timeIdx, const EqVector& value) const
    {
        assert(enableStorageCache_);
        storageCache_[timeIdx][globalIdx] = value;
    }
 
    Scalar globalResidual(GlobalEqVector& dest,
                          const SolutionVector& u) const
    {
        SolutionVector tmp(asImp_().solution(/*timeIdx=*/0));
        mutableSolution(/*timeIdx=*/0) = u;
        Scalar res = asImp_().globalResidual(dest);
        mutableSolution(/*timeIdx=*/0) = tmp;
        return res;
    }
 
    Scalar globalResidual(GlobalEqVector& dest) const
    {
        dest = 0;
 
        std::mutex mutex;
        ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView_);
#ifdef _OPENMP
#pragma omp parallel
#endif
        {
            // Attention: the variables below are thread specific and thus cannot be
            // moved in front of the #pragma!
            unsigned threadId = ThreadManager::threadId();
            ElementContext elemCtx(simulator_);
            ElementIterator elemIt = threadedElemIt.beginParallel();
            LocalEvalBlockVector residual, storageTerm;
 
            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
                const Element& elem = *elemIt;
                if (elem.partitionType() != Dune::InteriorEntity)
                    continue;
 
                elemCtx.updateAll(elem);
                residual.resize(elemCtx.numDof(/*timeIdx=*/0));
                storageTerm.resize(elemCtx.numPrimaryDof(/*timeIdx=*/0));
                asImp_().localResidual(threadId).eval(residual, elemCtx);
 
                size_t numPrimaryDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
                mutex.lock();
                for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx) {
                    unsigned globalI = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
                    for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
                        dest[globalI][eqIdx] += Toolbox::value(residual[dofIdx][eqIdx]);
                }
                mutex.unlock();
            }
        }
 
        // add up the residuals on the process borders
        const auto sumHandle =
            GridCommHandleFactory::template sumHandle<EqVector>(dest, asImp_().dofMapper());
        gridView_.communicate(*sumHandle,
                              Dune::InteriorBorder_InteriorBorder_Interface,
                              Dune::ForwardCommunication);
 
        // calculate the square norm of the residual. this is not
        // entirely correct, since the residual for the finite volumes
        // which are on the boundary are counted once for every
        // process. As often in life: shit happens (, we don't care)...
        Scalar result2 = dest.two_norm2();
        result2 = asImp_().gridView().comm().sum(result2);
 
        return std::sqrt(result2);
    }
 
    void globalStorage(EqVector& storage, unsigned timeIdx = 0) const
    {
        storage = 0;
 
        std::mutex mutex;
        ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
#ifdef _OPENMP
#pragma omp parallel
#endif
        {
            // Attention: the variables below are thread specific and thus cannot be
            // moved in front of the #pragma!
            unsigned threadId = ThreadManager::threadId();
            ElementContext elemCtx(simulator_);
            ElementIterator elemIt = threadedElemIt.beginParallel();
            LocalEvalBlockVector elemStorage;
 
            // in this method, we need to disable the storage cache because we want to
            // evaluate the storage term for other time indices than the most recent one
            elemCtx.setEnableStorageCache(false);
 
            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
                const Element& elem = *elemIt;
                if (elem.partitionType() != Dune::InteriorEntity)
                    continue; // ignore ghost and overlap elements
 
                elemCtx.updateStencil(elem);
                elemCtx.updatePrimaryIntensiveQuantities(timeIdx);
 
                size_t numPrimaryDof = elemCtx.numPrimaryDof(timeIdx);
                elemStorage.resize(numPrimaryDof);
 
                localResidual(threadId).evalStorage(elemStorage, elemCtx, timeIdx);
 
                mutex.lock();
                for (unsigned dofIdx = 0; dofIdx < numPrimaryDof; ++dofIdx)
                    for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
                        storage[eqIdx] += Toolbox::value(elemStorage[dofIdx][eqIdx]);
                mutex.unlock();
            }
        }
 
        storage = gridView_.comm().sum(storage);
    }
 
    void checkConservativeness([[maybe_unused]] Scalar tolerance = -1,
                               [[maybe_unused]] bool verbose=false) const
    {
#ifndef NDEBUG
        Scalar totalBoundaryArea(0.0);
        Scalar totalVolume(0.0);
        EvalEqVector totalRate(0.0);
 
        // take the newton tolerance times the total volume of the grid if we're not
        // given an explicit tolerance...
        if (tolerance <= 0) {
            tolerance =
                simulator_.model().newtonMethod().tolerance()
                * simulator_.model().gridTotalVolume()
                * 1000;
        }
 
        // we assume the implicit Euler time discretization for now...
        assert(historySize == 2);
 
        EqVector storageBeginTimeStep(0.0);
        globalStorage(storageBeginTimeStep, /*timeIdx=*/1);
 
        EqVector storageEndTimeStep(0.0);
        globalStorage(storageEndTimeStep, /*timeIdx=*/0);
 
        // calculate the rate at the boundary and the source rate
        ElementContext elemCtx(simulator_);
        elemCtx.setEnableStorageCache(false);
        auto eIt = simulator_.gridView().template begin</*codim=*/0>();
        const auto& elemEndIt = simulator_.gridView().template end</*codim=*/0>();
        for (; eIt != elemEndIt; ++eIt) {
            if (eIt->partitionType() != Dune::InteriorEntity)
                continue; // ignore ghost and overlap elements
 
            elemCtx.updateAll(*eIt);
 
            // handle the boundary terms
            if (elemCtx.onBoundary()) {
                BoundaryContext boundaryCtx(elemCtx);
 
                for (unsigned faceIdx = 0; faceIdx < boundaryCtx.numBoundaryFaces(/*timeIdx=*/0); ++faceIdx) {
                    BoundaryRateVector values;
                    simulator_.problem().boundary(values,
                                                         boundaryCtx,
                                                         faceIdx,
                                                         /*timeIdx=*/0);
                    Valgrind::CheckDefined(values);
 
                    unsigned dofIdx = boundaryCtx.interiorScvIndex(faceIdx, /*timeIdx=*/0);
                    const auto& insideIntQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
 
                    Scalar bfArea =
                        boundaryCtx.boundarySegmentArea(faceIdx, /*timeIdx=*/0)
                        * insideIntQuants.extrusionFactor();
 
                    for (unsigned i = 0; i < values.size(); ++i)
                        values[i] *= bfArea;
 
                    totalBoundaryArea += bfArea;
                    for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
                        totalRate[eqIdx] += values[eqIdx];
                }
            }
 
            // deal with the source terms
            for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++ dofIdx) {
                RateVector values;
                simulator_.problem().source(values,
                                            elemCtx,
                                            dofIdx,
                                            /*timeIdx=*/0);
                Valgrind::CheckDefined(values);
 
                const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
                Scalar dofVolume =
                    elemCtx.dofVolume(dofIdx, /*timeIdx=*/0)
                    * intQuants.extrusionFactor();
                for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx)
                    totalRate[eqIdx] += -dofVolume*Toolbox::value(values[eqIdx]);
                totalVolume += dofVolume;
            }
        }
 
        // summarize everything over all processes
        const auto& comm = simulator_.gridView().comm();
        totalRate = comm.sum(totalRate);
        totalBoundaryArea = comm.sum(totalBoundaryArea);
        totalVolume = comm.sum(totalVolume);
 
        if (comm.rank() == 0) {
            EqVector storageRate = storageBeginTimeStep;
            storageRate -= storageEndTimeStep;
            storageRate /= simulator_.timeStepSize();
            if (verbose) {
                std::cout << "storage at beginning of time step: " << storageBeginTimeStep << "\n";
                std::cout << "storage at end of time step: " << storageEndTimeStep << "\n";
                std::cout << "rate based on storage terms: " << storageRate << "\n";
                std::cout << "rate based on source and boundary terms: " << totalRate << "\n";
                std::cout << "difference in rates: ";
                for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx)
                    std::cout << (storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) << " ";
                std::cout << "\n";
            }
            for (unsigned eqIdx = 0; eqIdx < EqVector::dimension; ++eqIdx) {
                Scalar eps =
                    (std::abs(storageRate[eqIdx]) + Toolbox::value(totalRate[eqIdx]))*tolerance;
                eps = std::max(tolerance, eps);
                assert(std::abs(storageRate[eqIdx] - Toolbox::value(totalRate[eqIdx])) <= eps);
            }
        }
#endif // NDEBUG
    }
 
    Scalar dofTotalVolume(unsigned globalIdx) const
    { return dofTotalVolume_[globalIdx]; }
 
    bool isLocalDof(unsigned globalIdx) const
    { return isLocalDof_[globalIdx]; }
 
    Scalar gridTotalVolume() const
    { return gridTotalVolume_; }
 
    const SolutionVector& solution(unsigned timeIdx) const
    { return solution_[timeIdx]->blockVector(); }
 
    SolutionVector& solution(unsigned timeIdx)
    { return solution_[timeIdx]->blockVector(); }
 
  protected:
    SolutionVector& mutableSolution(unsigned timeIdx) const
    { return solution_[timeIdx]->blockVector(); }
 
  public:
    const Linearizer& linearizer() const
    { return *linearizer_; }
 
    Linearizer& linearizer()
    { return *linearizer_; }
 
    const LocalLinearizer& localLinearizer(unsigned openMpThreadId) const
    { return localLinearizer_[openMpThreadId]; }
    LocalLinearizer& localLinearizer(unsigned openMpThreadId)
    { return localLinearizer_[openMpThreadId]; }
 
    const LocalResidual& localResidual(unsigned openMpThreadId) const
    { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
    LocalResidual& localResidual(unsigned openMpThreadId)
    { return asImp_().localLinearizer(openMpThreadId).localResidual(); }
 
    Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
    {
        Scalar absPv = std::abs(asImp_().solution(/*timeIdx=*/1)[globalDofIdx][pvIdx]);
        return 1.0/std::max(absPv, 1.0);
    }
 
    Scalar eqWeight(unsigned, unsigned) const
    { return 1.0; }
 
    Scalar relativeDofError(unsigned vertexIdx,
                            const PrimaryVariables& pv1,
                            const PrimaryVariables& pv2) const
    {
        Scalar result = 0.0;
        for (unsigned j = 0; j < numEq; ++j) {
            Scalar weight = asImp_().primaryVarWeight(vertexIdx, j);
            Scalar eqErr = std::abs((pv1[j] - pv2[j])*weight);
            //Scalar eqErr = std::abs(pv1[j] - pv2[j]);
            //eqErr *= std::max<Scalar>(1.0, std::abs(pv1[j] + pv2[j])/2);
 
            result = std::max(result, eqErr);
        }
        return result;
    }
 
    bool update()
    {
        TimerGuard prePostProcessGuard(prePostProcessTimer_);
 
#ifndef NDEBUG
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            // Make sure that the primary variables are defined. Note that because of padding
            // bytes, we can't just simply ask valgrind to check the whole solution vectors
            // for definedness...
            for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
                asImp_().solution(timeIdx)[i].checkDefined();
            }
        }
#endif // NDEBUG
 
        // make sure all timers are prestine
        prePostProcessTimer_.halt();
        linearizeTimer_.halt();
        solveTimer_.halt();
        updateTimer_.halt();
 
        prePostProcessTimer_.start();
        asImp_().updateBegin();
        prePostProcessTimer_.stop();
 
        bool converged = false;
 
        try {
            converged = newtonMethod_.apply();
        }
        catch(...) {
            prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
            linearizeTimer_ += newtonMethod_.linearizeTimer();
            solveTimer_ += newtonMethod_.solveTimer();
            updateTimer_ += newtonMethod_.updateTimer();
 
            throw;
        }
 
#ifndef NDEBUG
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            // Make sure that the primary variables are defined. Note that because of padding
            // bytes, we can't just simply ask valgrind to check the whole solution vectors
            // for definedness...
            for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
                asImp_().solution(timeIdx)[i].checkDefined();
            }
        }
#endif // NDEBUG
 
        prePostProcessTimer_ += newtonMethod_.prePostProcessTimer();
        linearizeTimer_ += newtonMethod_.linearizeTimer();
        solveTimer_ += newtonMethod_.solveTimer();
        updateTimer_ += newtonMethod_.updateTimer();
 
        prePostProcessTimer_.start();
        if (converged)
            asImp_().updateSuccessful();
        else
            asImp_().updateFailed();
        prePostProcessTimer_.stop();
 
#ifndef NDEBUG
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            // Make sure that the primary variables are defined. Note that because of padding
            // bytes, we can't just simply ask valgrind to check the whole solution vectors
            // for definedness...
            for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i) {
                asImp_().solution(timeIdx)[i].checkDefined();
            }
        }
#endif // NDEBUG
 
        return converged;
    }
 
    void syncOverlap()
    { }
 
    void updateBegin()
    { }
 
    void updateSuccessful()
    { }
 
    void adaptGrid()
    {
        throw std::invalid_argument("Grid adaptation need to be implemented for "
                                    "specific settings of grid and function spaces");
    }
 
    void updateFailed()
    {
        // Reset the current solution to the one of the
        // previous time step so that we can start the next
        // update at a physically meaningful solution.
        solution(/*timeIdx=*/0) = solution(/*timeIdx=*/1);
        invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
 
#ifndef NDEBUG
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            // Make sure that the primary variables are defined. Note that because of padding
            // bytes, we can't just simply ask valgrind to check the whole solution vectors
            // for definedness...
            for (size_t i = 0; i < asImp_().solution(/*timeIdx=*/0).size(); ++i)
                asImp_().solution(timeIdx)[i].checkDefined();
        }
#endif // NDEBUG
    }
 
    void advanceTimeLevel()
    {
        // at this point we can adapt the grid
        if (this->enableGridAdaptation_) {
            asImp_().adaptGrid();
        }
 
        // make the current solution the previous one.
        solution(/*timeIdx=*/1) = solution(/*timeIdx=*/0);
 
        // shift the intensive quantities cache by one position in the
        // history
        asImp_().shiftIntensiveQuantityCache(/*numSlots=*/1);
    }
 
    template <class Restarter>
    void serialize(Restarter&)
    {
        throw std::runtime_error("Not implemented: The discretization chosen for this problem "
                                 "does not support restart files. (serialize() method unimplemented)");
    }
 
    template <class Restarter>
    void deserialize(Restarter&)
    {
        throw std::runtime_error("Not implemented: The discretization chosen for this problem "
                                 "does not support restart files. (deserialize() method unimplemented)");
    }
 
    template <class DofEntity>
    void serializeEntity(std::ostream& outstream,
                         const DofEntity& dof)
    {
        unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
 
        // write phase state
        if (!outstream.good()) {
            throw std::runtime_error("Could not serialize degree of freedom "
                                     +std::to_string(dofIdx));
        }
 
        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            outstream << solution(/*timeIdx=*/0)[dofIdx][eqIdx] << " ";
        }
    }
 
    template <class DofEntity>
    void deserializeEntity(std::istream& instream,
                           const DofEntity& dof)
    {
        unsigned dofIdx = static_cast<unsigned>(asImp_().dofMapper().index(dof));
 
        for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
            if (!instream.good())
                throw std::runtime_error("Could not deserialize degree of freedom "
                                         +std::to_string(dofIdx));
            instream >> solution(/*timeIdx=*/0)[dofIdx][eqIdx];
        }
    }
 
    size_t numGridDof() const
    { throw std::logic_error("The discretization class must implement the numGridDof() method!"); }
 
    size_t numAuxiliaryDof() const
    {
        size_t result = 0;
        auto auxModIt = auxEqModules_.begin();
        const auto& auxModEndIt = auxEqModules_.end();
        for (; auxModIt != auxModEndIt; ++auxModIt)
            result += (*auxModIt)->numDofs();
 
        return result;
    }
 
    size_t numTotalDof() const
    { return asImp_().numGridDof() + numAuxiliaryDof(); }
 
    const DofMapper& dofMapper() const
    { throw std::logic_error("The discretization class must implement the dofMapper() method!"); }
 
    const VertexMapper& vertexMapper() const
    { return vertexMapper_; }
 
    const ElementMapper& elementMapper() const
    { return elementMapper_; }
 
    void resetLinearizer ()
    {
        delete linearizer_;
        linearizer_ = new Linearizer;
        linearizer_->init(simulator_);
    }
 
    static std::string discretizationName()
    { return ""; }
 
    std::string primaryVarName(unsigned pvIdx) const
    {
        std::ostringstream oss;
        oss << "primary variable_" << pvIdx;
        return oss.str();
    }
 
    std::string eqName(unsigned eqIdx) const
    {
        std::ostringstream oss;
        oss << "equation_" << eqIdx;
        return oss.str();
    }
 
    void updatePVWeights(const ElementContext&) const
    { }
 
    void addOutputModule(BaseOutputModule<TypeTag>* newModule)
    { outputModules_.push_back(newModule); }
 
    template <class VtkMultiWriter>
    void addConvergenceVtkFields(VtkMultiWriter& writer,
                                 const SolutionVector& u,
                                 const GlobalEqVector& deltaU) const
    {
        using ScalarBuffer = std::vector<double>;
 
        GlobalEqVector globalResid(u.size());
        asImp_().globalResidual(globalResid, u);
 
        // create the required scalar fields
        size_t numGridDof = asImp_().numGridDof();
 
        // global defect of the two auxiliary equations
        ScalarBuffer* def[numEq];
        ScalarBuffer* delta[numEq];
        ScalarBuffer* priVars[numEq];
        ScalarBuffer* priVarWeight[numEq];
        ScalarBuffer* relError = writer.allocateManagedScalarBuffer(numGridDof);
        ScalarBuffer* normalizedRelError = writer.allocateManagedScalarBuffer(numGridDof);
        for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
            priVars[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
            priVarWeight[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
            delta[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
            def[pvIdx] = writer.allocateManagedScalarBuffer(numGridDof);
        }
 
        Scalar minRelErr = 1e30;
        Scalar maxRelErr = -1e30;
        for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx) {
            for (unsigned pvIdx = 0; pvIdx < numEq; ++pvIdx) {
                (*priVars[pvIdx])[globalIdx] = u[globalIdx][pvIdx];
                (*priVarWeight[pvIdx])[globalIdx] = asImp_().primaryVarWeight(globalIdx, pvIdx);
                (*delta[pvIdx])[globalIdx] = - deltaU[globalIdx][pvIdx];
                (*def[pvIdx])[globalIdx] = globalResid[globalIdx][pvIdx];
            }
 
            PrimaryVariables uOld(u[globalIdx]);
            PrimaryVariables uNew(uOld);
            uNew -= deltaU[globalIdx];
 
            Scalar err = asImp_().relativeDofError(globalIdx, uOld, uNew);
            (*relError)[globalIdx] = err;
            (*normalizedRelError)[globalIdx] = err;
            minRelErr = std::min(err, minRelErr);
            maxRelErr = std::max(err, maxRelErr);
        }
 
        // do the normalization of the relative error
        Scalar alpha = std::max(Scalar{1e-20},
                                std::max(std::abs(maxRelErr),
                                         std::abs(minRelErr)));
        for (unsigned globalIdx = 0; globalIdx < numGridDof; ++ globalIdx)
            (*normalizedRelError)[globalIdx] /= alpha;
 
        DiscBaseOutputModule::attachScalarDofData_(writer, *relError, "relative error");
        DiscBaseOutputModule::attachScalarDofData_(writer, *normalizedRelError, "normalized relative error");
 
        for (unsigned i = 0; i < numEq; ++i) {
            std::ostringstream oss;
            oss.str(""); oss << "priVar_" << asImp_().primaryVarName(i);
            DiscBaseOutputModule::attachScalarDofData_(writer,
                                                       *priVars[i],
                                                       oss.str());
 
            oss.str(""); oss << "delta_" << asImp_().primaryVarName(i);
            DiscBaseOutputModule::attachScalarDofData_(writer,
                                                       *delta[i],
                                                       oss.str());
 
            oss.str(""); oss << "weight_" << asImp_().primaryVarName(i);
            DiscBaseOutputModule::attachScalarDofData_(writer,
                                                       *priVarWeight[i],
                                                       oss.str());
 
            oss.str(""); oss << "defect_" << asImp_().eqName(i);
            DiscBaseOutputModule::attachScalarDofData_(writer,
                                                       *def[i],
                                                       oss.str());
        }
 
        asImp_().prepareOutputFields();
        asImp_().appendOutputFields(writer);
    }
 
    void prepareOutputFields() const
    {
        bool needFullContextUpdate = false;
        auto modIt = outputModules_.begin();
        const auto& modEndIt = outputModules_.end();
        for (; modIt != modEndIt; ++modIt) {
            (*modIt)->allocBuffers();
            needFullContextUpdate = needFullContextUpdate || (*modIt)->needExtensiveQuantities();
        }
 
        // iterate over grid
        ThreadedEntityIterator<GridView, /*codim=*/0> threadedElemIt(gridView());
#ifdef _OPENMP
#pragma omp parallel
#endif
        {
            ElementContext elemCtx(simulator_);
            ElementIterator elemIt = threadedElemIt.beginParallel();
            for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) {
                const auto& elem = *elemIt;
                if (elem.partitionType() != Dune::InteriorEntity)
                    // ignore non-interior entities
                    continue;
 
                if (needFullContextUpdate)
                    elemCtx.updateAll(elem);
                else {
                    elemCtx.updatePrimaryStencil(elem);
                    elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
                }
 
                // we cannot reuse the "modIt" variable here because the code here might
                // be threaded and "modIt" is is the same for all threads, i.e., if a
                // given thread modifies it, the changes affect all threads.
                auto modIt2 = outputModules_.begin();
                for (; modIt2 != modEndIt; ++modIt2)
                    (*modIt2)->processElement(elemCtx);
            }
        }
    }
 
    void appendOutputFields(BaseOutputWriter& writer) const
    {
        auto modIt = outputModules_.begin();
        const auto& modEndIt = outputModules_.end();
        for (; modIt != modEndIt; ++modIt)
            (*modIt)->commitBuffers(writer);
    }
 
    const GridView& gridView() const
    { return gridView_; }
 
    void addAuxiliaryModule(BaseAuxiliaryModule<TypeTag>* auxMod)
    {
        auxMod->setDofOffset(numTotalDof());
        auxEqModules_.push_back(auxMod);
 
        // resize the solutions
        if (enableGridAdaptation_
            && !std::is_same<DiscreteFunction, BlockVectorWrapper>::value)
        {
            throw std::invalid_argument("Problems which require auxiliary modules cannot be used in"
                                      " conjunction with dune-fem");
        }
 
        size_t numDof = numTotalDof();
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx)
            solution(timeIdx).resize(numDof);
 
        auxMod->applyInitial();
    }
 
    void clearAuxiliaryModules()
    {
        auxEqModules_.clear();
        linearizer_->eraseMatrix();
        newtonMethod_.eraseMatrix();
    }
 
    size_t numAuxiliaryModules() const
    { return auxEqModules_.size(); }
 
    BaseAuxiliaryModule<TypeTag>* auxiliaryModule(unsigned auxEqModIdx)
    { return auxEqModules_[auxEqModIdx]; }
 
    const BaseAuxiliaryModule<TypeTag>* auxiliaryModule(unsigned auxEqModIdx) const
    { return auxEqModules_[auxEqModIdx]; }
 
    bool storeIntensiveQuantities() const
    { return enableIntensiveQuantityCache_ || enableThermodynamicHints_; }
 
    const Timer& prePostProcessTimer() const
    { return prePostProcessTimer_; }
 
    const Timer& linearizeTimer() const
    { return linearizeTimer_; }
 
    const Timer& solveTimer() const
    { return solveTimer_; }
 
    const Timer& updateTimer() const
    { return updateTimer_; }
 
    template<class Serializer>
    void serializeOp(Serializer& serializer)
    {
        using BaseDiscretization = GetPropType<TypeTag, Properties::BaseDiscretizationType>;
        using Helper = typename BaseDiscretization::template SerializeHelper<Serializer>;
        Helper::serializeOp(serializer, solution_);
    }
 
    bool operator==(const FvBaseDiscretization& rhs) const
    {
        return std::equal(this->solution_.begin(), this->solution_.end(),
                          rhs.solution_.begin(), rhs.solution_.end(),
                          [](const auto& x, const auto& y)
                          {
                              return *x == *y;
                          });
    }
 
protected:
    void resizeAndResetIntensiveQuantitiesCache_()
    {
        // allocate the storage cache
        if (enableStorageCache()) {
            size_t numDof = asImp_().numGridDof();
            for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
                storageCache_[timeIdx].resize(numDof);
            }
        }
 
        // allocate the intensive quantities cache
        if (storeIntensiveQuantities()) {
            size_t numDof = asImp_().numGridDof();
            for(unsigned timeIdx=0; timeIdx<historySize; ++timeIdx) {
                intensiveQuantityCache_[timeIdx].resize(numDof);
                intensiveQuantityCacheUpToDate_[timeIdx].resize(numDof);
                invalidateIntensiveQuantitiesCache(timeIdx);
            }
        }
    }
    template <class Context>
    void supplementInitialSolution_(PrimaryVariables&,
                                    const Context&,
                                    unsigned,
                                    unsigned)
    { }
 
    void registerOutputModules_()
    {
        // add the output modules available on all model
        auto *mod = new VtkPrimaryVarsModule<TypeTag>(simulator_);
        this->outputModules_.push_back(mod);
    }
 
    LocalResidual& localResidual_()
    { return localLinearizer_.localResidual(); }
 
    bool verbose_() const
    { return gridView_.comm().rank() == 0; }
 
    Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }
    const Implementation& asImp_() const
    { return *static_cast<const Implementation*>(this); }
 
    // the problem we want to solve. defines the constitutive
    // relations, matxerial laws, etc.
    Simulator& simulator_;
 
    // the representation of the spatial domain of the problem
    GridView gridView_;
 
    // the mappers for element and vertex entities to global indices
    ElementMapper elementMapper_;
    VertexMapper vertexMapper_;
 
    // a vector with all auxiliary equations to be considered
    std::vector<BaseAuxiliaryModule<TypeTag>*> auxEqModules_;
 
    NewtonMethod newtonMethod_;
 
    Timer prePostProcessTimer_;
    Timer linearizeTimer_;
    Timer solveTimer_;
    Timer updateTimer_;
 
    // calculates the local jacobian matrix for a given element
    std::vector<LocalLinearizer> localLinearizer_;
    // Linearizes the problem at the current time step using the
    // local jacobian
    Linearizer *linearizer_;
 
    // cur is the current iterative solution, prev the converged
    // solution of the previous time step
    mutable IntensiveQuantitiesVector intensiveQuantityCache_[historySize];
    // while these are logically bools, concurrent writes to vector<bool> are not thread safe.
    mutable std::vector<unsigned char> intensiveQuantityCacheUpToDate_[historySize];
 
    mutable std::array< std::unique_ptr< DiscreteFunction >, historySize > solution_;
 
    std::list<BaseOutputModule<TypeTag>*> outputModules_;
 
    Scalar gridTotalVolume_;
    std::vector<Scalar> dofTotalVolume_;
    std::vector<bool> isLocalDof_;
 
    mutable GlobalEqVector storageCache_[historySize];
 
    bool enableGridAdaptation_;
    bool enableIntensiveQuantityCache_;
    bool enableStorageCache_;
    bool enableThermodynamicHints_;
};
 
template<class TypeTag>
class FvBaseDiscretizationNoAdapt : public FvBaseDiscretization<TypeTag>
{
    using ParentType = FvBaseDiscretization<TypeTag>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using DiscreteFunction = GetPropType<TypeTag, Properties::DiscreteFunction>;
 
    static constexpr unsigned historySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>();
 
public:
    template<class Serializer>
    struct SerializeHelper {
        template<class SolutionType>
        static void serializeOp(Serializer& serializer,
                                SolutionType& solution)
        {
            for (auto& sol : solution) {
                serializer(*sol);
            }
        }
    };
 
    FvBaseDiscretizationNoAdapt(Simulator& simulator)
        : ParentType(simulator)
    {
        if (this->enableGridAdaptation_) {
            throw std::invalid_argument("Grid adaptation need to use"
                                        " BaseDiscretization = FvBaseDiscretizationFemAdapt"
                                        " which currently requires the presence of the"
                                        " dune-fem module");
        }
        size_t numDof = this->asImp_().numGridDof();
        for (unsigned timeIdx = 0; timeIdx < historySize; ++timeIdx) {
            this->solution_[timeIdx] = std::make_unique<DiscreteFunction>("solution", numDof);
        }
    }
};
 
} // namespace Opm
 
#endif // EWOMS_FV_BASE_DISCRETIZATION_HH