changed: ewoms/models/immiscible -> opm/models/immiscible

This commit is contained in:
Arne Morten Kvarving 2019-09-16 12:07:20 +02:00
parent 72b5e42016
commit 9247935c8a
30 changed files with 1610 additions and 21 deletions

View File

@ -29,7 +29,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include "problems/co2injectionproblem.hh"

View File

@ -29,7 +29,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include "problems/co2injectionproblem.hh"

View File

@ -29,7 +29,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
#include "problems/co2injectionproblem.hh"

View File

@ -29,7 +29,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
#include "problems/co2injectionproblem.hh"

View File

@ -28,7 +28,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include "problems/fingerproblem.hh"

View File

@ -28,7 +28,7 @@
#include "config.h"
#include <ewoms/common/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/vcfv/vcfvdiscretization.hh>
#include "problems/fingerproblem.hh"

View File

@ -28,7 +28,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/groundwaterproblem.hh"
BEGIN_PROPERTIES

View File

@ -29,7 +29,7 @@
#ifndef EWOMS_LENS_IMMISCIBLE_ECFV_AD_HH
#define EWOMS_LENS_IMMISCIBLE_ECFV_AD_HH
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>
#include "problems/lensproblem.hh"

View File

@ -29,7 +29,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/lensproblem.hh"
BEGIN_PROPERTIES

View File

@ -29,7 +29,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/lensproblem.hh"
BEGIN_PROPERTIES

View File

@ -29,7 +29,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/obstacleproblem.hh"
BEGIN_PROPERTIES

View File

@ -28,7 +28,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/powerinjectionproblem.hh"
BEGIN_PROPERTIES

View File

@ -28,7 +28,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/powerinjectionproblem.hh"
BEGIN_PROPERTIES

View File

@ -28,7 +28,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/powerinjectionproblem.hh"
BEGIN_PROPERTIES

View File

@ -28,7 +28,7 @@
#include "config.h"
#include <opm/models/utils/start.hh>
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include "problems/powerinjectionproblem.hh"
BEGIN_PROPERTIES

View File

@ -28,7 +28,7 @@
#ifndef EWOMS_CO2_INJECTION_PROBLEM_HH
#define EWOMS_CO2_INJECTION_PROBLEM_HH
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <opm/simulators/linalg/parallelamgbackend.hh>
#include <opm/material/fluidsystems/H2ON2FluidSystem.hpp>

View File

@ -41,7 +41,7 @@
#include <opm/material/components/SimpleH2O.hpp>
#include <opm/material/components/Air.hpp>
#include <ewoms/models/immiscible/immiscibleproperties.hh>
#include <opm/models/immiscible/immiscibleproperties.hh>
#include <opm/models/discretization/common/restrictprolong.hh>
#if HAVE_DUNE_ALUGRID

View File

@ -28,7 +28,7 @@
#ifndef EWOMS_GROUND_WATER_PROBLEM_HH
#define EWOMS_GROUND_WATER_PROBLEM_HH
#include <ewoms/models/immiscible/immiscibleproperties.hh>
#include <opm/models/immiscible/immiscibleproperties.hh>
#include <opm/simulators/linalg/parallelistlbackend.hh>
#include <opm/material/components/SimpleH2O.hpp>

View File

@ -29,7 +29,7 @@
#define EWOMS_LENS_PROBLEM_HH
#include <ewoms/io/structuredgridvanguard.hh>
#include <ewoms/models/immiscible/immiscibleproperties.hh>
#include <opm/models/immiscible/immiscibleproperties.hh>
#include <opm/models/discretization/common/fvbaseadlocallinearizer.hh>
#include <opm/models/discretization/ecfv/ecfvdiscretization.hh>

View File

@ -28,7 +28,7 @@
#ifndef EWOMS_POWER_INJECTION_PROBLEM_HH
#define EWOMS_POWER_INJECTION_PROBLEM_HH
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
#include <ewoms/io/cubegridvanguard.hh>
#include <opm/material/fluidmatrixinteractions/RegularizedVanGenuchten.hpp>

View File

@ -29,7 +29,7 @@
#define EWOMS_TUTORIAL1_PROBLEM_HH /*@\label{tutorial1:guardian2}@*/
// The numerical model
#include <ewoms/models/immiscible/immisciblemodel.hh>
#include <opm/models/immiscible/immisciblemodel.hh>
// The spatial discretization (VCFV == Vertex-Centered Finite Volumes)
#include <opm/models/discretization/vcfv/vcfvdiscretization.hh> /*@\label{tutorial1:include-discretization}@*/

View File

