ncpmodel.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 Philipp Nuske

  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_NCP_MODEL_HH
#define EWOMS_NCP_MODEL_HH

#include <opm/material/localad/Math.hpp>

#include "ncpproperties.hh"
#include "ncplocalresidual.hh"
#include "ncpextensivequantities.hh"
#include "ncpprimaryvariables.hh"
#include "ncpboundaryratevector.hh"
#include "ncpratevector.hh"
#include "ncpintensivequantities.hh"
#include "ncpnewtonmethod.hh"
#include "ncpindices.hh"

#include <ewoms/models/common/multiphasebasemodel.hh>
#include <ewoms/models/common/energymodule.hh>
#include <ewoms/models/common/diffusionmodule.hh>
#include <ewoms/io/vtkcompositionmodule.hh>
#include <ewoms/io/vtkenergymodule.hh>
#include <ewoms/io/vtkdiffusionmodule.hh>

#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>

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

#include <sstream>
#include <string>
#include <vector>
#include <array>

namespace Ewoms {
template <class TypeTag>
class NcpModel;
}

namespace Ewoms {
namespace Properties {
NEW_TYPE_TAG(NcpModel, INHERITS_FROM(MultiPhaseBaseModel,
                                     VtkComposition,
                                     VtkEnergy,
                                     VtkDiffusion));

SET_TYPE_PROP(NcpModel,
              LocalResidual,
              Ewoms::NcpLocalResidual<TypeTag>);

SET_TYPE_PROP(NcpModel, NewtonMethod, Ewoms::NcpNewtonMethod<TypeTag>);

SET_TYPE_PROP(NcpModel, Model, Ewoms::NcpModel<TypeTag>);

SET_TYPE_PROP(NcpModel, BaseProblem, Ewoms::MultiPhaseBaseProblem<TypeTag>);

SET_BOOL_PROP(NcpModel, EnableEnergy, false);

SET_BOOL_PROP(NcpModel, EnableDiffusion, false);

SET_TYPE_PROP(NcpModel, RateVector, Ewoms::NcpRateVector<TypeTag>);

SET_TYPE_PROP(NcpModel, BoundaryRateVector, Ewoms::NcpBoundaryRateVector<TypeTag>);

SET_TYPE_PROP(NcpModel, PrimaryVariables, Ewoms::NcpPrimaryVariables<TypeTag>);

SET_TYPE_PROP(NcpModel, IntensiveQuantities, Ewoms::NcpIntensiveQuantities<TypeTag>);

SET_TYPE_PROP(NcpModel, ExtensiveQuantities, Ewoms::NcpExtensiveQuantities<TypeTag>);

SET_INT_PROP(NcpModel, NcpNewtonNumChoppedIterations, 3);

SET_TYPE_PROP(NcpModel, Indices, Ewoms::NcpIndices<TypeTag, 0>);

SET_SCALAR_PROP(NcpModel, NcpPressureBaseWeight, 1.0);
SET_SCALAR_PROP(NcpModel, NcpSaturationsBaseWeight, 1.0);
SET_SCALAR_PROP(NcpModel, NcpFugacitiesBaseWeight, 1.0);

} // namespace Properties

template <class TypeTag>
class NcpModel
    : public MultiPhaseBaseModel<TypeTag>
{
    typedef MultiPhaseBaseModel<TypeTag> ParentType;

    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
    typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
    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 { fugacity0Idx = Indices::fugacity0Idx };
    enum { pressure0Idx = Indices::pressure0Idx };
    enum { saturation0Idx = Indices::saturation0Idx };
    enum { conti0EqIdx = Indices::conti0EqIdx };
    enum { ncp0EqIdx = Indices::ncp0EqIdx };
    enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
    enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };

    typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;

    typedef Opm::MathToolbox<Evaluation> Toolbox;

    typedef Ewoms::EnergyModule<TypeTag, enableEnergy> EnergyModule;
    typedef Ewoms::DiffusionModule<TypeTag, enableDiffusion> DiffusionModule;

public:
    NcpModel(Simulator &simulator)
        : ParentType(simulator)
    {}

    static void registerParameters()
    {
        ParentType::registerParameters();

        DiffusionModule::registerParameters();
        EnergyModule::registerParameters();

        // register runtime parameters of the VTK output modules
        Ewoms::VtkCompositionModule<TypeTag>::registerParameters();

        if (enableDiffusion)
            Ewoms::VtkDiffusionModule<TypeTag>::registerParameters();

        if (enableEnergy)
            Ewoms::VtkEnergyModule<TypeTag>::registerParameters();
    }

    void finishInit()
    {
        ParentType::finishInit();

        minActivityCoeff_.resize(this->numGridDof());
        std::fill(minActivityCoeff_.begin(), minActivityCoeff_.end(), 1.0);
    }

    static std::string name()
    { return "ncp"; }

    std::string primaryVarName(int pvIdx) const
    {
        std::string s;
        if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
            return s;

        std::ostringstream oss;
        if (pvIdx == pressure0Idx)
            oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
        else if (saturation0Idx <= pvIdx && pvIdx < saturation0Idx + (numPhases - 1))
            oss << "saturation_" << FluidSystem::phaseName(/*phaseIdx=*/pvIdx - saturation0Idx);
        else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents)
            oss << "fugacity^" << FluidSystem::componentName(pvIdx - fugacity0Idx);
        else
            assert(false);

