diffusionmodule.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

  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_DIFFUSION_MODULE_HH
#define EWOMS_DIFFUSION_MODULE_HH

#include <ewoms/disc/common/fvbaseproperties.hh>
#include <ewoms/models/common/quantitycallbacks.hh>

#include <dune/common/fvector.hh>

namespace Ewoms {
namespace Properties {
NEW_PROP_TAG(Indices);
}

template <class TypeTag, bool enableDiffusion>
class DiffusionModule;

template <class TypeTag>
class DiffusionModule<TypeTag, /*enableDiffusion=*/false>
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
    typedef typename FluidSystem::ParameterCache ParameterCache;

public:
    static void registerParameters()
    {}

    template <class Context>
    static void addDiffusiveFlux(RateVector &flux, const Context &context,
                                 int spaceIdx, int timeIdx)
    {}
};

template <class TypeTag>
class DiffusionModule<TypeTag, /*enableDiffusion=*/true>
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;

    enum { numPhases = FluidSystem::numPhases };
    enum { numComponents = FluidSystem::numComponents };
    enum { conti0EqIdx = Indices::conti0EqIdx };

    typedef Opm::MathToolbox<Evaluation> Toolbox;

public:
    static void registerParameters()
    {}

    template <class Context>
    static void addDiffusiveFlux(RateVector &flux, const Context &context,
                                 int spaceIdx, int timeIdx)
    {
        const auto &extQuants = context.extensiveQuantities(spaceIdx, timeIdx);

        const auto &fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
        const auto &fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();

        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            // arithmetic mean of the phase's molar density
            Evaluation rhoMolar = fluidStateI.molarDensity(phaseIdx);
            rhoMolar += Toolbox::value(fluidStateJ.molarDensity(phaseIdx));
            rhoMolar /= 2;

            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
                // mass flux due to molecular diffusion
                flux[conti0EqIdx + compIdx] +=
                    -rhoMolar
                    * extQuants.moleFractionGradientNormal(phaseIdx, compIdx)
                    * extQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx);
        }
    }
};

template <class TypeTag, bool enableDiffusion>
class DiffusionIntensiveQuantities;

template <class TypeTag>
class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;

    typedef typename FluidSystem::ParameterCache ParameterCache;

public:
    Scalar tortuosity(int phaseIdx) const
    {
        OPM_THROW(std::logic_error, "Method tortuosity() does not make sense "
                                    "if diffusion is disabled");
    }

    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
    {
        OPM_THROW(std::logic_error, "Method diffusionCoefficient() does not "
                                    "make sense if diffusion is disabled");
    }

    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIdx) const
    {
        OPM_THROW(std::logic_error, "Method effectiveDiffusionCoefficient() "
                                    "does not make sense if diffusion is "
                                    "disabled");
    }

protected:
    template <class FluidState>
    void update_(FluidState &fs,
                 ParameterCache &paramCache,
                 const ElementContext &elemCtx,
                 int dofIdx,
                 int timeIdx)
    { }
};

template <class TypeTag>
class DiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;

    typedef typename FluidSystem::ParameterCache ParameterCache;

    enum { numPhases = FluidSystem::numPhases };
    enum { numComponents = FluidSystem::numComponents };

public:
    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
    { return diffusionCoefficient_[phaseIdx][compIdx]; }

    Scalar tortuosity(int phaseIdx) const
    { return tortuosity_[phaseIdx]; }

    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIdx) const
    { return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; }

protected:
    template <class FluidState>
    void update_(FluidState &fluidState,
                 ParameterCache &paramCache,
                 const ElementContext &elemCtx,
                 int dofIdx,
                 int timeIdx)
    {
        const auto &intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);

        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!elemCtx.model().phaseIsConsidered(phaseIdx))
                continue;

            // TODO: let the problem do this (this is a constitutive
            // relation of which the model should be free of from the
            // abstraction POV!)
            tortuosity_[phaseIdx] =
                1.0 / (intQuants.porosity() * intQuants.porosity())
                * std::pow(std::max(0.0001,
                                    intQuants.porosity()
                                    * intQuants.fluidState().saturation(phaseIdx)),
                           7.0 / 3);

            for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                diffusionCoefficient_[phaseIdx][compIdx] =
                    FluidSystem::diffusionCoefficient(fluidState,
                                                      paramCache,
                                                      phaseIdx,
                                                      compIdx);
            }
        }
    }

