fingerproblem.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:
/*
  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/>.
 
  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_FINGER_PROBLEM_HH
#define EWOMS_FINGER_PROBLEM_HH
 
#if HAVE_DUNE_ALUGRID
#include <dune/alugrid/grid.hh>
#endif
 
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/version.hh>
 
#include <dune/grid/utility/persistentcontainer.hh>
 
#include <opm/material/components/Air.hpp>
#include <opm/material/components/SimpleH2O.hpp>
 
#include <opm/material/fluidmatrixinteractions/EffToAbsLaw.hpp>
#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
#include <opm/material/fluidmatrixinteractions/ParkerLenhard.hpp>
#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>
 
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
 
#include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
 
#include <opm/models/common/multiphasebaseparameters.hh>
 
#include <opm/models/discretization/common/fvbasefdlocallinearizer.hh>
#include <opm/models/discretization/common/restrictprolong.hh>
 
#include <opm/models/immiscible/immiscibleproperties.hh>
 
#include <opm/models/io/structuredgridvanguard.hh>
 
#include <string>
 
namespace Opm {
template <class TypeTag>
class FingerProblem;
 
} // namespace Opm
 
namespace Opm::Properties {
 
// Create new type tags
namespace TTag {
struct FingerBaseProblem { using InheritsFrom = std::tuple<StructuredGridVanguard>; };
} // end namespace TTag
 
#if HAVE_DUNE_ALUGRID
// use dune-alugrid if available
template<class TypeTag>
struct Grid<TypeTag, TTag::FingerBaseProblem>
{ using type = Dune::ALUGrid</*dim=*/2,
                             /*dimWorld=*/2,
                             Dune::cube,
                             Dune::nonconforming>; };
#endif
 
// Set the problem property
template<class TypeTag>
struct Problem<TypeTag, TTag::FingerBaseProblem> { using type = Opm::FingerProblem<TypeTag>; };
 
// Set the wetting phase
template<class TypeTag>
struct WettingPhase<TypeTag, TTag::FingerBaseProblem>
{
private:
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
public:
    using type = Opm::LiquidPhase<Scalar, Opm::SimpleH2O<Scalar> >;
};
 
// Set the non-wetting phase
template<class TypeTag>
struct NonwettingPhase<TypeTag, TTag::FingerBaseProblem>
{
private:
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
 
public:
    using type = Opm::GasPhase<Scalar, Opm::Air<Scalar> >;
};
 
// Set the material Law
template<class TypeTag>
struct MaterialLaw<TypeTag, TTag::FingerBaseProblem>
{
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using Traits = Opm::TwoPhaseMaterialTraits<Scalar,
                                               /*wettingPhaseIdx=*/FluidSystem::wettingPhaseIdx,
                                               /*nonWettingPhaseIdx=*/FluidSystem::nonWettingPhaseIdx>;
 
    // use the parker-lenhard hysteresis law
    using ParkerLenhard = Opm::ParkerLenhard<Traits>;
    using type = ParkerLenhard;
};
 
// Enable constraints
template<class TypeTag>
struct EnableConstraints<TypeTag, TTag::FingerBaseProblem> { static constexpr int value = true; };
 
} // namespace Opm::Properties
 
namespace Opm::Parameters {
 
template<class Scalar>
struct InitialWaterSaturation { static constexpr Scalar value = 0.01;  };
 
} // namespace Opm::Parameters
 
namespace Opm {
 
template <class TypeTag>
class FingerProblem : public GetPropType<TypeTag, Properties::BaseProblem>
{
    using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
 
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using WettingPhase = GetPropType<TypeTag, Properties::WettingPhase>;
    using NonwettingPhase = GetPropType<TypeTag, Properties::NonwettingPhase>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
    using Constraints = GetPropType<TypeTag, Properties::Constraints>;
    using Model = GetPropType<TypeTag, Properties::Model>;
 
    enum {
        // number of phases
        numPhases = FluidSystem::numPhases,
 
        // phase indices
        wettingPhaseIdx = FluidSystem::wettingPhaseIdx,
        nonWettingPhaseIdx = FluidSystem::nonWettingPhaseIdx,
 
        // equation indices
        contiWettingEqIdx = Indices::conti0EqIdx + wettingPhaseIdx,
 
        // Grid and world dimension
        dim = GridView::dimension,
        dimWorld = GridView::dimensionworld
    };
 
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Stencil = GetPropType<TypeTag, Properties::Stencil> ;
    enum { codim = Stencil::Entity::codimension };
    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
 
    using ParkerLenhard = typename GetProp<TypeTag, Properties::MaterialLaw>::ParkerLenhard;
    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
    using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
 
    using CoordScalar = typename GridView::ctype;
    using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
    using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
 
    using Grid = typename GridView :: Grid;
 