        return oss.str();
    }

    std::string eqName(int eqIdx) const
    {
        std::string s;
        if (!(s = EnergyModule::eqName(eqIdx)).empty())
            return s;

        std::ostringstream oss;
        if (conti0EqIdx <= eqIdx && eqIdx < conti0EqIdx + numComponents)
            oss << "continuity^" << FluidSystem::componentName(eqIdx - conti0EqIdx);
        else if (ncp0EqIdx <= eqIdx && eqIdx < ncp0EqIdx + numPhases)
            oss << "ncp_" << FluidSystem::phaseName(/*phaseIdx=*/eqIdx - ncp0EqIdx);
        else
            assert(false);

        return oss.str();
    }

    void updateBegin()
    {
        ParentType::updateBegin();

        // find the a reference pressure. The first degree of freedom
        // might correspond to non-interior entities which would lead
        // to an undefined value, so we have to iterate...
        for (size_t dofIdx = 0; dofIdx < this->numGridDof(); ++ dofIdx) {
            if (this->isLocalDof(dofIdx)) {
                referencePressure_ =
                    this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
                break;
            }
        }
    }

    void updatePVWeights(const ElementContext &elemCtx) const
    {
        for (int dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) {
            int globalIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);

            for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
                minActivityCoeff_[globalIdx][compIdx] = 1e100;
                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
                    const auto &fs = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState();

                    minActivityCoeff_[globalIdx][compIdx] =
                        std::min(minActivityCoeff_[globalIdx][compIdx],
                                 Toolbox::value(fs.fugacityCoefficient(phaseIdx, compIdx))
                                 * Toolbox::value(fs.pressure(phaseIdx)));
                    Valgrind::CheckDefined(minActivityCoeff_[globalIdx][compIdx]);
                }
                if (minActivityCoeff_[globalIdx][compIdx] <= 0)
                    OPM_THROW(Opm::NumericalProblem,
                              "The minumum activity coefficient for component " << compIdx
                              << " on DOF " << globalIdx << " is negative or zero!");
            }
        }
    }

    Scalar primaryVarWeight(int globalDofIdx, int pvIdx) const
    {
        Scalar tmp = EnergyModule::primaryVarWeight(*this, globalDofIdx, pvIdx);
        Scalar result;
        if (tmp > 0)
            // energy related quantity
            result = tmp;
        else if (fugacity0Idx <= pvIdx && pvIdx < fugacity0Idx + numComponents) {
            // component fugacity
            int compIdx = pvIdx - fugacity0Idx;
            assert(0 <= compIdx && compIdx <= numComponents);

            Valgrind::CheckDefined(minActivityCoeff_[globalDofIdx][compIdx]);
            static const Scalar fugacityBaseWeight =
                GET_PROP_VALUE(TypeTag, NcpFugacitiesBaseWeight);
            result = fugacityBaseWeight / minActivityCoeff_[globalDofIdx][compIdx];
        }
        else if (Indices::pressure0Idx == pvIdx) {
            static const Scalar pressureBaseWeight = GET_PROP_VALUE(TypeTag, NcpPressureBaseWeight);
            result = pressureBaseWeight / referencePressure_;
        }
        else {
#ifndef NDEBUG
            int phaseIdx = pvIdx - saturation0Idx;
            assert(0 <= phaseIdx && phaseIdx < numPhases - 1);
#endif

            // saturation
            static const Scalar saturationsBaseWeight =
                GET_PROP_VALUE(TypeTag, NcpSaturationsBaseWeight);
            result = saturationsBaseWeight;
        }

        assert(std::isfinite(result));
        assert(result > 0);

        return result;
    }

    Scalar eqWeight(int globalDofIdx, int eqIdx) const
    {
        Scalar tmp = EnergyModule::eqWeight(*this, globalDofIdx, eqIdx);
        if (tmp > 0)
            // energy related equation
            return tmp;
        // an NCP
        else if (ncp0EqIdx <= eqIdx && eqIdx < Indices::ncp0EqIdx + numPhases)
            return 1.0;

        // mass conservation equation
        int compIdx = eqIdx - Indices::conti0EqIdx;
        assert(0 <= compIdx && compIdx <= numComponents);

        // make all kg equal
        return FluidSystem::molarMass(compIdx);
    }

    Scalar minActivityCoeff(int globalDofIdx, int compIdx) const
    { return minActivityCoeff_[globalDofIdx][compIdx]; }

    void registerOutputModules_()
    {
        ParentType::registerOutputModules_();

        this->addOutputModule(new Ewoms::VtkCompositionModule<TypeTag>(this->simulator_));
        if (enableDiffusion)
            this->addOutputModule(new Ewoms::VtkDiffusionModule<TypeTag>(this->simulator_));
        if (enableEnergy)
            this->addOutputModule(new Ewoms::VtkEnergyModule<TypeTag>(this->simulator_));
    }

    mutable Scalar referencePressure_;
    mutable std::vector<ComponentVector> minActivityCoeff_;
};

} // namespace Ewoms

#endif