private:
    Scalar tortuosity_[numPhases];
    Scalar diffusionCoefficient_[numPhases][numComponents];
};

template <class TypeTag, bool enableDiffusion>
class DiffusionExtensiveQuantities;

template <class TypeTag>
class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/false>
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;

protected:
    void update_(const ElementContext &elemCtx, int faceIdx, int timeIdx)
    {}

    template <class Context, class FluidState>
    void updateBoundary_(const Context &context, int bfIdx, int timeIdx,
                         const FluidState &fluidState)
    {}

public:
    Scalar moleFractionGradientNormal(int phaseIdx, int compIdx) const
    {
        OPM_THROW(std::logic_error, "Method moleFractionGradient() does not "
                                    "make sense if diffusion is disabled.");
    }

    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIdx) const
    {
        OPM_THROW(std::logic_error, "Method effectiveDiffusionCoefficient() "
                                    "does not make sense if diffusion is "
                                    "disabled.");
    }
};

template <class TypeTag>
class DiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;

    enum { dimWorld = GridView::dimensionworld };
    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
    enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };

    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;

protected:
    void update_(const ElementContext &elemCtx, int faceIdx, int timeIdx)
    {
        const auto& gradCalc = elemCtx.gradientCalculator();
        Ewoms::MoleFractionCallback<TypeTag> moleFractionCallback(elemCtx);

        DimVector temperatureGrad;

        const auto &face = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
        const auto &extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);

        const auto &intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
        const auto &intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);

        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!elemCtx.model().phaseIsConsidered(phaseIdx))
                continue;

            moleFractionCallback.setPhaseIndex(phaseIdx);
            for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                moleFractionCallback.setComponentIndex(compIdx);

                DimVector moleFractionGradient(0.0);
                gradCalc.calculateGradient(moleFractionGradient,
                                           elemCtx,
                                           faceIdx,
                                           moleFractionCallback);

                moleFractionGradientNormal_[phaseIdx][compIdx] =
                    face.normal()*moleFractionGradient;
                Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);

                // use the arithmetic average for the effective
                // diffusion coefficients.
                effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
                    (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
                     +
                     intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
                    / 2;

                Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
            }
        }
    }

    template <class Context, class FluidState>
    void updateBoundary_(const Context &context, int bfIdx, int timeIdx,
                         const FluidState &fluidState)
    {
        const auto &stencil = context.stencil(timeIdx);
        const auto &face = stencil.boundaryFace(bfIdx);

        const auto &elemCtx = context.elementContext();
        int insideScvIdx = face.interiorIndex();
        const auto &insideScv = stencil.subControlVolume(insideScvIdx);

        const auto &intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
        const auto &fluidStateInside = intQuantsInside.fluidState();

        // distance between the center of the SCV and center of the boundary face
        DimVector distVec = face.integrationPos();
        distVec -= context.element().geometry().global(insideScv.localGeometry().center());

        Scalar dist = distVec * face.normal();

        // if the following assertation triggers, the center of the
        // center of the interior SCV was not inside the element!
        assert(dist > 0);

        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!elemCtx.model().phaseIsConsidered(phaseIdx))
                continue;

            for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                DimVector moleFractionGradient(0.0);

                // calculate mole fraction gradient using two-point
                // gradients
                moleFractionGradientNormal_[phaseIdx][compIdx] =
                    (fluidState.moleFraction(phaseIdx, compIdx)
                     -
                     fluidStateInside.moleFraction(phaseIdx, compIdx))
                    / dist;
                Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);

                // use effective diffusion coefficients of the interior finite
                // volume.
                effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
                    intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx);

                Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
            }
        }
    }

public:
    Scalar moleFractionGradientNormal(int phaseIdx, int compIdx) const
    { return moleFractionGradientNormal_[phaseIdx][compIdx]; }

    Scalar effectiveDiffusionCoefficient(int phaseIdx, int compIdx) const
    { return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }

private:
    Scalar moleFractionGradientNormal_[numPhases][numComponents];
    Scalar effectiveDiffusionCoefficient_[numPhases][numComponents];
};

} // namespace Ewoms

#endif