forchheimerfluxmodule.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_FORCHHEIMER_FLUX_MODULE_HH
#define EWOMS_FORCHHEIMER_FLUX_MODULE_HH

#include <ewoms/disc/common/fvbaseproperties.hh>
#include "darcyfluxmodule.hh"


#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <cmath>

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

namespace Ewoms {
template <class TypeTag>
class ForchheimerIntensiveQuantities;

template <class TypeTag>
class ForchheimerExtensiveQuantities;

template <class TypeTag>
class ForchheimerBaseProblem;

template <class TypeTag>
struct ForchheimerFluxModule
{
    typedef ForchheimerIntensiveQuantities<TypeTag> FluxIntensiveQuantities;
    typedef ForchheimerExtensiveQuantities<TypeTag> FluxExtensiveQuantities;
    typedef ForchheimerBaseProblem<TypeTag> FluxBaseProblem;

    static void registerParameters()
    {}
};

template <class TypeTag>
class ForchheimerBaseProblem
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;

public:
    template <class Context>
    Scalar ergunCoefficient(Context &context, int spaceIdx, int timeIdx) const
    {
        OPM_THROW(std::logic_error,
                  "Not implemented: Problem::ergunCoefficient()");
    }

    template <class Context>
    Scalar mobilityPassabilityRatio(Context &context, int spaceIdx, int timeIdx,
                                    int phaseIdx) const
    {
        return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(
                         phaseIdx);
    }
};

template <class TypeTag>
class ForchheimerIntensiveQuantities
{
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;

    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };

public:
    Scalar ergunCoefficient() const
    {
        return ergunCoefficient_;
    }

    Scalar mobilityPassabilityRatio(int phaseIdx) const
    {
        return mobilityPassabilityRatio_[phaseIdx];
    }

protected:
    void update_(const ElementContext &elemCtx, int dofIdx, int timeIdx)
    {
        const auto &problem = elemCtx.problem();
        ergunCoefficient_ = problem.ergunCoefficient(elemCtx, dofIdx, timeIdx);

        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
            mobilityPassabilityRatio_[phaseIdx] =
                problem.mobilityPassabilityRatio(elemCtx,
                                                 dofIdx,
                                                 timeIdx,
                                                 phaseIdx);
    }

private:
    Scalar ergunCoefficient_;
    Scalar mobilityPassabilityRatio_[numPhases];
};

