test_mpfa2pproblem.hh

Go to the documentation of this file.

/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program 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 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
 
#ifndef DUMUX_TEST_MPFA2P_PROBLEM_HH
#define DUMUX_TEST_MPFA2P_PROBLEM_HH
 
#if HAVE_UG
#include <dune/grid/uggrid.hh>
#endif
#include <dune/grid/yaspgrid.hh>
 
#include <dumux/material/fluidsystems/2pimmiscible.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/simpleh2o.hh>
#include <dumux/material/components/trichloroethene.hh>
 
#include <dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressureproperties.hh>
#include <dumux/porousmediumflow/2p/sequential/diffusion/cellcentered/pressurepropertiesadaptive.hh>
#include <dumux/porousmediumflow/2p/sequential/diffusion/mpfa/omethod/2dpressureproperties.hh>
#include <dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressureproperties.hh>
#include <dumux/porousmediumflow/2p/sequential/diffusion/mpfa/lmethod/2dpressurepropertiesadaptive.hh>
#include <dumux/porousmediumflow/2p/sequential/transport/cellcentered/properties.hh>
#include <dumux/porousmediumflow/2p/sequential/impes/problem.hh>
 
 
#include <dumux/porousmediumflow/2p/sequential/transport/cellcentered/evalcflfluxcoats.hh>
#include <dumux/porousmediumflow/2p/sequential/impes/gridadaptionindicatorlocal.hh>
 
#include "test_mpfa2pspatialparams.hh"
#include "buckleyleverettanalyticsolution.hh"
#include "mcwhorteranalyticsolution.hh"
 
namespace Dumux
{
 
template<class TypeTag>
class MPFATwoPTestProblem;
 
// Specify the properties
namespace Properties
{
 
NEW_TYPE_TAG(MPFATwoPTest, INHERITS_FROM(Test2PSpatialParams));
 
// Set the grid type
#if HAVE_UG
SET_TYPE_PROP(MPFATwoPTest, Grid, Dune::UGGrid<2>);
#else
SET_TYPE_PROP(MPFATwoPTest, Grid, Dune::YaspGrid<2>);
#endif
 
// Set the problem property
SET_TYPE_PROP(MPFATwoPTest, Problem, MPFATwoPTestProblem<TypeTag>);
 
// Set the fluid system
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::MPFATwoPTest>
{
    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
    using WettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
#if PROBLEM == 2
    using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::Trichloroethene<Scalar> >;
#else
    using NonwettingPhase = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
#endif
    using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
};
 
#if PROBLEM == 1
SET_INT_PROP(MPFATwoPTest, Formulation, SequentialTwoPCommonIndices::pnsw);
#endif
 
SET_TYPE_PROP(MPFATwoPTest, EvalCflFluxFunction, EvalCflFluxCoats<TypeTag>);
SET_TYPE_PROP(MPFATwoPTest, AdaptionIndicator, GridAdaptionIndicator2PLocal<TypeTag>);
 
NEW_TYPE_TAG(FVTwoPTest, INHERITS_FROM(FVPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTest));
NEW_TYPE_TAG(FVAdaptiveTwoPTest, INHERITS_FROM(FVPressureTwoPAdaptive, FVTransportTwoP, IMPESTwoPAdaptive, MPFATwoPTest));
NEW_TYPE_TAG(MPFAOTwoPTest, INHERITS_FROM(FvMpfaO2dPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTest));
NEW_TYPE_TAG(MPFALTwoPTest, INHERITS_FROM(FvMpfaL2dPressureTwoP, FVTransportTwoP, IMPESTwoP, MPFATwoPTest));
NEW_TYPE_TAG(MPFALAdaptiveTwoPTest, INHERITS_FROM(FvMpfaL2dPressureTwoPAdaptive, FVTransportTwoP, IMPESTwoPAdaptive, MPFATwoPTest));
 
}
 
template<class TypeTag>
class MPFATwoPTestProblem: public IMPESProblem2P<TypeTag>
{
using ParentType = IMPESProblem2P<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Grid = typename GridView::Grid;
 
using Indices = typename GET_PROP_TYPE(TypeTag, ModelTraits)::Indices;
 
using WettingPhase = typename GET_PROP(TypeTag, FluidSystem)::WettingPhase;
 
using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
 
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using SolutionTypes = typename GET_PROP(TypeTag, SolutionTypes);
using PrimaryVariables = typename SolutionTypes::PrimaryVariables;
 
enum
{
    dim = GridView::dimension, dimWorld = GridView::dimensionworld
};
 
enum
{
    nPhaseIdx = Indices::nPhaseIdx,
#if PROBLEM == 1
    pnIdx = Indices::pnIdx,
#else
    pwIdx = Indices::pwIdx,
#endif
    swIdx = Indices::swIdx,
    eqIdxPress = Indices::pressureEqIdx,
    eqIdxSat = Indices::satEqIdx
};
 
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
 
using Element = typename GridView::Traits::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
 
public:
MPFATwoPTestProblem(TimeManager &timeManager, Grid &grid) :
ParentType(timeManager, grid)
#if PROBLEM != 2
, analyticSolution_(*this)
#endif
{
    int refinementFactor = getParam<Scalar>("Grid.RefinementFactor", 0);
    this->grid().globalRefine(refinementFactor);
 
    Scalar inletWidth = getParam<Scalar>("Problem.InletWidth", 1.0);
    GlobalPosition inletCenter = this->bBoxMax();
    inletCenter[0] *= 0.5;
 
    inletLeftCoord_ = inletCenter;
    inletLeftCoord_[0] -=0.5*inletWidth;
    inletRightCoord_ = inletCenter;
    inletRightCoord_[0] +=0.5*inletWidth;
 
    inFlux_ = getParam<Scalar>("Problem.InjectionFlux", 1e-4);
 
    int outputInterval = getParam<int>("Problem.OutputInterval", 0);
    this->setOutputInterval(outputInterval);
 
    Scalar outputTimeInterval = getParam<Scalar>("Problem.OutputTimeInterval", 1e6);
    this->setOutputTimeInterval(outputTimeInterval);
}
 
#if PROBLEM != 2
void init()
{
    ParentType::init();
 
#if PROBLEM == 0
    Scalar vTot = 3e-6;
    analyticSolution_.initialize(vTot);
#endif
#if PROBLEM == 1
    analyticSolution_.initialize();
#endif
}
 
void addOutputVtkFields()
{
    analyticSolution_.calculateAnalyticSolution();
 
    ParentType::addOutputVtkFields();
    analyticSolution_.addOutputVtkFields(this->resultWriter());
}
#endif
 
// \{
 
std::string name() const
{
    return getParam<std::string>("Problem.Name");
}
 
bool shouldWriteRestartFile() const
{
    return false;
}
 
 
Scalar temperatureAtPos(const GlobalPosition& globalPos) const
{
    return 273.15 + 10; // -> 10°C
}
 
// \}
 
Scalar referencePressureAtPos(const GlobalPosition& globalPos) const
{
    return 1e5; // -> 10°C
}
 
void source(PrimaryVariables &values,const Element& element) const
{
    values = 0;
}
 
void boundaryTypesAtPos(BoundaryTypes &bcTypes, const GlobalPosition& globalPos) const
{
#if PROBLEM == 2
    if (isInlet(globalPos))
    {
        bcTypes.setNeumann(eqIdxPress);
        bcTypes.setDirichlet(swIdx);
    }
    else if (isBottom(globalPos) || isTop(globalPos))
    {
        bcTypes.setAllNeumann();
    }
    else
    {
        bcTypes.setDirichlet(pwIdx);
        bcTypes.setOutflow(eqIdxSat);
    }
#elif  PROBLEM == 0
    if (globalPos[0] < eps_)
    {
        bcTypes.setAllDirichlet();
    }
    else if (globalPos[0] > this->bBoxMax()[0] - eps_)
    {
        bcTypes.setNeumann(eqIdxPress);
        bcTypes.setOutflow(eqIdxSat);
    }
    else
    {
        bcTypes.setAllNeumann();
    }
#else
    if (globalPos[0] < eps_)
    {
        bcTypes.setAllDirichlet();
    }
    else
    {
        bcTypes.setAllNeumann();
    }
#endif
}
 
void dirichletAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
    values = 0;
#if PROBLEM == 2
    Scalar pRef = referencePressureAtPos(globalPos);
    Scalar temp = temperatureAtPos(globalPos);
 