@ -0,0 +1,218 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmiscibleBoundaryRateVector
*/
#ifndef EWOMS_IMMISCIBLE_BOUNDARY_RATE_VECTOR_HH
#define EWOMS_IMMISCIBLE_BOUNDARY_RATE_VECTOR_HH
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/constraintsolvers/NcpFlash.hpp>
#include "immiscibleintensivequantities.hh"
namespace Opm {
/*!
* \ingroup ImmiscibleModel
*
* \brief Implements a boundary vector for the fully implicit
* multi-phase model which assumes immiscibility.
*/
template <class TypeTag>
class ImmiscibleBoundaryRateVector : public GET_PROP_TYPE(TypeTag, RateVector)
{
typedef typename GET_PROP_TYPE(TypeTag, RateVector) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
public:
ImmiscibleBoundaryRateVector()
: ParentType()
{}
/*!
* \brief Constructor that assigns all entries to a scalar value.
*
* \param value The scalar value to which all components of the
* boundary rate vector will be set.
*/
ImmiscibleBoundaryRateVector(const Evaluation& value)
: ParentType(value)
{}
/*!
* \brief Copy constructor
*
* \param value The boundary rate vector to be duplicated.
*/
ImmiscibleBoundaryRateVector(const ImmiscibleBoundaryRateVector& value) = default;
ImmiscibleBoundaryRateVector& operator=(const ImmiscibleBoundaryRateVector& value) = default;
/*!
* \brief Specify a free-flow boundary
*
* \param context The execution context for which the boundary rate should
* be specified.
* \param bfIdx The local space index of the boundary segment.
* \param timeIdx The index used by the time discretization.
* \param fluidState The repesentation of the thermodynamic state
* of the system on the integration point of the
* boundary segment.
*/
template <class Context, class FluidState>
void setFreeFlow(const Context& context, unsigned bfIdx, unsigned timeIdx, const FluidState& fluidState)
{
ExtensiveQuantities extQuants;
extQuants.updateBoundary(context, bfIdx, timeIdx, fluidState);
const auto& insideIntQuants = context.intensiveQuantities(bfIdx, timeIdx);
unsigned focusDofIdx = context.focusDofIndex();
unsigned interiorDofIdx = context.interiorScvIndex(bfIdx, timeIdx);
////////
// advective fluxes of all components in all phases
////////
(*this) = Evaluation(0.0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
const auto& pBoundary = fluidState.pressure(phaseIdx);
const Evaluation& pInside = insideIntQuants.fluidState().pressure(phaseIdx);
// mass conservation
Evaluation density;
if (pBoundary > pInside) {
if (focusDofIdx == interiorDofIdx)
density = fluidState.density(phaseIdx);
else
density = Opm::getValue(fluidState.density(phaseIdx));
}
else if (focusDofIdx == interiorDofIdx)
density = insideIntQuants.fluidState().density(phaseIdx);
else
density = Opm::getValue(insideIntQuants.fluidState().density(phaseIdx));
Opm::Valgrind::CheckDefined(density);
Opm::Valgrind::CheckDefined(extQuants.volumeFlux(phaseIdx));
(*this)[conti0EqIdx + phaseIdx] += extQuants.volumeFlux(phaseIdx)*density;
// energy conservation
if (enableEnergy) {
Evaluation specificEnthalpy;
if (pBoundary > pInside) {
if (focusDofIdx == interiorDofIdx)
specificEnthalpy = fluidState.enthalpy(phaseIdx);
else
specificEnthalpy = Opm::getValue(fluidState.enthalpy(phaseIdx));
}
else if (focusDofIdx == interiorDofIdx)
specificEnthalpy = insideIntQuants.fluidState().enthalpy(phaseIdx);
else
specificEnthalpy = Opm::getValue(insideIntQuants.fluidState().enthalpy(phaseIdx));
Evaluation enthalpyRate = density*extQuants.volumeFlux(phaseIdx)*specificEnthalpy;
EnergyModule::addToEnthalpyRate(*this, enthalpyRate);
}
}
// thermal conduction
EnergyModule::addToEnthalpyRate(*this, EnergyModule::thermalConductionRate(extQuants));
#ifndef NDEBUG
for (unsigned i = 0; i < numEq; ++i)
Opm::Valgrind::CheckDefined((*this)[i]);
Opm::Valgrind::CheckDefined(*this);
#endif
}
/*!
* \brief Specify an inflow boundary
*
* An inflow boundary condition is basically a free flow boundary
* condition that is not prevented from specifying a flow out of
* the domain.
*
* \copydetails setFreeFlow
*/
template <class Context, class FluidState>
void setInFlow(const Context& context,
unsigned bfIdx,
unsigned timeIdx,
const FluidState& fluidState)
{
this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
// we only allow fluxes in the direction opposite to the outer unit normal
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
Evaluation& val = this->operator[](eqIdx);
val = Toolbox::min(0.0, val);
}
}
/*!
* \brief Specify an outflow boundary
*
* An outflow boundary condition is basically a free flow boundary
* condition that is not prevented from specifying a flow into
* the domain.
*
* \copydetails setFreeFlow
*/
template <class Context, class FluidState>
void setOutFlow(const Context& context,
unsigned bfIdx,
unsigned timeIdx,
const FluidState& fluidState)
{
this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
// we only allow fluxes in the same direction as the outer unit normal
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
Evaluation& val = this->operator[](eqIdx);
val = Toolbox::max(0.0, val);
}
}
/*!
* \brief Specify a no-flow boundary for all conserved quantities.
*/
void setNoFlow()
{ (*this) = Evaluation(0.0); }
};
} // namespace Opm
#endif

View File

