fvbaseelementcontext.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_ELEMENT_CONTEXT_HH
#define EWOMS_FV_BASE_ELEMENT_CONTEXT_HH

#include <dune/common/fvector.hh>

#include <opm/models/discretization/common/fvbaseparameters.hh>
#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/discretization/common/linearizationtype.hh>

#include <opm/models/utils/alignedallocator.hh>

#include <array>
#include <cassert>
#include <cstddef>
#include <stdexcept>
#include <vector>

namespace Opm {

template<class TypeTag>
class FvBaseElementContext
{
    using Implementation = GetPropType<TypeTag, Properties::ElementContext>;

    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;

    // the history size of the time discretization in number of steps
    enum { timeDiscHistorySize = getPropValue<TypeTag, Properties::TimeDiscHistorySize>() };

    struct DofStore_ {
        std::array<IntensiveQuantities, timeDiscHistorySize> intensiveQuantities;
        std::array<const PrimaryVariables*, timeDiscHistorySize> priVars;
        std::array<const IntensiveQuantities*, timeDiscHistorySize> thermodynamicHint;
    };
    using ExtensiveQuantitiesVector = std::vector<ExtensiveQuantities>;

    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using Model = GetPropType<TypeTag, Properties::Model>;
    using Stencil = GetPropType<TypeTag, Properties::Stencil>;
    using GradientCalculator = GetPropType<TypeTag, Properties::GradientCalculator>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;

    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Element = typename GridView::template Codim<0>::Entity;

    static const unsigned dimWorld = GridView::dimensionworld;

    using CoordScalar = typename GridView::ctype;
    using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;

    // we don't allow copies of element contexts!
    FvBaseElementContext(const FvBaseElementContext& ) = delete;

public:
    explicit FvBaseElementContext(const Simulator& simulator)
        : gridView_(simulator.gridView())
        , stencil_(gridView_, simulator.model().dofMapper() )
    {
        // remember the simulator object
        simulatorPtr_ = &simulator;
        enableStorageCache_ = Parameters::Get<Parameters::EnableStorageCache>();
    }

    static void *operator new(std::size_t size)
    { return aligned_alloc(alignof(FvBaseElementContext), size); }

    static void operator delete(void *ptr)
    { aligned_free(ptr); }

