From 9247935c8a3e6bd22c4bf87461e021606e94df58 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Mon, 16 Sep 2019 12:07:20 +0200 Subject: [PATCH] changed: ewoms/models/immiscible -> opm/models/immiscible --- examples/co2injection_immiscible_ecfv.cpp | 2 +- examples/co2injection_immiscible_ni_ecfv.cpp | 2 +- examples/co2injection_immiscible_ni_vcfv.cpp | 2 +- examples/co2injection_immiscible_vcfv.cpp | 2 +- examples/finger_immiscible_ecfv.cpp | 2 +- examples/finger_immiscible_vcfv.cpp | 2 +- examples/groundwater_immiscible.cpp | 2 +- examples/lens_immiscible_ecfv_ad.hh | 2 +- examples/lens_immiscible_vcfv_ad.cpp | 2 +- examples/lens_immiscible_vcfv_fd.cpp | 2 +- examples/obstacle_immiscible.cpp | 2 +- examples/powerinjection_darcy_ad.cpp | 2 +- examples/powerinjection_darcy_fd.cpp | 2 +- examples/powerinjection_forchheimer_ad.cpp | 2 +- examples/powerinjection_forchheimer_fd.cpp | 2 +- examples/problems/co2injectionproblem.hh | 2 +- examples/problems/fingerproblem.hh | 2 +- examples/problems/groundwaterproblem.hh | 2 +- examples/problems/lensproblem.hh | 2 +- examples/problems/powerinjectionproblem.hh | 2 +- examples/tutorial1problem.hh | 2 +- .../immiscibleboundaryratevector.hh | 218 +++++++++++ .../immiscibleextensivequantities.hh | 97 +++++ opm/models/immiscible/immiscibleindices.hh | 69 ++++ .../immiscibleintensivequantities.hh | 199 ++++++++++ .../immiscible/immisciblelocalresidual.hh | 199 ++++++++++ opm/models/immiscible/immisciblemodel.hh | 354 ++++++++++++++++++ .../immiscible/immiscibleprimaryvariables.hh | 214 +++++++++++ opm/models/immiscible/immiscibleproperties.hh | 58 +++ opm/models/immiscible/immiscibleratevector.hh | 181 +++++++++ 30 files changed, 1610 insertions(+), 21 deletions(-) create mode 100644 opm/models/immiscible/immiscibleboundaryratevector.hh create mode 100644 opm/models/immiscible/immiscibleextensivequantities.hh create mode 100644 opm/models/immiscible/immiscibleindices.hh create mode 100644 opm/models/immiscible/immiscibleintensivequantities.hh create mode 100644 opm/models/immiscible/immisciblelocalresidual.hh create mode 100644 opm/models/immiscible/immisciblemodel.hh create mode 100644 opm/models/immiscible/immiscibleprimaryvariables.hh create mode 100644 opm/models/immiscible/immiscibleproperties.hh create mode 100644 opm/models/immiscible/immiscibleratevector.hh diff --git a/examples/co2injection_immiscible_ecfv.cpp b/examples/co2injection_immiscible_ecfv.cpp index c04883ab4..80995da23 100644 --- a/examples/co2injection_immiscible_ecfv.cpp +++ b/examples/co2injection_immiscible_ecfv.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_immiscible_ni_ecfv.cpp b/examples/co2injection_immiscible_ni_ecfv.cpp index 0e8ba9f92..8c740e385 100644 --- a/examples/co2injection_immiscible_ni_ecfv.cpp +++ b/examples/co2injection_immiscible_ni_ecfv.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_immiscible_ni_vcfv.cpp b/examples/co2injection_immiscible_ni_vcfv.cpp index cf7ccdba2..04356b102 100644 --- a/examples/co2injection_immiscible_ni_vcfv.cpp +++ b/examples/co2injection_immiscible_ni_vcfv.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/co2injection_immiscible_vcfv.cpp b/examples/co2injection_immiscible_vcfv.cpp index 71f6b81e1..dd2d22714 100644 --- a/examples/co2injection_immiscible_vcfv.cpp +++ b/examples/co2injection_immiscible_vcfv.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include #include "problems/co2injectionproblem.hh" diff --git a/examples/finger_immiscible_ecfv.cpp b/examples/finger_immiscible_ecfv.cpp index 09101557d..404789fe1 100644 --- a/examples/finger_immiscible_ecfv.cpp +++ b/examples/finger_immiscible_ecfv.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include #include "problems/fingerproblem.hh" diff --git a/examples/finger_immiscible_vcfv.cpp b/examples/finger_immiscible_vcfv.cpp index e3fef9247..36a8463f2 100644 --- a/examples/finger_immiscible_vcfv.cpp +++ b/examples/finger_immiscible_vcfv.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include #include "problems/fingerproblem.hh" diff --git a/examples/groundwater_immiscible.cpp b/examples/groundwater_immiscible.cpp index c9a634896..9eb930a0a 100644 --- a/examples/groundwater_immiscible.cpp +++ b/examples/groundwater_immiscible.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include "problems/groundwaterproblem.hh" BEGIN_PROPERTIES diff --git a/examples/lens_immiscible_ecfv_ad.hh b/examples/lens_immiscible_ecfv_ad.hh index cdbadc036..022f5e6c2 100644 --- a/examples/lens_immiscible_ecfv_ad.hh +++ b/examples/lens_immiscible_ecfv_ad.hh @@ -29,7 +29,7 @@ #ifndef EWOMS_LENS_IMMISCIBLE_ECFV_AD_HH #define EWOMS_LENS_IMMISCIBLE_ECFV_AD_HH -#include +#include #include #include "problems/lensproblem.hh" diff --git a/examples/lens_immiscible_vcfv_ad.cpp b/examples/lens_immiscible_vcfv_ad.cpp index 7816ad4c8..a8ac5a0d8 100644 --- a/examples/lens_immiscible_vcfv_ad.cpp +++ b/examples/lens_immiscible_vcfv_ad.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include "problems/lensproblem.hh" BEGIN_PROPERTIES diff --git a/examples/lens_immiscible_vcfv_fd.cpp b/examples/lens_immiscible_vcfv_fd.cpp index efde5b6b0..2d6ef3e9c 100644 --- a/examples/lens_immiscible_vcfv_fd.cpp +++ b/examples/lens_immiscible_vcfv_fd.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include "problems/lensproblem.hh" BEGIN_PROPERTIES diff --git a/examples/obstacle_immiscible.cpp b/examples/obstacle_immiscible.cpp index 426cc6e77..dc2d5eb5a 100644 --- a/examples/obstacle_immiscible.cpp +++ b/examples/obstacle_immiscible.cpp @@ -29,7 +29,7 @@ #include "config.h" #include -#include +#include #include "problems/obstacleproblem.hh" BEGIN_PROPERTIES diff --git a/examples/powerinjection_darcy_ad.cpp b/examples/powerinjection_darcy_ad.cpp index da1b8ecee..3c8fa91de 100644 --- a/examples/powerinjection_darcy_ad.cpp +++ b/examples/powerinjection_darcy_ad.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include "problems/powerinjectionproblem.hh" BEGIN_PROPERTIES diff --git a/examples/powerinjection_darcy_fd.cpp b/examples/powerinjection_darcy_fd.cpp index 5724e9c15..dff65a232 100644 --- a/examples/powerinjection_darcy_fd.cpp +++ b/examples/powerinjection_darcy_fd.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include "problems/powerinjectionproblem.hh" BEGIN_PROPERTIES diff --git a/examples/powerinjection_forchheimer_ad.cpp b/examples/powerinjection_forchheimer_ad.cpp index c514d3ece..e73d19432 100644 --- a/examples/powerinjection_forchheimer_ad.cpp +++ b/examples/powerinjection_forchheimer_ad.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include "problems/powerinjectionproblem.hh" BEGIN_PROPERTIES diff --git a/examples/powerinjection_forchheimer_fd.cpp b/examples/powerinjection_forchheimer_fd.cpp index fcebe1d42..ca8c1555c 100644 --- a/examples/powerinjection_forchheimer_fd.cpp +++ b/examples/powerinjection_forchheimer_fd.cpp @@ -28,7 +28,7 @@ #include "config.h" #include -#include +#include #include "problems/powerinjectionproblem.hh" BEGIN_PROPERTIES diff --git a/examples/problems/co2injectionproblem.hh b/examples/problems/co2injectionproblem.hh index 0a2e26989..1a2b53e1a 100644 --- a/examples/problems/co2injectionproblem.hh +++ b/examples/problems/co2injectionproblem.hh @@ -28,7 +28,7 @@ #ifndef EWOMS_CO2_INJECTION_PROBLEM_HH #define EWOMS_CO2_INJECTION_PROBLEM_HH -#include +#include #include #include diff --git a/examples/problems/fingerproblem.hh b/examples/problems/fingerproblem.hh index f679012dc..c6d405432 100644 --- a/examples/problems/fingerproblem.hh +++ b/examples/problems/fingerproblem.hh @@ -41,7 +41,7 @@ #include #include -#include +#include #include #if HAVE_DUNE_ALUGRID diff --git a/examples/problems/groundwaterproblem.hh b/examples/problems/groundwaterproblem.hh index f833d0d0d..4f62c59ee 100644 --- a/examples/problems/groundwaterproblem.hh +++ b/examples/problems/groundwaterproblem.hh @@ -28,7 +28,7 @@ #ifndef EWOMS_GROUND_WATER_PROBLEM_HH #define EWOMS_GROUND_WATER_PROBLEM_HH -#include +#include #include #include diff --git a/examples/problems/lensproblem.hh b/examples/problems/lensproblem.hh index 17ea6eb4c..c0ca86bd3 100644 --- a/examples/problems/lensproblem.hh +++ b/examples/problems/lensproblem.hh @@ -29,7 +29,7 @@ #define EWOMS_LENS_PROBLEM_HH #include -#include +#include #include #include diff --git a/examples/problems/powerinjectionproblem.hh b/examples/problems/powerinjectionproblem.hh index 18e7bd9aa..b88c06eef 100644 --- a/examples/problems/powerinjectionproblem.hh +++ b/examples/problems/powerinjectionproblem.hh @@ -28,7 +28,7 @@ #ifndef EWOMS_POWER_INJECTION_PROBLEM_HH #define EWOMS_POWER_INJECTION_PROBLEM_HH -#include +#include #include #include diff --git a/examples/tutorial1problem.hh b/examples/tutorial1problem.hh index 294a64fb5..a9fa5a71e 100644 --- a/examples/tutorial1problem.hh +++ b/examples/tutorial1problem.hh @@ -29,7 +29,7 @@ #define EWOMS_TUTORIAL1_PROBLEM_HH /*@\label{tutorial1:guardian2}@*/ // The numerical model -#include +#include // The spatial discretization (VCFV == Vertex-Centered Finite Volumes) #include /*@\label{tutorial1:include-discretization}@*/ diff --git a/opm/models/immiscible/immiscibleboundaryratevector.hh b/opm/models/immiscible/immiscibleboundaryratevector.hh new file mode 100644 index 000000000..c976dfc3e --- /dev/null +++ b/opm/models/immiscible/immiscibleboundaryratevector.hh @@ -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 . + + 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 +#include + +#include "immiscibleintensivequantities.hh" + +namespace Opm { + +/*! + * \ingroup ImmiscibleModel + * + * \brief Implements a boundary vector for the fully implicit + * multi-phase model which assumes immiscibility. + */ +template +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 Toolbox; + typedef Opm::EnergyModule 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 + 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 + 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 + 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 diff --git a/opm/models/immiscible/immiscibleextensivequantities.hh b/opm/models/immiscible/immiscibleextensivequantities.hh new file mode 100644 index 000000000..9de296242 --- /dev/null +++ b/opm/models/immiscible/immiscibleextensivequantities.hh @@ -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 . + + 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 +#include + +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 ImmiscibleExtensiveQuantities + : public MultiPhaseBaseExtensiveQuantities + , public EnergyExtensiveQuantities +{ + typedef MultiPhaseBaseExtensiveQuantities 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 ParameterCache; + typedef Opm::EnergyExtensiveQuantities 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 + 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 diff --git a/opm/models/immiscible/immiscibleindices.hh b/opm/models/immiscible/immiscibleindices.hh new file mode 100644 index 000000000..f1da4f7d9 --- /dev/null +++ b/opm/models/immiscible/immiscibleindices.hh @@ -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 . + + 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 + +namespace Opm { + +/*! + * \ingroup ImmiscibleModel + * + * \brief The indices for the isothermal multi-phase model. + */ +template +struct ImmiscibleIndices + : public EnergyIndices +{ + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) }; + typedef Opm::EnergyIndices 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 diff --git a/opm/models/immiscible/immiscibleintensivequantities.hh b/opm/models/immiscible/immiscibleintensivequantities.hh new file mode 100644 index 000000000..5edd531b3 --- /dev/null +++ b/opm/models/immiscible/immiscibleintensivequantities.hh @@ -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 . + + 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 + +#include +#include + +#include +#include + +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 ImmiscibleIntensiveQuantities + : public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) + , public EnergyIntensiveQuantities + , 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 Toolbox; + typedef Dune::FieldMatrix DimMatrix; + typedef Dune::FieldVector PhaseVector; + typedef Dune::FieldVector EvalPhaseVector; + + typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities; + typedef Opm::EnergyIntensiveQuantities EnergyIntensiveQuantities; + typedef Opm::ImmiscibleFluidState 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 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 diff --git a/opm/models/immiscible/immisciblelocalresidual.hh b/opm/models/immiscible/immisciblelocalresidual.hh new file mode 100644 index 000000000..a209357e8 --- /dev/null +++ b/opm/models/immiscible/immisciblelocalresidual.hh @@ -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 . + + 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 + +#include + +namespace Opm { +/*! + * \ingroup ImmiscibleModel + * + * \brief Calculates the local residual of the immiscible multi-phase + * model. + */ +template +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 EnergyModule; + typedef Opm::MathToolbox 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 + void addPhaseStorage(Dune::FieldVector& 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(intQuants.porosity()) + * Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(fs.density(phaseIdx)); + + EnergyModule::addPhaseStorage(storage, intQuants, phaseIdx); + } + + /*! + * \copydoc FvBaseLocalResidual::computeStorage + */ + template + void computeStorage(Dune::FieldVector& 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(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(this); } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/immiscible/immisciblemodel.hh b/opm/models/immiscible/immisciblemodel.hh new file mode 100644 index 000000000..4193f7ce9 --- /dev/null +++ b/opm/models/immiscible/immisciblemodel.hh @@ -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 . + + 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 +#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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Opm { +template +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); + +//! the Model property +SET_TYPE_PROP(ImmiscibleModel, Model, Opm::ImmiscibleModel); + +//! the RateVector property +SET_TYPE_PROP(ImmiscibleModel, RateVector, Opm::ImmiscibleRateVector); + +//! the BoundaryRateVector property +SET_TYPE_PROP(ImmiscibleModel, BoundaryRateVector, Opm::ImmiscibleBoundaryRateVector); + +//! the PrimaryVariables property +SET_TYPE_PROP(ImmiscibleModel, PrimaryVariables, Opm::ImmisciblePrimaryVariables); + +//! the IntensiveQuantities property +SET_TYPE_PROP(ImmiscibleModel, IntensiveQuantities, Opm::ImmiscibleIntensiveQuantities); + +//! the ExtensiveQuantities property +SET_TYPE_PROP(ImmiscibleModel, ExtensiveQuantities, Opm::ImmiscibleExtensiveQuantities); + +//! The indices required by the isothermal immiscible multi-phase model +SET_TYPE_PROP(ImmiscibleModel, Indices, Opm::ImmiscibleIndices); + +//! 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 type; +}; + +SET_PROP(ImmiscibleSinglePhaseModel, Fluid) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + +public: + typedef Opm::LiquidPhase > 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 > type; +}; + +SET_PROP(ImmiscibleTwoPhaseModel, NonwettingPhase) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + +public: + typedef Opm::LiquidPhase > 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 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); + * \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 ImmiscibleModel + : public Opm::MultiPhaseBaseModel +{ + typedef Opm::MultiPhaseBaseModel 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 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::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(this->simulator_)); + } + +private: + const Implementation& asImp_() const + { return *static_cast(this); } + + mutable Scalar referencePressure_; +}; +} // namespace Opm + +#endif diff --git a/opm/models/immiscible/immiscibleprimaryvariables.hh b/opm/models/immiscible/immiscibleprimaryvariables.hh new file mode 100644 index 000000000..916f3228f --- /dev/null +++ b/opm/models/immiscible/immiscibleprimaryvariables.hh @@ -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 . + + 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 +#include + +#include +#include +#include + +#include + +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 ImmisciblePrimaryVariables : public FvBasePrimaryVariables +{ + typedef FvBasePrimaryVariables 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 Toolbox; + typedef Dune::FieldVector ComponentVector; + typedef Opm::ImmiscibleFlash ImmiscibleFlash; + typedef Opm::EnergyModule 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 + 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 paramCache; + Opm::ImmiscibleFluidState 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(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 + 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(this); } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/immiscible/immiscibleproperties.hh b/opm/models/immiscible/immiscibleproperties.hh new file mode 100644 index 000000000..273d1ca0e --- /dev/null +++ b/opm/models/immiscible/immiscibleproperties.hh @@ -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 . + + 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 +#include + +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 diff --git a/opm/models/immiscible/immiscibleratevector.hh b/opm/models/immiscible/immiscibleratevector.hh new file mode 100644 index 000000000..9e0d907f2 --- /dev/null +++ b/opm/models/immiscible/immiscibleratevector.hh @@ -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 . + + 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 + +#include +#include + +#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 ImmiscibleRateVector + : public Dune::FieldVector +{ + 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 ParentType; + typedef Opm::EnergyModule 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 + 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 + 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 + 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