    values[pwIdx] = (1e5 - (this->bBoxMax()- globalPos) * this->gravity() * WettingPhase::density(temp, pRef));
    values[swIdx] = 1.0;
 
    if (isInlet(globalPos))
    {
        values[swIdx] = 0.0;
    }
#elif PROBLEM == 0
    if (globalPos[0] < eps_)
    {
        values[swIdx] = 0.8;
        values[pwIdx] = 1;
    }
#else
    if (globalPos[0] < eps_)
    {
        values[swIdx] = 1.0;
        values[pnIdx] = 1e5;
    }
#endif
}
 
void neumannAtPos(PrimaryVariables &values, const GlobalPosition& globalPos) const
{
    values = 0;
#if PROBLEM == 2
    if (isInlet(globalPos))
    {
        values[nPhaseIdx] = -inFlux_;
    }
#elif PROBLEM == 0
    if (globalPos[0] > this->bBoxMax()[0] - eps_)
    {
        values[nPhaseIdx] = 3e-3;
    }
#endif
}
 
void initialAtPos(PrimaryVariables &values,
        const GlobalPosition& globalPos) const
{
#if PROBLEM == 2
    values[pwIdx] = 0;
    values[swIdx] = 1.0;
#elif PROBLEM == 0
    values[pwIdx] = 0;
    values[swIdx] = 0.2;
#else
    values[pnIdx] = 0;
    values[swIdx] = 0.0;
#endif
}
 
private:
 
bool isInlet(const GlobalPosition& globalPos) const
{
        if (!isTop(globalPos))
            return false;
 
    for (int i = 0; i < dim; i++)
    {
        if (globalPos[i] < inletLeftCoord_[i] - eps_)
            return false;
        if (globalPos[i] > inletRightCoord_[i] + eps_)
            return false;
    }
    return true;
}
 
bool isTop(const GlobalPosition& globalPos) const
{
    if (dim == 2)
    {
        if (globalPos[1] > this->bBoxMax()[1] - eps_)
            return true;
    }
    if (dim == 3)
    {
        if (globalPos[2] > this->bBoxMax()[2] - eps_)
            return true;
    }
    return false;
}
 
bool isBottom(const GlobalPosition& globalPos) const
{
    if (dim == 2)
    {
        if (globalPos[1] < eps_)
            return true;
    }
    if (dim == 3)
    {
        if (globalPos[2] < eps_)
            return true;
    }
    return false;
}
 
Scalar inFlux_;
GlobalPosition inletLeftCoord_;
GlobalPosition inletRightCoord_;
static constexpr Scalar eps_ = 1e-6;
#if PROBLEM == 0
BuckleyLeverettAnalytic<TypeTag> analyticSolution_;
#endif
#if PROBLEM == 1
McWhorterAnalytic<TypeTag> analyticSolution_;
#endif
};
} //end namespace
 
#endif // DUMUX_TEST_MPFA2P_PROBLEM_HH