@ -0,0 +1,97 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmiscibleExtensiveQuantities
*/
#ifndef EWOMS_IMMISCIBLE_EXTENSIVE_QUANTITIES_HH
#define EWOMS_IMMISCIBLE_EXTENSIVE_QUANTITIES_HH
#include "immiscibleproperties.hh"
#include <opm/models/common/multiphasebaseextensivequantities.hh>
#include <opm/models/common/energymodule.hh>
namespace Opm {
/*!
* \ingroup ImmiscibleModel
* \ingroup ExtensiveQuantities
*
* \brief This class provides the data all quantities that are required to
* calculate the fluxes of the fluid phases over a face of a
* finite volume for the immiscible multi-phase model.
*
* This means pressure and concentration gradients, phase densities at
* the intergration point, etc.
*/
template <class TypeTag>
class ImmiscibleExtensiveQuantities
: public MultiPhaseBaseExtensiveQuantities<TypeTag>
, public EnergyExtensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy)>
{
typedef MultiPhaseBaseExtensiveQuantities<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef typename FluidSystem::template ParameterCache<Evaluation> ParameterCache;
typedef Opm::EnergyExtensiveQuantities<TypeTag, enableEnergy> EnergyExtensiveQuantities;
public:
/*!
* \copydoc MultiPhaseBaseExtensiveQuantities::registerParameters()
*/
static void registerParameters()
{
ParentType::registerParameters();
}
/*!
* \copydoc MultiPhaseBaseExtensiveQuantities::update()
*/
void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
{
ParentType::update(elemCtx, scvfIdx, timeIdx);
EnergyExtensiveQuantities::update_(elemCtx, scvfIdx, timeIdx);
}
/*!
* \copydoc MultiPhaseBaseExtensiveQuantities::updateBoundary()
*/
template <class Context, class FluidState>
void updateBoundary(const Context& context,
unsigned bfIdx,
unsigned timeIdx,
const FluidState& fluidState)
{
ParentType::updateBoundary(context, bfIdx, timeIdx, fluidState);
EnergyExtensiveQuantities::updateBoundary_(context, bfIdx, timeIdx, fluidState);
}
};
} // namespace Opm
#endif

View File

@ -0,0 +1,69 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmiscibleIndices
*/
#ifndef EWOMS_IMMISCIBLE_INDICES_HH
#define EWOMS_IMMISCIBLE_INDICES_HH
#include "immiscibleproperties.hh"
#include <opm/models/common/energymodule.hh>
namespace Opm {
/*!
* \ingroup ImmiscibleModel
*
* \brief The indices for the isothermal multi-phase model.
*/
template <class TypeTag, int PVOffset>
struct ImmiscibleIndices
: public EnergyIndices<PVOffset + GET_PROP_VALUE(TypeTag, NumPhases),
GET_PROP_VALUE(TypeTag, EnableEnergy)>
{
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef Opm::EnergyIndices<PVOffset + numPhases, enableEnergy> EnergyIndices;
public:
// number of equations/primary variables
static const int numEq = numPhases + EnergyIndices::numEq_;
// Primary variable indices
//! Index for wetting/non-wetting phase pressure
//! (depending on formulation) in a solution vector
static const int pressure0Idx = PVOffset + 0;
//! Index of the saturation of the non-wetting/wetting phase
static const int saturation0Idx = PVOffset + 1;
// indices of the equations
//! Index of the continuity equation of the first phase
static const int conti0EqIdx = PVOffset + 0;
};
} // namespace Opm
#endif

View File

@ -0,0 +1,199 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmiscibleIntensiveQuantities
*/
#ifndef EWOMS_IMMISCIBLE_INTENSIVE_QUANTITIES_HH
#define EWOMS_IMMISCIBLE_INTENSIVE_QUANTITIES_HH
#include "immiscibleproperties.hh"
#include <opm/models/common/energymodule.hh>
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
namespace Opm {
/*!
* \ingroup ImmiscibleModel
* \ingroup IntensiveQuantities
*
* \brief Contains the quantities which are are constant within a
* finite volume for the immiscible multi-phase model.
*/
template <class TypeTag>
class ImmiscibleIntensiveQuantities
: public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
, public EnergyIntensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy)>
, public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
{
typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
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, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { pressure0Idx = Indices::pressure0Idx };
enum { saturation0Idx = Indices::saturation0Idx };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
enum { dimWorld = GridView::dimensionworld };
typedef Opm::MathToolbox<Evaluation> Toolbox;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
typedef Dune::FieldVector<Evaluation, numPhases> EvalPhaseVector;
typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
typedef Opm::EnergyIntensiveQuantities<TypeTag, enableEnergy> EnergyIntensiveQuantities;
typedef Opm::ImmiscibleFluidState<Evaluation, FluidSystem,
/*storeEnthalpy=*/enableEnergy> FluidState;
public:
ImmiscibleIntensiveQuantities()
{ }
ImmiscibleIntensiveQuantities(const ImmiscibleIntensiveQuantities& other) = default;
ImmiscibleIntensiveQuantities& operator=(const ImmiscibleIntensiveQuantities& other) = default;
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
{
ParentType::update(elemCtx, dofIdx, timeIdx);
EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
// material law parameters
const auto& problem = elemCtx.problem();
const typename MaterialLaw::Params& materialParams =
problem.materialLawParams(elemCtx, dofIdx, timeIdx);
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
Opm::Valgrind::CheckDefined(priVars);
Evaluation sumSat = 0.0;
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
const Evaluation& Salpha = priVars.makeEvaluation(saturation0Idx + phaseIdx, timeIdx);
fluidState_.setSaturation(phaseIdx, Salpha);
sumSat += Salpha;
}
fluidState_.setSaturation(numPhases - 1, 1 - sumSat);
EvalPhaseVector pC;
MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
Opm::Valgrind::CheckDefined(pC);
// calculate relative permeabilities
MaterialLaw::relativePermeabilities(relativePermeability_, materialParams, fluidState_);
Opm::Valgrind::CheckDefined(relativePermeability_);
const Evaluation& p0 = priVars.makeEvaluation(pressure0Idx, timeIdx);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
fluidState_.setPressure(phaseIdx, p0 + (pC[phaseIdx] - pC[0]));
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.updateAll(fluidState_);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// compute and set the viscosity
const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
fluidState_.setViscosity(phaseIdx, mu);
// compute and set the density
const Evaluation& rho = FluidSystem::density(fluidState_, paramCache, phaseIdx);
fluidState_.setDensity(phaseIdx, rho);
mobility_[phaseIdx] = relativePermeability_[phaseIdx]/mu;
}
// porosity
porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
// intrinsic permeability
intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
// energy related quantities
EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
// update the quantities specific for the velocity model
FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
}
/*!
* \brief Returns the phase state for the control-volume.
*/
const FluidState& fluidState() const
{ return fluidState_; }
/*!
* \brief Returns the intrinsic permeability tensor a degree of freedom.
*/
const DimMatrix& intrinsicPermeability() const
{ return intrinsicPerm_; }
/*!
* \brief Returns the relative permeability of a given phase
* within the control volume.
*
* \copydetails Doxygen::phaseIdxParam
*/
const Evaluation& relativePermeability(unsigned phaseIdx) const
{ return relativePermeability_[phaseIdx]; }
/*!
* \brief Returns the effective mobility of a given phase within
* the control volume.
*
* \copydetails Doxygen::phaseIdxParam
*/
const Evaluation& mobility(unsigned phaseIdx) const
{ return mobility_[phaseIdx]; }
/*!
* \brief Returns the average porosity within the control volume.
*/
const Evaluation& porosity() const
{ return porosity_; }
protected:
FluidState fluidState_;
Evaluation porosity_;
DimMatrix intrinsicPerm_;
Evaluation relativePermeability_[numPhases];
Evaluation mobility_[numPhases];
};
} // namespace Opm
#endif