    using MaterialLawParamsContainer = Dune::PersistentContainer< Grid, std::shared_ptr< MaterialLawParams > >  ;
 
public:
    using RestrictProlongOperator = CopyRestrictProlong< Grid, MaterialLawParamsContainer >;
 
    FingerProblem(Simulator& simulator)
        : ParentType(simulator),
          materialParams_( simulator.vanguard().grid(), codim )
    {
    }
 
 
    RestrictProlongOperator restrictProlongOperator()
    {
        return RestrictProlongOperator( materialParams_ );
    }
 
    std::string name() const
    { return
            std::string("finger") +
            "_" + Model::name() +
            "_" + Model::discretizationName() +
            (this->model().enableGridAdaptation()?"_adaptive":"");
    }
 
    static void registerParameters()
    {
        ParentType::registerParameters();
 
        Parameters::Register<Parameters::InitialWaterSaturation<Scalar>>
            ("The initial saturation in the domain [] of the wetting phase");
 
        Parameters::SetDefault<Parameters::CellsX>(20);
        Parameters::SetDefault<Parameters::DomainSizeX<Scalar>>(0.1);
 
        if constexpr (dim > 1) {
            Parameters::SetDefault<Parameters::CellsY>(70);
            Parameters::SetDefault<Parameters::DomainSizeY<Scalar>>(0.3);
        }
        if constexpr (dim == 3) {
            Parameters::SetDefault<Parameters::CellsZ>(1);
            Parameters::SetDefault<Parameters::DomainSizeZ<Scalar>>(0.1);
        }
 
        // Use forward differences
        Parameters::SetDefault<Parameters::NumericDifferenceMethod>(+1);
 
        Parameters::SetDefault<Parameters::EndTime<Scalar>>(215);
        Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(10);
        Parameters::SetDefault<Parameters::EnableGravity>(true);
    }
 
    void finishInit()
    {
        ParentType::finishInit();
 
        eps_ = 3e-6;
 
        temperature_ = 273.15 + 20; // -> 20°C
 
        FluidSystem::init();
 
        // parameters for the Van Genuchten law of the main imbibition
        // and the main drainage curves.
        micParams_.setVgAlpha(0.0037);
        micParams_.setVgN(4.7);
        micParams_.finalize();
 
        mdcParams_.setVgAlpha(0.0037);
        mdcParams_.setVgN(4.7);
        mdcParams_.finalize();
 
        // initialize the material parameter objects of the individual
        // finite volumes, resize will resize the container to the number of elements
        materialParams_.resize();
 
        for (auto it = materialParams_.begin(),
                 end = materialParams_.end(); it != end; ++it ) {
            std::shared_ptr< MaterialLawParams >& materialParams = *it ;
            if( ! materialParams )
            {
                materialParams.reset( new MaterialLawParams() );
                materialParams->setMicParams(&micParams_);
                materialParams->setMdcParams(&mdcParams_);
                materialParams->setSwr(0.0);
                materialParams->setSnr(0.1);
                materialParams->finalize();
                ParkerLenhard::reset(*materialParams);
            }
        }
 
        K_ = this->toDimMatrix_(4.6e-10);
 
        setupInitialFluidState_();
    }
 
    void endTimeStep()
    {
#ifndef NDEBUG
        // checkConservativeness() does not include the effect of constraints, so we
        // disable it for this problem...
        //this->model().checkConservativeness();
 
        // Calculate storage terms
        EqVector storage;
        this->model().globalStorage(storage);
 
        // Write mass balance information for rank 0
        if (this->gridView().comm().rank() == 0) {
            std::cout << "Storage: " << storage << std::endl << std::flush;
        }
#endif // NDEBUG
 
        // update the history of the hysteresis law
        ElementContext elemCtx(this->simulator());
 
        for (const auto& elem : elements(this->gridView())) {
            elemCtx.updateAll(elem);
            size_t numDofs = elemCtx.numDof(/*timeIdx=*/0);
            for (unsigned scvIdx = 0; scvIdx < numDofs; ++scvIdx)
            {
                MaterialLawParams& materialParam = materialLawParams( elemCtx, scvIdx, /*timeIdx=*/0 );
                const auto& fs = elemCtx.intensiveQuantities(scvIdx, /*timeIdx=*/0).fluidState();
                ParkerLenhard::update(materialParam, fs);
            }
        }
    }
 
 
 
