stokesextensivequantities.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) 2011-2013 by Andreas Lauser
  Copyright (C) 2012 by Klaus Mosthaf

  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_STOKES_EXTENSIVE_QUANTITIES_HH
#define EWOMS_STOKES_EXTENSIVE_QUANTITIES_HH

#include "stokesproperties.hh"

#include <ewoms/models/common/energymodule.hh>
#include <ewoms/models/common/quantitycallbacks.hh>

#include <opm/material/common/Valgrind.hpp>

#include <dune/common/fvector.hh>

namespace Ewoms {

template <class TypeTag>
class StokesExtensiveQuantities
    : public EnergyExtensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy)>
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;

    enum { dimWorld = GridView::dimensionworld };
    enum { phaseIdx = GET_PROP_VALUE(TypeTag, StokesPhaseIndex) };
    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };

    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
    typedef Ewoms::EnergyExtensiveQuantities<TypeTag, enableEnergy> EnergyExtensiveQuantities;

public:
    static void registerParameters()
    {}

    void update(const ElementContext &elemCtx, int scvfIdx, int timeIdx, bool isBoundaryFace=false)
    {
        const auto &stencil = elemCtx.stencil(timeIdx);
        const auto &scvf =
            isBoundaryFace?
            stencil.boundaryFace(scvfIdx):
            stencil.interiorFace(scvfIdx);

        insideIdx_ = scvf.interiorIndex();
        outsideIdx_ = scvf.exteriorIndex();

        onBoundary_ = isBoundaryFace;
        normal_ = scvf.normal();
        Valgrind::CheckDefined(normal_);

        // calculate gradients and secondary variables at IPs
        const auto& gradCalc = elemCtx.gradientCalculator();
        PressureCallback<TypeTag> pressureCallback(elemCtx, phaseIdx);
        DensityCallback<TypeTag> densityCallback(elemCtx, phaseIdx);
        MolarDensityCallback<TypeTag> molarDensityCallback(elemCtx, phaseIdx);
        ViscosityCallback<TypeTag> viscosityCallback(elemCtx, phaseIdx);
        VelocityCallback<TypeTag> velocityCallback(elemCtx);
        VelocityComponentCallback<TypeTag> velocityComponentCallback(elemCtx);

        pressure_ = gradCalc.calculateValue(elemCtx, scvfIdx, pressureCallback);
        gradCalc.calculateGradient(pressureGrad_, elemCtx, scvfIdx, pressureCallback);
        density_ = gradCalc.calculateValue(elemCtx, scvfIdx, densityCallback);
        molarDensity_ = gradCalc.calculateValue(elemCtx, scvfIdx, molarDensityCallback);
        viscosity_ = gradCalc.calculateValue(elemCtx, scvfIdx, viscosityCallback);
        velocity_ = gradCalc.calculateValue(elemCtx, scvfIdx, velocityCallback);

        for (int dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            velocityComponentCallback.setDimIndex(dimIdx);
            gradCalc.calculateGradient(velocityGrad_[dimIdx],
                                       elemCtx,
                                       scvfIdx,
                                       velocityComponentCallback);
        }

        volumeFlux_ = (velocity_ * normal_);
        Valgrind::CheckDefined(volumeFlux_);

        // set the upstream and downstream DOFs
        upstreamIdx_ = scvf.interiorIndex();
        downstreamIdx_ = scvf.exteriorIndex();
        if (volumeFlux_ < 0)
            std::swap(upstreamIdx_, downstreamIdx_);

        EnergyExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx);

        Valgrind::CheckDefined(density_);
        Valgrind::CheckDefined(viscosity_);
        Valgrind::CheckDefined(velocity_);
        Valgrind::CheckDefined(pressureGrad_);
        Valgrind::CheckDefined(velocityGrad_);
    }

    template <class Context, class FluidState>
    void updateBoundary(const Context &context, int bfIdx, int timeIdx,
                        const FluidState &fluidState,
                        typename FluidSystem::ParameterCache &paramCache)
    {
        update(context, bfIdx, timeIdx, fluidState, paramCache,
               /*isOnBoundary=*/true);
    }

    Scalar pressure() const
    { return pressure_; }

    Scalar density() const
    { return density_; }

    Scalar molarDensity() const
    { return molarDensity_; }

    Scalar viscosity() const
    { return viscosity_; }

    const DimVector &pressureGrad() const
    { return pressureGrad_; }

    const DimVector &velocity() const
    { return velocity_; }

    const DimVector &velocityGrad(int axisIdx) const
    { return velocityGrad_[axisIdx]; }

    Scalar eddyViscosity() const
    { return 0; }

    Scalar eddyDiffusivity() const
    { return 0; }

    Scalar volumeFlux(int phaseIdx) const
    { return volumeFlux_; }

    Scalar upstreamWeight(int phaseIdx) const
    { return 1.0; }

    Scalar downstreamWeight(int phaseIdx) const
    { return 0.0; }

    int upstreamIndex(int phaseIdx) const
    { return upstreamIdx_; }

    int downstreamIndex(int phaseIdx) const
    { return downstreamIdx_; }

    int interiorIndex() const
    { return insideIdx_; }

    int exteriorIndex() const
    { return outsideIdx_; }

    bool onBoundary() const
    { return onBoundary_; }

    Scalar extrusionFactor() const
    { return 1.0; }

    const DimVector &normal() const
    { return normal_; }

private:
    bool onBoundary_;

    // values at the integration point
    Scalar density_;
    Scalar molarDensity_;
    Scalar viscosity_;
    Scalar pressure_;
    Scalar volumeFlux_;
    DimVector velocity_;
    DimVector normal_;

    // gradients at the IPs
    DimVector pressureGrad_;
    DimVector velocityGrad_[dimWorld];

    // local index of the upstream dof
    int upstreamIdx_;
    // local index of the downstream dof
    int downstreamIdx_;

    int insideIdx_;
    int outsideIdx_;
};

} // namespace Ewoms

#endif