template <class TypeTag>
class ForchheimerExtensiveQuantities
    : public DarcyExtensiveQuantities<TypeTag>
{
    typedef DarcyExtensiveQuantities<TypeTag> DarcyExtQuants;
    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
    typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;


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

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

public:
    Scalar ergunCoefficient() const
    { return ergunCoefficient_; }

    Scalar mobilityPassabilityRatio(int phaseIdx) const
    { return mobilityPassabilityRatio_[phaseIdx]; }

protected:
    void calculateGradients_(const ElementContext &elemCtx,
                             int faceIdx,
                             int timeIdx)
    {
        DarcyExtQuants::calculateGradients_(elemCtx, faceIdx, timeIdx);

        const auto &intQuantsIn = elemCtx.intensiveQuantities(this->interiorDofIdx_, timeIdx);
        const auto &intQuantsEx = elemCtx.intensiveQuantities(this->exteriorDofIdx_, timeIdx);

        // calculate the square root of the intrinsic permeability
        assert(isDiagonal_(this->K_));
        sqrtK_ = 0.0;
        for (int i = 0; i < dimWorld; ++i)
            sqrtK_[i][i] = std::sqrt(this->K_[i][i]);

        // obtain the Ergun coefficient. Lacking better ideas, we use its the arithmetic mean.
        ergunCoefficient_ = (intQuantsIn.ergunCoefficient() + intQuantsEx.ergunCoefficient())/2;

        // obtain the mobility to passability ratio for each phase.
        for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
            if (!elemCtx.model().phaseIsConsidered(phaseIdx))
                continue;

            const auto &up = elemCtx.intensiveQuantities(this->upstreamIndex_(phaseIdx), timeIdx);

            density_[phaseIdx] = up.fluidState().density(phaseIdx);
            mobilityPassabilityRatio_[phaseIdx] = up.mobilityPassabilityRatio(phaseIdx);
        }
    }

    template <class FluidState>
    void calculateBoundaryGradients_(const ElementContext &elemCtx,
                                     int boundaryFaceIdx,
                                     int timeIdx,
                                     const FluidState& fluidState,
                                     const typename FluidSystem::ParameterCache& paramCache)
    {
        DarcyExtQuants::calculateBoundaryGradients_(elemCtx,
                                                    boundaryFaceIdx,
                                                    timeIdx,
                                                    fluidState,
                                                    paramCache);

        const auto &intQuantsIn = elemCtx.intensiveQuantities(this->interiorDofIdx_, timeIdx);

        // obtain the Ergun coefficient. Because we are on the boundary here, we will
        // take the Ergun coefficient of the interior
        ergunCoefficient_ = intQuantsIn.ergunCoefficient();

        // calculate the square root of the intrinsic permeability
        assert(isDiagonal_(this->K_));
        sqrtK_ = 0.0;
        for (int i = 0; i < dimWorld; ++i)
            sqrtK_[i][i] = std::sqrt(this->K_[i][i]);

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

            density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx);
            mobilityPassabilityRatio_[phaseIdx] = intQuantsIn.mobilityPassabilityRatio(phaseIdx);
        }
    }

    void calculateFluxes_(const ElementContext &elemCtx, int scvfIdx, int timeIdx)
    {
        const auto &intQuantsI = elemCtx.intensiveQuantities(asImp_().interiorIndex(), timeIdx);
        const auto &intQuantsJ = elemCtx.intensiveQuantities(asImp_().exteriorIndex(), timeIdx);

        const auto &scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
        const auto &normal = scvf.normal();
        Valgrind::CheckDefined(normal);

        // obtain the Ergun coefficient from the intensive quantity object. Until a
        // better method comes along, we use arithmetic averaging.
        ergunCoefficient_ = (intQuantsI.ergunCoefficient() + intQuantsJ.ergunCoefficient()) / 2;

        // calculate the weights of the upstream and the downstream control volumes
        for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
            if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
                this->filterVelocity_[phaseIdx] = 0;
                this->volumeFlux_[phaseIdx] = 0;
                continue;
            }

            calculateForchheimerFlux_(phaseIdx);
            this->volumeFlux_[phaseIdx] = (this->filterVelocity_[phaseIdx] * normal);
        }
    }

    void calculateBoundaryFluxes_(const ElementContext& elemCtx,
                                      int bfIdx,
                                      int timeIdx)
    {
        const auto &boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx);
        const auto &normal = boundaryFace.normal();

        // calculate the weights of the upstream and the downstream degrees of freedom
        for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) {
            if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
                this->filterVelocity_[phaseIdx] = 0;
                this->volumeFlux_[phaseIdx] = 0;
                continue;
            }

            calculateForchheimerFlux_(phaseIdx);
            this->volumeFlux_[phaseIdx] = (this->filterVelocity_[phaseIdx] * normal);
        }
    }

    void calculateForchheimerFlux_(int phaseIdx)
    {
        // initial guess: filter velocity is zero
        DimVector &velocity = this->filterVelocity_[phaseIdx];
        velocity = 0;

        // the change of velocity between two consecutive Newton iterations
        DimVector deltaV(1e5);
        // the function value that is to be minimized of the equation that is to be
        // fulfilled
        DimVector residual;
        // derivative of equation that is to be solved
        DimMatrix gradResid;

        // search by means of the Newton method for a root of Forchheimer equation
        int newtonIter = 0;
        while (deltaV.two_norm() > 1e-11) {
            if (newtonIter >= 50)
                OPM_THROW(Opm::NumericalProblem,
                          "Could not determine Forchheimer velocity within "
                          << newtonIter << " iterations");
            ++newtonIter;

            // calculate the residual and its Jacobian matrix
            gradForchheimerResid_(residual, gradResid, phaseIdx);

            // newton method
            gradResid.solve(deltaV, residual);
            velocity -= deltaV;
        }
    }

    void forchheimerResid_(DimVector &residual, int phaseIdx) const
    {
        const DimVector &velocity = this->filterVelocity_[phaseIdx];

        // Obtaining the upstreamed quantities
        Scalar mobility = this->mobility_[phaseIdx];
        Scalar density = density_[phaseIdx];
        Scalar mobilityPassabilityRatio = mobilityPassabilityRatio_[phaseIdx];

        // optain the quantites for the integration point
        const auto &pGrad = this->potentialGrad_[phaseIdx];

        // residual = v_\alpha
        residual = velocity;

        // residual += mobility_\alpha K(\grad p_\alpha - \rho_\alpha g)
        this->K_.usmv(mobility, pGrad, residual);

        // Forchheimer turbulence correction:
        //
        // residual +=
        //   \rho_\alpha
        //   * mobility_\alpha
        //   * C_E / \eta_{r,\alpha}
        //   * abs(v_\alpha) * sqrt(K)*v_\alpha
        sqrtK_.usmv(density*mobilityPassabilityRatio*ergunCoefficient_*velocity.two_norm(),
                    velocity,
                    residual);

        Valgrind::CheckDefined(residual);
    }

    void gradForchheimerResid_(DimVector &residual, DimMatrix &gradResid,
                               int phaseIdx)
    {
        DimVector &velocity = this->filterVelocity_[phaseIdx];
        forchheimerResid_(residual, phaseIdx);

        Scalar eps = 1e-11;
        DimVector tmp;
        for (int i = 0; i < dimWorld; ++i) {
            Scalar coordEps = std::max(eps, velocity[i] * (1 + eps));
            velocity[i] += coordEps;
            forchheimerResid_(tmp, phaseIdx);
            tmp -= residual;
            tmp /= coordEps;
            gradResid[i] = tmp;
            velocity[i] -= coordEps;
        }
    }

    bool isDiagonal_(const DimMatrix &K) const
    {
        for (int i = 0; i < dimWorld; i++) {
            for (int j = 0; j < dimWorld; j++) {
                if (i == j)
                    continue;

                if (std::abs(K[i][j]) > 1e-25)
                    return false;
            }
        }
        return true;
    }

private:
    Implementation &asImp_()
    { return *static_cast<Implementation *>(this); }

    const Implementation &asImp_() const
    { return *static_cast<const Implementation *>(this); }

protected:
    // intrinsic permeability tensor and its square root
    DimMatrix sqrtK_;

    // Ergun coefficient of all phases at the integration point
    Scalar ergunCoefficient_;

    // Passability of all phases at the integration point
    Scalar mobilityPassabilityRatio_[numPhases];

    // Density of all phases at the integration point
    Scalar density_[numPhases];
};

} // namespace Ewoms

#endif