View File

@ -0,0 +1,199 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmiscibleLocalResidual
*/
#ifndef EWOMS_IMMISCIBLE_LOCAL_RESIDUAL_BASE_HH
#define EWOMS_IMMISCIBLE_LOCAL_RESIDUAL_BASE_HH
#include "immiscibleproperties.hh"
#include <opm/models/common/energymodule.hh>
#include <opm/material/common/Valgrind.hpp>
namespace Opm {
/*!
* \ingroup ImmiscibleModel
*
* \brief Calculates the local residual of the immiscible multi-phase
* model.
*/
template <class TypeTag>
class ImmiscibleLocalResidual : public GET_PROP_TYPE(TypeTag, DiscLocalResidual)
{
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
typedef Opm::MathToolbox<Evaluation> Toolbox;
public:
/*!
* \brief Adds the amount all conservation quantities (e.g. phase
* mass) within a single fluid phase
*
* \copydetails Doxygen::storageParam
* \copydetails Doxygen::dofCtxParams
* \copydetails Doxygen::phaseIdxParam
*/
template <class LhsEval>
void addPhaseStorage(Dune::FieldVector<LhsEval, numEq>& storage,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx,
unsigned phaseIdx) const
{
// retrieve the intensive quantities for the SCV at the specified
// point in time
const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
const auto& fs = intQuants.fluidState();
storage[conti0EqIdx + phaseIdx] =
Toolbox::template decay<LhsEval>(intQuants.porosity())
* Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
* Toolbox::template decay<LhsEval>(fs.density(phaseIdx));
EnergyModule::addPhaseStorage(storage, intQuants, phaseIdx);
}
/*!
* \copydoc FvBaseLocalResidual::computeStorage
*/
template <class LhsEval>
void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx) const
{
storage = 0.0;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
asImp_().addPhaseStorage(storage, elemCtx, dofIdx, timeIdx, phaseIdx);
EnergyModule::addSolidEnergyStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx));
}
/*!
* \copydoc FvBaseLocalResidual::computeFlux
*/
void computeFlux(RateVector& flux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx) const
{
flux = 0.0;
asImp_().addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
asImp_().addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
}
/*!
* \brief Add the advective mass flux at a given flux integration point
*
* \copydetails computeFlux
*/
void addAdvectiveFlux(RateVector& flux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx) const
{
const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
////////
// advective fluxes of all components in all phases
////////
unsigned focusDofIdx = elemCtx.focusDofIndex();
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
// data attached to upstream DOF of the current phase.
unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, /*timeIdx=*/0);
// add advective flux of current component in current phase.
const Evaluation& rho = up.fluidState().density(phaseIdx);
if (focusDofIdx == upIdx)
flux[conti0EqIdx + phaseIdx] += extQuants.volumeFlux(phaseIdx)*rho;
else
flux[conti0EqIdx + phaseIdx] += extQuants.volumeFlux(phaseIdx)*Toolbox::value(rho);
}
EnergyModule::addAdvectiveFlux(flux, elemCtx, scvfIdx, timeIdx);
}
/*!
* \brief Adds the diffusive flux at a given flux integration point.
*
* For the immiscible model, this is a no-op for mass fluxes. For energy it adds the
* contribution of thermal conduction to the enthalpy flux.
*
* \copydetails computeFlux
*/
void addDiffusiveFlux(RateVector& flux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx) const
{
// no diffusive mass fluxes for the immiscible model
// thermal conduction
EnergyModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
}
/*!
* \copydoc FvBaseLocalResidual::computeSource
*
* By default, this method only asks the problem to specify a
* source term.
*/
void computeSource(RateVector& source,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx) const
{
Opm::Valgrind::SetUndefined(source);
elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
Opm::Valgrind::CheckDefined(source);
}
private:
const Implementation& asImp_() const
{ return *static_cast<const Implementation *>(this); }
};
} // namespace Opm
#endif

