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:
/*
  Copyright (C) 2008-2013 by Andreas Lauser

  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/>.
*/
#ifndef EWOMS_FV_BASE_ELEMENT_CONTEXT_HH
#define EWOMS_FV_BASE_ELEMENT_CONTEXT_HH

#include "fvbaseproperties.hh"

#include <dune/common/fvector.hh>

#include <vector>

namespace Ewoms {

template<class TypeTag>
class FvBaseElementContext
{
    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 int dim = GridView::dimension;
    static const int numEq = GET_PROP_VALUE(TypeTag, NumEq);
    static const int requireScvCenterGradients =
        GET_PROP_VALUE(TypeTag, RequireScvCenterGradients);

    typedef typename GridView::ctype CoordScalar;
    typedef Dune::FieldVector<CoordScalar, dim> GlobalPosition;

    // we don't allow copies of element contexts!
    FvBaseElementContext(const FvBaseElementContext &context)
    {}

public:
    explicit FvBaseElementContext(const Simulator &simulator)
        : gridView_(simulator.gridView())
        , stencil_(gridView_)
    {
        // remember the simulator object
        simulatorPtr_ = &simulator;
    }

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

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

        // update the stencil. the center gradients are kind of expensive to calculate
        // and most models don't need them, so that we only do this if the model needs
        // them
        stencil_.update(elem);
        if (requireScvCenterGradients)
            stencil_.updateCenterGradients();

        // resize the arrays containing the flux and the volume variables
        dofVars_.resize(stencil_.numDof());
        extensiveQuantities_.resize(stencil_.numInteriorFaces());
    }

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

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

    void updateAllIntensiveQuantities()
    {
        for (int timeIdx = 0; timeIdx < timeDiscHistorySize; ++ timeIdx)
            updateIntensiveQuantities(timeIdx);
    }

    void updateIntensiveQuantities(int timeIdx)
    { updateIntensiveQuantities_(timeIdx, numDof(/*timeIdx=*/0)); }

    void updatePrimaryIntensiveQuantities(int timeIdx)
    { updateIntensiveQuantities_(timeIdx, numPrimaryDof(/*timeIdx=*/0)); }

    void updateIntensiveQuantities(const PrimaryVariables &priVars, int dofIdx, int timeIdx)
    {
        updateSingleIntQuants_(priVars, dofIdx, timeIdx);

        // update gradients inside a sub control volume
        int nDof = numDof(/*timeIdx=*/0);
        for (int gradDofIdx = 0; gradDofIdx < nDof; gradDofIdx++) {
            dofVars_[gradDofIdx].intensiveQuantities[timeIdx].updateScvGradients(/*context=*/*this,
                                                                                 gradDofIdx,
                                                                                 timeIdx);
        }
    }

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

    void updateExtensiveQuantities(int timeIdx)
    {
        gradientCalculator_.prepare(/*context=*/*this, timeIdx);

        for (int fluxIdx = 0; fluxIdx < numInteriorFaces(timeIdx); fluxIdx++) {
            extensiveQuantities_[fluxIdx].update(/*context=*/ *this,
                                                 /*localIndex=*/fluxIdx,
                                                 timeIdx);
        }
    }

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

    int numDof(int timeIdx) const
    { return stencil(timeIdx).numDof(); }

    int numPrimaryDof(int timeIdx) const
    { return stencil(timeIdx).numPrimaryDof(); }

    int numInteriorFaces(int timeIdx) const
    { return stencil(timeIdx).numInteriorFaces(); }

    int numBoundaryFaces(int timeIdx) const
    { return stencil(timeIdx).numBoundaryFaces(); }

    const Stencil &stencil(int timeIdx) const
    { return stencil_; }

    const GlobalPosition &pos(int dofIdx, int timeIdx) const
    { return stencil_.subControlVolume(dofIdx).globalPos(); }

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


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

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

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

    const IntensiveQuantities &intensiveQuantities(int dofIdx, int timeIdx) const
    {
        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
        return dofVars_[dofIdx].intensiveQuantities[timeIdx];
    }

    const IntensiveQuantities *thermodynamicHint(int dofIdx, int timeIdx) const
    {
        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
        return dofVars_[dofIdx].thermodynamicHint[timeIdx];
    }
    IntensiveQuantities &intensiveQuantities(int dofIdx, int timeIdx)
    {
        assert(0 <= dofIdx && dofIdx < numDof(timeIdx));
        return dofVars_[dofIdx].intensiveQuantities[timeIdx];
    }

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

    void saveIntensiveQuantities(int dofIdx)
    {
        assert(0 <= dofIdx && dofIdx < numDof(/*timeIdx=*/0));

        intensiveQuantitiesSaved_ = dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0];
        priVarsSaved_ = dofVars_[dofIdx].priVars[/*timeIdx=*/0];
    }

    void restoreIntensiveQuantities(int dofIdx)
    {
        dofVars_[dofIdx].priVars[/*timeIdx=*/0] = priVarsSaved_;
        dofVars_[dofIdx].intensiveQuantities[/*timeIdx=*/0] = intensiveQuantitiesSaved_;
    }

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

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

protected:
    void updateIntensiveQuantities_(int timeIdx, int numDof)
    {
        // update the intensive quantities for the whole history
        const SolutionVector &globalSol = model().solution(timeIdx);

        // update the non-gradient quantities
        for (int dofIdx = 0; dofIdx < numDof; dofIdx++) {
            int globalIdx = globalSpaceIndex(dofIdx, timeIdx);
            const PrimaryVariables& dofSol = globalSol[globalIdx];

            dofVars_[dofIdx].thermodynamicHint[timeIdx] =
                model().thermodynamicHint(globalIdx, timeIdx);

            const auto *cachedIntQuants = model().cachedIntensiveQuantities(globalIdx, timeIdx);
            if (cachedIntQuants) {
                dofVars_[dofIdx].priVars[timeIdx] = dofSol;
                dofVars_[dofIdx].intensiveQuantities[timeIdx] = *cachedIntQuants;
            }
            else {
                updateSingleIntQuants_(dofSol, dofIdx, timeIdx);
                model().updateCachedIntensiveQuantities(dofVars_[dofIdx].intensiveQuantities[timeIdx],
                                                        globalIdx,
                                                        timeIdx);
            }
        }

        // update gradients
        for (int dofIdx = 0; dofIdx < numDof; dofIdx++) {
            dofVars_[dofIdx].intensiveQuantities[timeIdx].updateScvGradients(/*context=*/*this,
                                                                             dofIdx,
                                                                             timeIdx);
        }
    }

    void updateSingleIntQuants_(const PrimaryVariables &priVars, int dofIdx, int timeIdx)
    {
        dofVars_[dofIdx].priVars[timeIdx] = priVars;
        dofVars_[dofIdx].intensiveQuantities[timeIdx].update(/*context=*/*this, dofIdx, timeIdx);
    }

    DofVarsVector dofVars_;

    IntensiveQuantities intensiveQuantitiesSaved_;
    PrimaryVariables priVarsSaved_;

    GradientCalculator gradientCalculator_;

    ExtensiveQuantitiesVector extensiveQuantities_;

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

} // namespace Ewoms

#endif