    void updateAll(const Element& elem)
    {
        asImp_().updateStencil(elem);
        asImp_().updateAllIntensiveQuantities();
        asImp_().updateAllExtensiveQuantities();
    }

    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());
    }

    void updatePrimaryStencil(const Element& elem)
    {
        // remember the current element
        elemPtr_ = &elem;

        // update the finite element geometry
        stencil_.updatePrimaryTopology(elem);

        dofVars_.resize(stencil_.numPrimaryDof());
    }

    void updateStencilTopology(const Element& elem)
    {
        // remember the current element
        elemPtr_ = &elem;

        // update the finite element geometry
        stencil_.updateTopology(elem);
    }

    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);
        }
    }

    void updateIntensiveQuantities(unsigned timeIdx)
    { updateIntensiveQuantities_(timeIdx, numDof(timeIdx)); }

    void updatePrimaryIntensiveQuantities(unsigned timeIdx)
    { updateIntensiveQuantities_(timeIdx, numPrimaryDof(timeIdx)); }

    void updateIntensiveQuantities(const PrimaryVariables& priVars, unsigned dofIdx, unsigned timeIdx)
    { asImp_().updateSingleIntQuants_(priVars, dofIdx, timeIdx); }

    void updateAllExtensiveQuantities()
    { asImp_().updateExtensiveQuantities(/*timeIdx=*/0); }

    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);
        }
    }

    void setFocusDofIndex(unsigned dofIdx)
    { focusDofIdx_ = dofIdx; }

    unsigned focusDofIndex() const
    { return focusDofIdx_; }

    LinearizationType linearizationType() const
    { return this->model().linearizer().getLinearizationType(); }

    const Simulator& simulator() const
    { return *simulatorPtr_; }

    const Problem& problem() const
    { return simulatorPtr_->problem(); }

    const Model& model() const
    { return simulatorPtr_->model(); }

    const GridView& gridView() const
    { return gridView_; }

    const Element& element() const
    { return *elemPtr_; }

    std::size_t numDof(unsigned timeIdx) const
    { return stencil(timeIdx).numDof(); }

    std::size_t numPrimaryDof(unsigned timeIdx) const
    { return stencil(timeIdx).numPrimaryDof(); }

    std::size_t numInteriorFaces(unsigned timeIdx) const
    { return stencil(timeIdx).numInteriorFaces(); }

    std::size_t numBoundaryFaces(unsigned timeIdx) const
    { return stencil(timeIdx).numBoundaryFaces(); }

    const Stencil& stencil(unsigned) const
    { return stencil_; }

    decltype(auto) pos(unsigned dofIdx, unsigned) const
    { return stencil_.subControlVolume(dofIdx).globalPos(); }

    unsigned globalSpaceIndex(unsigned dofIdx, unsigned timeIdx) const
    { return stencil(timeIdx).globalSpaceIndex(dofIdx); }

    Scalar dofVolume(unsigned dofIdx, unsigned timeIdx) const
    { return stencil(timeIdx).subControlVolume(dofIdx).volume();  }

    Scalar dofTotalVolume(unsigned dofIdx, unsigned timeIdx) const
    { return model().dofTotalVolume(globalSpaceIndex(dofIdx, timeIdx)); }

    bool onBoundary() const
    { return element().hasBoundaryIntersections(); }

    const IntensiveQuantities& intensiveQuantities(unsigned dofIdx, unsigned timeIdx) const
    {
#ifndef NDEBUG
        assert(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];
    }

    const IntensiveQuantities *thermodynamicHint(unsigned dofIdx, unsigned timeIdx) const
    {
        assert(dofIdx < numDof(timeIdx));
        return dofVars_[dofIdx].thermodynamicHint[timeIdx];
    }

    IntensiveQuantities& intensiveQuantities(unsigned dofIdx, unsigned timeIdx)
    {
        assert(dofIdx < numDof(timeIdx));
        return dofVars_[dofIdx].intensiveQuantities[timeIdx];
    }

    const PrimaryVariables& primaryVars(unsigned dofIdx, unsigned timeIdx) const
    {
        assert(dofIdx < numDof(timeIdx));
        return *dofVars_[dofIdx].priVars[timeIdx];
    }

    bool haveStashedIntensiveQuantities() const
    { return stashedDofIdx_ != -1; }

    int stashedDofIdx() const
    { return stashedDofIdx_; }

    void stashIntensiveQuantities(unsigned dofIdx)
    {
        assert(dofIdx < numDof(/*timeIdx=*/0));

        intensiveQuantitiesStashed_ = dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0];
        priVarsStashed_ = *dofVars_[dofIdx].priVars[/*timeIdx=*/0];
        stashedDofIdx_ = static_cast<int>(dofIdx);
    }

    void restoreIntensiveQuantities(unsigned dofIdx)
    {
        dofVars_[dofIdx].priVars[/*timeIdx=*/0] = &priVarsStashed_;
        dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0] = intensiveQuantitiesStashed_;
        stashedDofIdx_ = -1;
    }

    const GradientCalculator& gradientCalculator() const
    { return gradientCalculator_; }

    const ExtensiveQuantities& extensiveQuantities(unsigned fluxIdx, unsigned) const
    { return extensiveQuantities_[fluxIdx]; }

    bool enableStorageCache() const
    { return enableStorageCache_; }

    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:
    void updateIntensiveQuantities_(unsigned timeIdx, std::size_t numDof)
    {
        // update the intensive quantities for the whole history
        const SolutionVector& globalSol = model().solution(timeIdx);

        const unsigned intensiveHistorySize = model().cachedIntensiveQuantityHistorySize();

        // 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;

            if (timeIdx >= intensiveHistorySize) {
                // If we do not have a cache for this time index, it does not make sense to try to fetch it or
                // update the intensive quantities cache for this time index.
                updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
            }
            else {
                // If we have a cache for this time index, we can fetch the intensive quantities from the cache.
                // We can then also fetch the thermodynamic hint from the cache.
                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_;

    template<class T> using AlignedVector = std::vector<T, aligned_allocator<T, alignof(T)>>;

    AlignedVector<DofStore_> dofVars_;
    AlignedVector<ExtensiveQuantities> extensiveQuantities_;

    const Simulator* simulatorPtr_{};
    const Element* elemPtr_{};
    const GridView gridView_;
    Stencil stencil_;

    int stashedDofIdx_{-1};
    int focusDofIdx_{-1};
    bool enableStorageCache_{false};
};

} // namespace Opm

#endif