View File

@ -0,0 +1,354 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmiscibleModel
*/
#ifndef EWOMS_IMMISCIBLE_MODEL_HH
#define EWOMS_IMMISCIBLE_MODEL_HH
#include <opm/material/densead/Math.hpp>
#include "immiscibleproperties.hh"
#include "immiscibleindices.hh"
#include "immiscibleextensivequantities.hh"
#include "immiscibleprimaryvariables.hh"
#include "immiscibleintensivequantities.hh"
#include "immiscibleratevector.hh"
#include "immiscibleboundaryratevector.hh"
#include "immisciblelocalresidual.hh"
#include <opm/models/common/multiphasebasemodel.hh>
#include <opm/models/common/energymodule.hh>
#include <ewoms/io/vtkenergymodule.hh>
#include <opm/material/components/NullComponent.hpp>
#include <opm/material/fluidsystems/GasPhase.hpp>
#include <opm/material/fluidsystems/LiquidPhase.hpp>
#include <opm/material/fluidsystems/SinglePhaseFluidSystem.hpp>
#include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
#include <sstream>
#include <string>
namespace Opm {
template <class TypeTag>
class ImmiscibleModel;
}
BEGIN_PROPERTIES
//! The generic type tag for problems using the immiscible multi-phase model
NEW_TYPE_TAG(ImmiscibleModel, INHERITS_FROM(MultiPhaseBaseModel, VtkEnergy));
//! The type tag for single-phase immiscible problems
NEW_TYPE_TAG(ImmiscibleSinglePhaseModel, INHERITS_FROM(ImmiscibleModel));
//! The type tag for two-phase immiscible problems
NEW_TYPE_TAG(ImmiscibleTwoPhaseModel, INHERITS_FROM(ImmiscibleModel));
//! Use the immiscible multi-phase local jacobian operator for the immiscible multi-phase model
SET_TYPE_PROP(ImmiscibleModel, LocalResidual,
Opm::ImmiscibleLocalResidual<TypeTag>);
//! the Model property
SET_TYPE_PROP(ImmiscibleModel, Model, Opm::ImmiscibleModel<TypeTag>);
//! the RateVector property
SET_TYPE_PROP(ImmiscibleModel, RateVector, Opm::ImmiscibleRateVector<TypeTag>);
//! the BoundaryRateVector property
SET_TYPE_PROP(ImmiscibleModel, BoundaryRateVector, Opm::ImmiscibleBoundaryRateVector<TypeTag>);
//! the PrimaryVariables property
SET_TYPE_PROP(ImmiscibleModel, PrimaryVariables, Opm::ImmisciblePrimaryVariables<TypeTag>);
//! the IntensiveQuantities property
SET_TYPE_PROP(ImmiscibleModel, IntensiveQuantities, Opm::ImmiscibleIntensiveQuantities<TypeTag>);
//! the ExtensiveQuantities property
SET_TYPE_PROP(ImmiscibleModel, ExtensiveQuantities, Opm::ImmiscibleExtensiveQuantities<TypeTag>);
//! The indices required by the isothermal immiscible multi-phase model
SET_TYPE_PROP(ImmiscibleModel, Indices, Opm::ImmiscibleIndices<TypeTag, /*PVOffset=*/0>);
//! Disable the energy equation by default
SET_BOOL_PROP(ImmiscibleModel, EnableEnergy, false);
/////////////////////
// set slightly different properties for the single-phase case
/////////////////////
//! The fluid system to use by default
SET_PROP(ImmiscibleSinglePhaseModel, FluidSystem)
{ private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid;
public:
typedef Opm::SinglePhaseFluidSystem<Scalar , Fluid> type;
};
SET_PROP(ImmiscibleSinglePhaseModel, Fluid)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Opm::LiquidPhase<Scalar, Opm::NullComponent<Scalar> > type;
};
// disable output of a few quantities which make sense in a
// multi-phase but not in a single-phase context
SET_BOOL_PROP(ImmiscibleSinglePhaseModel, VtkWriteSaturations, false);
SET_BOOL_PROP(ImmiscibleSinglePhaseModel, VtkWriteMobilities, false);
SET_BOOL_PROP(ImmiscibleSinglePhaseModel, VtkWriteRelativePermeabilities, false);
/////////////////////
// set slightly different properties for the two-phase case
/////////////////////
SET_PROP(ImmiscibleTwoPhaseModel, WettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Opm::LiquidPhase<Scalar, Opm::NullComponent<Scalar> > type;
};
SET_PROP(ImmiscibleTwoPhaseModel, NonwettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Opm::LiquidPhase<Scalar, Opm::NullComponent<Scalar> > type;
};
SET_PROP(ImmiscibleTwoPhaseModel, FluidSystem)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, WettingPhase) WettingPhase;
typedef typename GET_PROP_TYPE(TypeTag, NonwettingPhase) NonwettingPhase;
public:
typedef Opm::TwoPhaseImmiscibleFluidSystem<Scalar, WettingPhase, NonwettingPhase> type;
};
END_PROPERTIES
namespace Opm {
/*!
* \ingroup ImmiscibleModel
* \brief A fully-implicit multi-phase flow model which assumes
* immiscibility of the phases.
*
* This model implements multi-phase flow of \f$M > 0\f$ immiscible
* fluids \f$\alpha\f$. By default, the standard multi-phase Darcy
* approach is used to determine the velocity, i.e.
* \f[
* \mathbf{v}_\alpha =
* - \frac{k_{r\alpha}}{\mu_\alpha}
* \mathbf{K}\left(\mathbf{grad}\, p_\alpha -
* \varrho_{\alpha} \mathbf{g} \right) \;,
* \f]
* although the actual approach which is used can be specified via the
* \c FluxModule property. For example, the velocity model can by
* changed to the Forchheimer approach by
* \code
* SET_TYPE_PROP(MyProblemTypeTag, FluxModule, Opm::ForchheimerFluxModule<TypeTag>);
* \endcode
*
* The core of the model is the conservation mass of each component by
* means of the equation
* \f[
* \frac{\partial\;\phi S_\alpha \rho_\alpha }{\partial t}
* - \mathrm{div} \left\{ \rho_\alpha \mathbf{v}_\alpha \right\}
* - q_\alpha = 0 \;.
* \f]
*
* The model uses the following primary variables:
* - The pressure \f$p_0\f$ in Pascal of the phase with the lowest index
* - The saturations \f$S_\alpha\f$ of the \f$M - 1\f$ phases that
* exhibit the lowest indices
* - The absolute temperature \f$T\f$ in Kelvin if energy is conserved
* via the energy equation
*/
template <class TypeTag>
class ImmiscibleModel
: public Opm::MultiPhaseBaseModel<TypeTag>
{
typedef Opm::MultiPhaseBaseModel<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
enum { numComponents = FluidSystem::numComponents };
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
public:
ImmiscibleModel(Simulator& simulator)
: ParentType(simulator)
{}
/*!
* \brief Register all run-time parameters for the immiscible model.
*/
static void registerParameters()
{
ParentType::registerParameters();
if (enableEnergy)
Opm::VtkEnergyModule<TypeTag>::registerParameters();
}
/*!
* \copydoc FvBaseDiscretization::name
*/
static std::string name()
{ return "immiscible"; }
/*!
* \copydoc FvBaseDiscretization::primaryVarName
*/
std::string primaryVarName(unsigned pvIdx) const
{
std::string s;
if (!(s = EnergyModule::primaryVarName(pvIdx)).empty())
return s;
std::ostringstream oss;
if (pvIdx == Indices::pressure0Idx) {
oss << "pressure_" << FluidSystem::phaseName(/*phaseIdx=*/0);
}
else if (Indices::saturation0Idx <= pvIdx
&& pvIdx < Indices::saturation0Idx + numPhases - 1) {
unsigned phaseIdx = pvIdx - Indices::saturation0Idx;
oss << "saturation_" << FluidSystem::phaseName(phaseIdx);
}
else
assert(false);
return oss.str();
}
/*!
* \copydoc FvBaseDiscretization::eqName
*/
std::string eqName(unsigned eqIdx) const
{
std::string s;
if (!(s = EnergyModule::eqName(eqIdx)).empty())
return s;
std::ostringstream oss;
if (Indices::conti0EqIdx <= eqIdx && eqIdx < Indices::conti0EqIdx + numComponents)
oss << "conti_" << FluidSystem::phaseName(eqIdx - Indices::conti0EqIdx);
else
assert(false);
return oss.str();
}
/*!
* \copydoc FvBaseDiscretization::updateBegin
*/
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...
size_t nDof = this->numTotalDof();
for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
if (this->isLocalDof(dofIdx)) {
referencePressure_ =
this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressure0Idx];
break;
}
}
}
/*!
* \copydetails FvBaseDiscretization::primaryVarWeight
*/
Scalar primaryVarWeight(unsigned globalDofIdx, unsigned pvIdx) const
{
assert(referencePressure_ > 0);
Scalar tmp = EnergyModule::primaryVarWeight(asImp_(), globalDofIdx, pvIdx);
if (tmp > 0)
// energy related quantity
return tmp;
if (Indices::pressure0Idx == pvIdx) {
return 10 / referencePressure_;
}
return 1.0;
}
/*!
* \copydetails FvBaseDiscretization::eqWeight
*/
Scalar eqWeight(unsigned globalDofIdx, unsigned eqIdx) const
{
Scalar tmp = EnergyModule::eqWeight(asImp_(), globalDofIdx, eqIdx);
if (tmp > 0)
// energy related equation
return tmp;
#ifndef NDEBUG
unsigned compIdx = eqIdx - Indices::conti0EqIdx;
assert(0 <= compIdx && compIdx <= numPhases);
#endif
// make all kg equal
return 1.0;
}
void registerOutputModules_()
{
ParentType::registerOutputModules_();
if (enableEnergy)
this->addOutputModule(new Opm::VtkEnergyModule<TypeTag>(this->simulator_));
}
private:
const Implementation& asImp_() const
{ return *static_cast<const Implementation *>(this); }
mutable Scalar referencePressure_;
};
} // namespace Opm
#endif