    template <class Context>
    Scalar temperature(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
    { return temperature_; }
 
    template <class Context>
    const DimMatrix& intrinsicPermeability(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
    { return K_; }
 
    template <class Context>
    Scalar porosity(const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
    { return 0.4; }
 
    template <class Context>
    MaterialLawParams& materialLawParams(const Context& context,
                                         unsigned spaceIdx, unsigned timeIdx)
    {
        const auto& entity = context.stencil(timeIdx).entity(spaceIdx);
        assert(materialParams_[entity]);
        return *materialParams_[entity];
    }
 
    template <class Context>
    const MaterialLawParams& materialLawParams(const Context& context,
                                               unsigned spaceIdx, unsigned timeIdx) const
    {
        const auto& entity = context.stencil(timeIdx).entity( spaceIdx );
        assert(materialParams_[entity]);
        return *materialParams_[entity];
    }
 
 
 
    template <class Context>
    void boundary(BoundaryRateVector& values, const Context& context,
                  unsigned spaceIdx, unsigned timeIdx) const
    {
        const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
 
        if (onLeftBoundary_(pos) || onRightBoundary_(pos) || onLowerBoundary_(pos))
            values.setNoFlow();
        else {
            assert(onUpperBoundary_(pos));
 
            values.setFreeFlow(context, spaceIdx, timeIdx, initialFluidState_);
        }
 
        // override the value for the liquid phase by forced
        // imbibition of water on inlet boundary segments
        if (onInlet_(pos)) {
            values[contiWettingEqIdx] = -0.001; // [kg/(m^2 s)]
        }
    }
 
 
 
    template <class Context>
    void initial(PrimaryVariables& values, const Context& /*context*/, unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
    {
        // assign the primary variables
        values.assignNaive(initialFluidState_);
    }
 
    template <class Context>
    void constraints(Constraints& constraints, const Context& context,
                     unsigned spaceIdx, unsigned timeIdx) const
    {
        const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
 
        if (onUpperBoundary_(pos) && !onInlet_(pos)) {
            constraints.setActive(true);
            constraints.assignNaive(initialFluidState_);
        }
        else if (onLowerBoundary_(pos)) {
            constraints.setActive(true);
            constraints.assignNaive(initialFluidState_);
        }
    }
 
    template <class Context>
    void source(RateVector& rate, const Context& /*context*/,
                unsigned /*spaceIdx*/, unsigned /*timeIdx*/) const
    { rate = Scalar(0.0); }
 
private:
    bool onLeftBoundary_(const GlobalPosition& pos) const
    { return pos[0] < this->boundingBoxMin()[0] + eps_; }
 
    bool onRightBoundary_(const GlobalPosition& pos) const
    { return pos[0] > this->boundingBoxMax()[0] - eps_; }
 
    bool onLowerBoundary_(const GlobalPosition& pos) const
    { return pos[1] < this->boundingBoxMin()[1] + eps_; }
 
    bool onUpperBoundary_(const GlobalPosition& pos) const
    { return pos[1] > this->boundingBoxMax()[1] - eps_; }
 
    bool onInlet_(const GlobalPosition& pos) const
    {
        Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
        Scalar lambda = (this->boundingBoxMax()[0] - pos[0]) / width;
 
        if (!onUpperBoundary_(pos))
            return false;
 
        Scalar xInject[] = { 0.25, 0.75 };
        Scalar injectLen[] = { 0.1, 0.1 };
        for (unsigned i = 0; i < sizeof(xInject) / sizeof(Scalar); ++i) {
            if (xInject[i] - injectLen[i] / 2 < lambda
                && lambda < xInject[i] + injectLen[i] / 2)
                return true;
        }
        return false;
    }
 
    void setupInitialFluidState_()
    {
        auto& fs = initialFluidState_;
        fs.setPressure(wettingPhaseIdx, /*pressure=*/1e5);
 
        Scalar Sw = Parameters::Get<Parameters::InitialWaterSaturation<Scalar>>();
        fs.setSaturation(wettingPhaseIdx, Sw);
        fs.setSaturation(nonWettingPhaseIdx, 1 - Sw);
 
        fs.setTemperature(temperature_);
 
        // set the absolute pressures
        Scalar pn = 1e5;
        fs.setPressure(nonWettingPhaseIdx, pn);
        fs.setPressure(wettingPhaseIdx, pn);
 
        typename FluidSystem::template ParameterCache<Scalar> paramCache;
        paramCache.updateAll(fs);
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
            fs.setDensity(phaseIdx, FluidSystem::density(fs, paramCache, phaseIdx));
            fs.setViscosity(phaseIdx, FluidSystem::viscosity(fs, paramCache, phaseIdx));
        }
 
    }
 
    DimMatrix K_;
 
    typename MaterialLawParams::VanGenuchtenParams micParams_;
    typename MaterialLawParams::VanGenuchtenParams mdcParams_;
 
    MaterialLawParamsContainer materialParams_;
 
    Opm::ImmiscibleFluidState<Scalar, FluidSystem> initialFluidState_;
 
    Scalar temperature_;
    Scalar eps_;
};
 
} // namespace Opm
 
#endif