View File

@ -0,0 +1,214 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmisciblePrimaryVariables
*/
#ifndef EWOMS_IMMISCIBLE_PRIMARY_VARIABLES_HH
#define EWOMS_IMMISCIBLE_PRIMARY_VARIABLES_HH
#include "immiscibleproperties.hh"
#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
#include <opm/models/common/energymodule.hh>
#include <opm/material/constraintsolvers/ImmiscibleFlash.hpp>
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <dune/common/fvector.hh>
namespace Opm {
/*!
* \ingroup ImmiscibleModel
*
* \brief Represents the primary variables used by the immiscible
* multi-phase, model.
*
* This class is basically a Dune::FieldVector which can retrieve its
* contents from an aribitatry fluid state.
*/
template <class TypeTag>
class ImmisciblePrimaryVariables : public FvBasePrimaryVariables<TypeTag>
{
typedef FvBasePrimaryVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
// primary variable indices
enum { pressure0Idx = Indices::pressure0Idx };
enum { saturation0Idx = Indices::saturation0Idx };
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
typedef typename Opm::MathToolbox<Evaluation> Toolbox;
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef Opm::ImmiscibleFlash<Scalar, FluidSystem> ImmiscibleFlash;
typedef Opm::EnergyModule<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy)> EnergyModule;
public:
/*!
* \brief Default constructor
*/
ImmisciblePrimaryVariables() : ParentType()
{ Opm::Valgrind::SetUndefined(*this); }
/*!
* \brief Constructor with assignment from scalar
*
* \param value The scalar value to which all entries of the vector will be set.
*/
ImmisciblePrimaryVariables(Scalar value) : ParentType(value)
{}
/*!
* \brief Copy constructor
*
* \param value The primary variables that will be duplicated.
*/
ImmisciblePrimaryVariables(const ImmisciblePrimaryVariables& value) = default;
/*!
* \brief Assignment operator
*
* \param value The primary variables that will be duplicated.
*/
ImmisciblePrimaryVariables& operator=(const ImmisciblePrimaryVariables& value) = default;
/*!
* \brief Set the primary variables from an arbitrary fluid state
* in a mass conservative way.
*
* If an energy equation is included, the fluid temperatures are
* the same as the one given in the fluid state, *not* the
* enthalpy.
*
* \param fluidState The fluid state which should be represented
* by the primary variables. The temperatures,
* pressures, compositions and densities of all
* phases must be defined.
* \param matParams The capillary pressure law parameters
* \param isInEquilibrium If true, the fluid state expresses
* thermodynamic equilibrium assuming the
* relations expressed by the fluid
* system. This implies that in addition to
* the quantities mentioned above, the
* fugacities are also defined.
*/
template <class FluidState>
void assignMassConservative(const FluidState& fluidState,
const MaterialLawParams& matParams,
bool isInEquilibrium = false)
{
#ifndef NDEBUG
// make sure the temperature is the same in all fluid phases
for (unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
assert(std::abs(fluidState.temperature(0) - fluidState.temperature(phaseIdx)) < 1e-30);
}
#endif // NDEBUG
// for the equilibrium case, we don't need complicated
// computations.
if (isInEquilibrium) {
assignNaive(fluidState);
return;
}
// use a flash calculation to calculate a fluid state in
// thermodynamic equilibrium
typename FluidSystem::template ParameterCache<Scalar> paramCache;
Opm::ImmiscibleFluidState<Scalar, FluidSystem> fsFlash;
// use the externally given fluid state as initial value for
// the flash calculation
fsFlash.assign(fluidState);
// calculate the phase densities
paramCache.updateAll(fsFlash);
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
Scalar rho = FluidSystem::density(fsFlash, paramCache, phaseIdx);
fsFlash.setDensity(phaseIdx, rho);
}
// calculate the "global molarities"
ComponentVector globalMolarities(0.0);
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
globalMolarities[compIdx] +=
fsFlash.saturation(phaseIdx) * fsFlash.molarity(phaseIdx, compIdx);
}
}
// run the flash calculation
ImmiscibleFlash::template solve<MaterialLaw>(fsFlash, matParams, paramCache, globalMolarities);
// use the result to assign the primary variables
assignNaive(fsFlash);
}
/*!
* \brief Directly retrieve the primary variables from an
* arbitrary fluid state.
*
* This method retrieves all primary variables from an abitrary
* fluid state without careing whether the state which is
* represented by the resulting primary variables features the
* equivalent mass as the given fluid state. This method is
* massively cheaper and simpler than assignMassConservative() but
* it should be used with care!
*
* \param fluidState The fluid state which should be represented
* by the primary variables. The temperatures,
* pressures, compositions and densities of all
* phases must be defined.
*/
template <class FluidState>
void assignNaive(const FluidState& fluidState)
{
// assign the phase temperatures. this is out-sourced to
// the energy module
EnergyModule::setPriVarTemperatures(asImp_(), fluidState);
(*this)[pressure0Idx] = fluidState.pressure(/*phaseIdx=*/0);
for (unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
(*this)[saturation0Idx + phaseIdx] = fluidState.saturation(phaseIdx);
}
private:
Implementation& asImp_()
{ return *static_cast<Implementation *>(this); }
};
} // namespace Opm
#endif

View File

@ -0,0 +1,58 @@
// -*- 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.
*/
/*!
* \file
* \ingroup ImmiscibleModel
*
* \brief Defines the properties required for the immiscible
* multi-phase model.
*/
#ifndef EWOMS_IMMISCIBLE_PROPERTIES_HH
#define EWOMS_IMMISCIBLE_PROPERTIES_HH
#include <opm/models/common/multiphasebaseproperties.hh>
#include <ewoms/io/vtkenergymodule.hh>
BEGIN_PROPERTIES
//!The fluid systems including the information about the phases
NEW_PROP_TAG(FluidSystem);
//! Specify whether energy should be considered as a conservation quantity or not
NEW_PROP_TAG(EnableEnergy);
// these properties only make sense for the ImmiscibleTwoPhase type tag
//! The wetting phase for two-phase models
NEW_PROP_TAG(WettingPhase);
//! The non-wetting phase for two-phase models
NEW_PROP_TAG(NonwettingPhase);
// these properties only make sense for the ImmiscibleSinglePhase type tag
//! The fluid used by the model
NEW_PROP_TAG(Fluid);
END_PROPERTIES
#endif

View File

@ -0,0 +1,181 @@
// -*- 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.
*/
/*!
* \file
*
* \copydoc Opm::ImmiscibleRateVector
*/
#ifndef EWOMS_IMMISCIBLE_RATE_VECTOR_HH
#define EWOMS_IMMISCIBLE_RATE_VECTOR_HH
#include <dune/common/fvector.hh>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/constraintsolvers/NcpFlash.hpp>
#include "immiscibleintensivequantities.hh"
namespace Opm {
/*!
* \ingroup ImmiscibleModel
*
* \brief Implements a vector representing rates of conserved quantities.
*
* This class is basically a Dune::FieldVector which can be set using
* either mass, molar or volumetric rates.
*/
template <class TypeTag>
class ImmiscibleRateVector
: public Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Evaluation),
GET_PROP_VALUE(TypeTag, NumEq)>
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { conti0EqIdx = Indices::conti0EqIdx };
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
typedef Dune::FieldVector<Evaluation, numEq> ParentType;
typedef Opm::EnergyModule<TypeTag, enableEnergy> EnergyModule;
public:
/*!
* \brief Default constructor
*/
ImmiscibleRateVector() : ParentType()
{ Opm::Valgrind::SetUndefined(*this); }
/*!
* \brief Constructor with assignment from scalar
*
* \param value The scalar value to which all entries of the vector will be set.
*/
ImmiscibleRateVector(const Evaluation& value)
: ParentType(value)
{}
/*!
* \brief Copy constructor
*
* \param value The rate vector that will be duplicated.
*/
ImmiscibleRateVector(const ImmiscibleRateVector& value)
: ParentType(value)
{}
/*!
* \brief Set a mass rate of the conservation quantities.
*
* Enthalpy is _not_ taken into account seperately here. This
* means that it must be set to the desired value in the
* parameter.
*
* \param value The mass rate in \f$[kg/(m^2\,s)]\f$ (unit for areal fluxes)
*/
void setMassRate(const ParentType& value)
{ ParentType::operator=(value); }
/*!
* \brief Set a molar rate of the conservation quantities.
*
* Enthalpy is _not_ taken into account seperately here. This
* means that it must be set to the desired value in the
* parameter.
*
* \param value The new molar rate in \f$[mol/(m^2\,s)]\f$
*/
void setMolarRate(const ParentType& value)
{
// convert to mass rates
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
ParentType::operator[](conti0EqIdx + compIdx) =
value[conti0EqIdx + compIdx]*FluidSystem::molarMass(compIdx);
}
/*!
* \brief Set an enthalpy rate [J/As] where \f$A \in \{m^2, m^3\}\f$
*
* If the energy equation is not enabled, this method is a no-op.
*
* \param rate The enthalpy rate in \f$[J/(m^2\,s)]\f$
*/
template <class RhsEval>
void setEnthalpyRate(const RhsEval& rate)
{ EnergyModule::setEnthalpyRate(*this, rate); }
/*!
* \brief Set a volumetric rate of a phase.
*
* The enthalpy transported into the domain is taken into account
* by this method.
*
* \param fluidState The thermodynamic state of the fluids which
* should be considered. The density and the
* composition of the considered phase must be
* specified before calling this method.
* \param phaseIdx The index of the fluid phase for which the
* given amount of volume should be specified.
* \param volume The volumetric rate of the fluid phase in
* \f$[m^3/(m^2\,s)]\f$ (unit for areal fluxes)
*/
template <class FluidState, class RhsEval>
void setVolumetricRate(const FluidState& fluidState, unsigned phaseIdx, const RhsEval& volume)
{
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
(*this)[conti0EqIdx + compIdx] =
fluidState.density(phaseIdx, compIdx)
* fluidState.massFraction(phaseIdx, compIdx)
* volume;
EnergyModule::setEnthalpyRate(*this, fluidState, phaseIdx, volume);
}
/*!
* \brief Assignment operator from a scalar or a function evaluation
*/
template <class RhsEval>
ImmiscibleRateVector& operator=(const RhsEval& value)
{
for (unsigned i=0; i < this->size(); ++i)
(*this)[i] = value;
return *this;
}
/*!
* \brief Assignment operator from another rate vector
*/
ImmiscibleRateVector& operator=(const ImmiscibleRateVector& other)
{
for (unsigned i=0; i < this->size(); ++i)
(*this)[i] = other[i];
return *this;
}
};
} // namespace Opm
#endif