diff --git a/opm/models/common/darcyfluxmodule.hh b/opm/models/common/darcyfluxmodule.hh new file mode 100644 index 000000000..410cbbed5 --- /dev/null +++ b/opm/models/common/darcyfluxmodule.hh @@ -0,0 +1,569 @@ +// -*- 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 + * + * \brief This file contains the necessary classes to calculate the + * volumetric fluxes out of a pressure potential gradient using the + * Darcy relation. + */ +#ifndef EWOMS_DARCY_FLUX_MODULE_HH +#define EWOMS_DARCY_FLUX_MODULE_HH + +#include "multiphasebaseproperties.hh" +#include + +#include +#include +#include + +#include +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(MaterialLaw); + +END_PROPERTIES + +namespace Opm { + +template +class DarcyIntensiveQuantities; + +template +class DarcyExtensiveQuantities; + +template +class DarcyBaseProblem; + +/*! + * \ingroup FluxModules + * \brief Specifies a flux module which uses the Darcy relation. + */ +template +struct DarcyFluxModule +{ + typedef DarcyIntensiveQuantities FluxIntensiveQuantities; + typedef DarcyExtensiveQuantities FluxExtensiveQuantities; + typedef DarcyBaseProblem FluxBaseProblem; + + /*! + * \brief Register all run-time parameters for the flux module. + */ + static void registerParameters() + { } +}; + +/*! + * \ingroup FluxModules + * \brief Provides the defaults for the parameters required by the + * Darcy velocity approach. + */ +template +class DarcyBaseProblem +{ }; + +/*! + * \ingroup FluxModules + * \brief Provides the intensive quantities for the Darcy flux module + */ +template +class DarcyIntensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; +protected: + void update_(const ElementContext& elemCtx OPM_UNUSED, + unsigned dofIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + { } +}; + +/*! + * \ingroup FluxModules + * \brief Provides the Darcy flux module + * + * The commonly used Darcy relation looses its validity for Reynolds numbers \f$ Re < + * 1\f$. If one encounters flow velocities in porous media above this threshold, the + * Forchheimer approach can be used. + * + * The Darcy equation is given by the following relation: + * + * \f[ + \vec{v}_\alpha = + \left( \nabla p_\alpha - \rho_\alpha \vec{g}\right) + \frac{\mu_\alpha}{k_{r,\alpha} K} + \f] + */ +template +class DarcyExtensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw; + + enum { dimWorld = GridView::dimensionworld }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + + typedef typename Opm::MathToolbox Toolbox; + typedef typename FluidSystem::template ParameterCache ParameterCache; + typedef Dune::FieldVector EvalDimVector; + typedef Dune::FieldVector DimVector; + typedef Dune::FieldMatrix DimMatrix; + +public: + /*! + * \brief Returns the intrinsic permeability tensor for a given + * sub-control volume face. + */ + const DimMatrix& intrinsicPermability() const + { return K_; } + + /*! + * \brief Return the pressure potential gradient of a fluid phase + * at the face's integration point [Pa/m] + * + * \param phaseIdx The index of the fluid phase + */ + const EvalDimVector& potentialGrad(unsigned phaseIdx) const + { return potentialGrad_[phaseIdx]; } + + /*! + * \brief Return the filter velocity of a fluid phase at the + * face's integration point [m/s] + * + * \param phaseIdx The index of the fluid phase + */ + const EvalDimVector& filterVelocity(unsigned phaseIdx) const + { return filterVelocity_[phaseIdx]; } + + /*! + * \brief Return the volume flux of a fluid phase at the face's integration point + * \f$[m^3/s / m^2]\f$ + * + * This is the fluid volume of a phase per second and per square meter of face + * area. + * + * \param phaseIdx The index of the fluid phase + */ + const Evaluation& volumeFlux(unsigned phaseIdx) const + { return volumeFlux_[phaseIdx]; } + +protected: + short upstreamIndex_(unsigned phaseIdx) const + { return upstreamDofIdx_[phaseIdx]; } + + short downstreamIndex_(unsigned phaseIdx) const + { return downstreamDofIdx_[phaseIdx]; } + + /*! + * \brief Calculate the gradients which are required to determine the volumetric fluxes + * + * The the upwind directions is also determined by method. + */ + void calculateGradients_(const ElementContext& elemCtx, + unsigned faceIdx, + unsigned timeIdx) + { + const auto& gradCalc = elemCtx.gradientCalculator(); + Opm::PressureCallback pressureCallback(elemCtx); + + const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(faceIdx); + const auto& faceNormal = scvf.normal(); + + unsigned i = scvf.interiorIndex(); + unsigned j = scvf.exteriorIndex(); + interiorDofIdx_ = static_cast(i); + exteriorDofIdx_ = static_cast(j); + unsigned focusDofIdx = elemCtx.focusDofIndex(); + + // calculate the "raw" pressure gradient + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) { + Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]); + continue; + } + + pressureCallback.setPhaseIndex(phaseIdx); + gradCalc.calculateGradient(potentialGrad_[phaseIdx], + elemCtx, + faceIdx, + pressureCallback); + Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]); + } + + // correct the pressure gradients by the gravitational acceleration + if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) { + // estimate the gravitational acceleration at a given SCV face + // using the arithmetic mean + const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx); + const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx); + + const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx); + const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx); + + const auto& posIn = elemCtx.pos(i, timeIdx); + const auto& posEx = elemCtx.pos(j, timeIdx); + const auto& posFace = scvf.integrationPos(); + + // the distance between the centers of the control volumes + DimVector distVecIn(posIn); + DimVector distVecEx(posEx); + DimVector distVecTotal(posEx); + + distVecIn -= posFace; + distVecEx -= posFace; + distVecTotal -= posIn; + Scalar absDistTotalSquared = distVecTotal.two_norm2(); + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + // calculate the hydrostatic pressure at the integration point of the face + Evaluation pStatIn; + + if (std::is_same::value || + interiorDofIdx_ == static_cast(focusDofIdx)) + { + const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx); + pStatIn = - rhoIn*(gIn*distVecIn); + } + else { + Scalar rhoIn = Toolbox::value(intQuantsIn.fluidState().density(phaseIdx)); + pStatIn = - rhoIn*(gIn*distVecIn); + } + + // the quantities on the exterior side of the face do not influence the + // result for the TPFA scheme, so they can be treated as scalar values. + Evaluation pStatEx; + + if (std::is_same::value || + exteriorDofIdx_ == static_cast(focusDofIdx)) + { + const Evaluation& rhoEx = intQuantsEx.fluidState().density(phaseIdx); + pStatEx = - rhoEx*(gEx*distVecEx); + } + else { + Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx)); + pStatEx = - rhoEx*(gEx*distVecEx); + } + + // compute the hydrostatic gradient between the two control volumes (this + // gradient exhibitis the same direction as the vector between the two + // control volume centers and the length (pStaticExterior - + // pStaticInterior)/distanceInteriorToExterior + Dune::FieldVector f(distVecTotal); + f *= (pStatEx - pStatIn)/absDistTotalSquared; + + // calculate the final potential gradient + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + potentialGrad_[phaseIdx][dimIdx] += f[dimIdx]; + + for (unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) { + if (!Opm::isfinite(potentialGrad_[phaseIdx][dimIdx])) { + throw Opm::NumericalIssue("Non-finite potential gradient for phase '" + +std::string(FluidSystem::phaseName(phaseIdx))+"'"); + } + } + } + } + + Opm::Valgrind::SetUndefined(K_); + elemCtx.problem().intersectionIntrinsicPermeability(K_, elemCtx, faceIdx, timeIdx); + Opm::Valgrind::CheckDefined(K_); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) { + Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]); + continue; + } + + // determine the upstream and downstream DOFs + Evaluation tmp = 0.0; + for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx) + tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx]; + + if (tmp > 0) { + upstreamDofIdx_[phaseIdx] = exteriorDofIdx_; + downstreamDofIdx_[phaseIdx] = interiorDofIdx_; + } + else { + upstreamDofIdx_[phaseIdx] = interiorDofIdx_; + downstreamDofIdx_[phaseIdx] = exteriorDofIdx_; + } + + // we only carry the derivatives along if the upstream DOF is the one which + // we currently focus on + const auto& up = elemCtx.intensiveQuantities(upstreamDofIdx_[phaseIdx], timeIdx); + if (upstreamDofIdx_[phaseIdx] == static_cast(focusDofIdx)) + mobility_[phaseIdx] = up.mobility(phaseIdx); + else + mobility_[phaseIdx] = Toolbox::value(up.mobility(phaseIdx)); + } + } + + /*! + * \brief Calculate the gradients at the grid boundary which are required to + * determine the volumetric fluxes + * + * The the upwind directions is also determined by method. + */ + template + void calculateBoundaryGradients_(const ElementContext& elemCtx, + unsigned boundaryFaceIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + const auto& gradCalc = elemCtx.gradientCalculator(); + Opm::BoundaryPressureCallback pressureCallback(elemCtx, fluidState); + + // calculate the pressure gradient + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) { + Opm::Valgrind::SetUndefined(potentialGrad_[phaseIdx]); + continue; + } + + pressureCallback.setPhaseIndex(phaseIdx); + gradCalc.calculateBoundaryGradient(potentialGrad_[phaseIdx], + elemCtx, + boundaryFaceIdx, + pressureCallback); + Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]); + } + + const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx); + auto i = scvf.interiorIndex(); + interiorDofIdx_ = static_cast(i); + exteriorDofIdx_ = -1; + int focusDofIdx = elemCtx.focusDofIndex(); + + // calculate the intrinsic permeability + const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx); + K_ = intQuantsIn.intrinsicPermeability(); + + // correct the pressure gradients by the gravitational acceleration + if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) { + // estimate the gravitational acceleration at a given SCV face + // using the arithmetic mean + const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx); + const auto& posIn = elemCtx.pos(i, timeIdx); + const auto& posFace = scvf.integrationPos(); + + // the distance between the face center and the center of the control volume + DimVector distVecIn(posIn); + distVecIn -= posFace; + Scalar absDist = distVecIn.two_norm(); + Scalar gTimesDist = gIn*distVecIn; + + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + // calculate the hydrostatic pressure at the integration point of the face + Evaluation rhoIn = intQuantsIn.fluidState().density(phaseIdx); + Evaluation pStatIn = - gTimesDist*rhoIn; + + Opm::Valgrind::CheckDefined(pStatIn); + + // compute the hydrostatic gradient between the two control volumes (this + // gradient exhibitis the same direction as the vector between the two + // control volume centers and the length (pStaticExterior - + // pStaticInterior)/distanceInteriorToExterior. Note that for the + // boundary, 'pStaticExterior' is zero as the boundary pressure is + // defined on boundary face's integration point... + EvalDimVector f(distVecIn); + f *= - pStatIn/absDist; + + // calculate the final potential gradient + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + potentialGrad_[phaseIdx][dimIdx] += f[dimIdx]; + + Opm::Valgrind::CheckDefined(potentialGrad_[phaseIdx]); + for (unsigned dimIdx = 0; dimIdx < potentialGrad_[phaseIdx].size(); ++dimIdx) { + if (!Opm::isfinite(potentialGrad_[phaseIdx][dimIdx])) { + throw Opm::NumericalIssue("Non finite potential gradient for phase '" + +std::string(FluidSystem::phaseName(phaseIdx))+"'"); + } + } + } + } + + // determine the upstream and downstream DOFs + const auto& faceNormal = scvf.normal(); + + const auto& matParams = elemCtx.problem().materialLawParams(elemCtx, i, timeIdx); + + Scalar kr[numPhases]; + MaterialLaw::relativePermeabilities(kr, matParams, fluidState); + Opm::Valgrind::CheckDefined(kr); + + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + Evaluation tmp = 0.0; + for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx) + tmp += potentialGrad_[phaseIdx][dimIdx]*faceNormal[dimIdx]; + + if (tmp > 0) { + upstreamDofIdx_[phaseIdx] = exteriorDofIdx_; + downstreamDofIdx_[phaseIdx] = interiorDofIdx_; + } + else { + upstreamDofIdx_[phaseIdx] = interiorDofIdx_; + downstreamDofIdx_[phaseIdx] = exteriorDofIdx_; + } + + // take the phase mobility from the DOF in upstream direction + if (upstreamDofIdx_[phaseIdx] < 0) { + if (interiorDofIdx_ == focusDofIdx) + mobility_[phaseIdx] = + kr[phaseIdx] / fluidState.viscosity(phaseIdx); + else + mobility_[phaseIdx] = + Toolbox::value(kr[phaseIdx]) + / Toolbox::value(fluidState.viscosity(phaseIdx)); + } + else if (upstreamDofIdx_[phaseIdx] != focusDofIdx) + mobility_[phaseIdx] = Toolbox::value(intQuantsIn.mobility(phaseIdx)); + else + mobility_[phaseIdx] = intQuantsIn.mobility(phaseIdx); + Opm::Valgrind::CheckDefined(mobility_[phaseIdx]); + } + } + + /*! + * \brief Calculate the volumetric fluxes of all phases + * + * The pressure potentials and upwind directions must already be + * determined before calling this method! + */ + void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) + { + const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx); + const DimVector& normal = scvf.normal(); + Opm::Valgrind::CheckDefined(normal); + + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + filterVelocity_[phaseIdx] = 0.0; + volumeFlux_[phaseIdx] = 0.0; + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + asImp_().calculateFilterVelocity_(phaseIdx); + Opm::Valgrind::CheckDefined(filterVelocity_[phaseIdx]); + + volumeFlux_[phaseIdx] = 0.0; + for (unsigned i = 0; i < normal.size(); ++i) + volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i]; + } + } + + /*! + * \brief Calculate the volumetric fluxes at a boundary face of all fluid phases + * + * The pressure potentials and upwind directions must already be determined before + * calling this method! + */ + void calculateBoundaryFluxes_(const ElementContext& elemCtx, + unsigned boundaryFaceIdx, + unsigned timeIdx) + { + const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(boundaryFaceIdx); + const DimVector& normal = scvf.normal(); + Opm::Valgrind::CheckDefined(normal); + + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) { + filterVelocity_[phaseIdx] = 0.0; + volumeFlux_[phaseIdx] = 0.0; + continue; + } + + asImp_().calculateFilterVelocity_(phaseIdx); + Opm::Valgrind::CheckDefined(filterVelocity_[phaseIdx]); + volumeFlux_[phaseIdx] = 0.0; + for (unsigned i = 0; i < normal.size(); ++i) + volumeFlux_[phaseIdx] += filterVelocity_[phaseIdx][i] * normal[i]; + } + } + + void calculateFilterVelocity_(unsigned phaseIdx) + { +#ifndef NDEBUG + assert(Opm::isfinite(mobility_[phaseIdx])); + for (unsigned i = 0; i < K_.M(); ++ i) + for (unsigned j = 0; j < K_.N(); ++ j) + assert(std::isfinite(K_[i][j])); +#endif + + K_.mv(potentialGrad_[phaseIdx], filterVelocity_[phaseIdx]); + filterVelocity_[phaseIdx] *= - mobility_[phaseIdx]; + +#ifndef NDEBUG + for (unsigned i = 0; i < filterVelocity_[phaseIdx].size(); ++ i) + assert(Opm::isfinite(filterVelocity_[phaseIdx][i])); +#endif + } + +private: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + +protected: + // intrinsic permeability tensor and its square root + DimMatrix K_; + + // mobilities of all fluid phases [1 / (Pa s)] + Evaluation mobility_[numPhases]; + + // filter velocities of all phases [m/s] + EvalDimVector filterVelocity_[numPhases]; + + // the volumetric flux of all fluid phases over the control + // volume's face [m^3/s / m^2] + Evaluation volumeFlux_[numPhases]; + + // pressure potential gradients of all phases [Pa / m] + EvalDimVector potentialGrad_[numPhases]; + + // upstream, downstream, interior and exterior DOFs + short upstreamDofIdx_[numPhases]; + short downstreamDofIdx_[numPhases]; + short interiorDofIdx_; + short exteriorDofIdx_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/common/diffusionmodule.hh b/opm/models/common/diffusionmodule.hh new file mode 100644 index 000000000..502267156 --- /dev/null +++ b/opm/models/common/diffusionmodule.hh @@ -0,0 +1,496 @@ +// -*- 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 + * + * \brief Classes required for molecular diffusion. + */ +#ifndef EWOMS_DIFFUSION_MODULE_HH +#define EWOMS_DIFFUSION_MODULE_HH + +#include +#include + +#include +#include +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Indices); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup Diffusion + * \class Opm::DiffusionModule + * \brief Provides the auxiliary methods required for consideration of the + * diffusion equation. + */ +template +class DiffusionModule; + +/*! + * \copydoc Opm::DiffusionModule + */ +template +class DiffusionModule +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; + +public: + /*! + * \brief Register all run-time parameters for the diffusion module. + */ + static void registerParameters() + {} + + /*! + * \brief Adds the diffusive mass flux flux to the flux vector over a flux + * integration point. + */ + template + static void addDiffusiveFlux(RateVector& flux OPM_UNUSED, + const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + {} +}; + +/*! + * \copydoc Opm::DiffusionModule + */ +template +class DiffusionModule +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + enum { conti0EqIdx = Indices::conti0EqIdx }; + + typedef Opm::MathToolbox Toolbox; + +public: + /*! + * \brief Register all run-time parameters for the diffusion module. + */ + static void registerParameters() + {} + + /*! + * \brief Adds the mass flux due to molecular diffusion to the flux vector over the + * flux integration point. + */ + template + static void addDiffusiveFlux(RateVector& flux, const Context& context, + unsigned spaceIdx, unsigned timeIdx) + { + const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx); + + const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState(); + const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState(); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + // arithmetic mean of the phase's molar density + Evaluation rhoMolar = fluidStateI.molarDensity(phaseIdx); + rhoMolar += Toolbox::value(fluidStateJ.molarDensity(phaseIdx)); + rhoMolar /= 2; + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) + // mass flux due to molecular diffusion + flux[conti0EqIdx + compIdx] += + -rhoMolar + * extQuants.moleFractionGradientNormal(phaseIdx, compIdx) + * extQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx); + } + } +}; + +/*! + * \ingroup Diffusion + * \class Opm::DiffusionIntensiveQuantities + * + * \brief Provides the volumetric quantities required for the + * calculation of molecular diffusive fluxes. + */ +template +class DiffusionIntensiveQuantities; + +/*! + * \copydoc Opm::DiffusionIntensiveQuantities + */ +template +class DiffusionIntensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + +public: + /*! + * \brief Returns the tortuousity of the sub-domain of a fluid + * phase in the porous medium. + */ + Scalar tortuosity(unsigned phaseIdx OPM_UNUSED) const + { + throw std::logic_error("Method tortuosity() does not make sense " + "if diffusion is disabled"); + } + + /*! + * \brief Returns the molecular diffusion coefficient for a + * component in a phase. + */ + Scalar diffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const + { + throw std::logic_error("Method diffusionCoefficient() does not " + "make sense if diffusion is disabled"); + } + + /*! + * \brief Returns the effective molecular diffusion coefficient of + * the porous medium for a component in a phase. + */ + Scalar effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const + { + throw std::logic_error("Method effectiveDiffusionCoefficient() " + "does not make sense if diffusion is disabled"); + } + +protected: + /*! + * \brief Update the quantities required to calculate diffusive + * mass fluxes. + */ + template + void update_(FluidState& fs OPM_UNUSED, + typename FluidSystem::template ParameterCache& paramCache OPM_UNUSED, + const ElementContext& elemCtx OPM_UNUSED, + unsigned dofIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + { } +}; + +/*! + * \copydoc Opm::DiffusionIntensiveQuantities + */ +template +class DiffusionIntensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + enum { numPhases = FluidSystem::numPhases }; + enum { numComponents = FluidSystem::numComponents }; + +public: + /*! + * \brief Returns the molecular diffusion coefficient for a + * component in a phase. + */ + Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const + { return diffusionCoefficient_[phaseIdx][compIdx]; } + + /*! + * \brief Returns the tortuousity of the sub-domain of a fluid + * phase in the porous medium. + */ + Evaluation tortuosity(unsigned phaseIdx) const + { return tortuosity_[phaseIdx]; } + + /*! + * \brief Returns the effective molecular diffusion coefficient of + * the porous medium for a component in a phase. + */ + Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const + { return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; } + +protected: + /*! + * \brief Update the quantities required to calculate diffusive + * mass fluxes. + */ + template + void update_(FluidState& fluidState, + typename FluidSystem::template ParameterCache& paramCache, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) + { + typedef Opm::MathToolbox Toolbox; + + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + // TODO: let the problem do this (this is a constitutive + // relation of which the model should be free of from the + // abstraction POV!) + const Evaluation& base = + Toolbox::max(0.0001, + intQuants.porosity() + * intQuants.fluidState().saturation(phaseIdx)); + tortuosity_[phaseIdx] = + 1.0 / (intQuants.porosity() * intQuants.porosity()) + * Toolbox::pow(base, 7.0/3.0); + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + diffusionCoefficient_[phaseIdx][compIdx] = + FluidSystem::diffusionCoefficient(fluidState, + paramCache, + phaseIdx, + compIdx); + } + } + } + +private: + Evaluation tortuosity_[numPhases]; + Evaluation diffusionCoefficient_[numPhases][numComponents]; +}; + +/*! + * \ingroup Diffusion + * \class Opm::DiffusionExtensiveQuantities + * + * \brief Provides the quantities required to calculate diffusive mass fluxes. + */ +template +class DiffusionExtensiveQuantities; + +/*! + * \copydoc Opm::DiffusionExtensiveQuantities + */ +template +class DiffusionExtensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + +protected: + /*! + * \brief Update the quantities required to calculate + * the diffusive mass fluxes. + */ + void update_(const ElementContext& elemCtx OPM_UNUSED, + unsigned faceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + {} + + template + void updateBoundary_(const Context& context OPM_UNUSED, + unsigned bfIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED, + const FluidState& fluidState OPM_UNUSED) + {} + +public: + /*! + * \brief The the gradient of the mole fraction times the face normal. + * + * \copydoc Doxygen::phaseIdxParam + * \copydoc Doxygen::compIdxParam + */ + const Evaluation& moleFractionGradientNormal(unsigned phaseIdx OPM_UNUSED, + unsigned compIdx OPM_UNUSED) const + { + throw std::logic_error("The method moleFractionGradient() does not " + "make sense if diffusion is disabled."); + } + + /*! + * \brief The effective diffusion coeffcient of a component in a + * fluid phase at the face's integration point + * + * \copydoc Doxygen::phaseIdxParam + * \copydoc Doxygen::compIdxParam + */ + const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, + unsigned compIdx OPM_UNUSED) const + { + throw std::logic_error("The method effectiveDiffusionCoefficient() " + "does not make sense if diffusion is disabled."); + } +}; + +/*! + * \copydoc Opm::DiffusionExtensiveQuantities + */ +template +class DiffusionExtensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + enum { dimWorld = GridView::dimensionworld }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) }; + + typedef Dune::FieldVector DimVector; + typedef Dune::FieldVector DimEvalVector; + +protected: + /*! + * \brief Update the quantities required to calculate + * the diffusive mass fluxes. + */ + void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx) + { + const auto& gradCalc = elemCtx.gradientCalculator(); + Opm::MoleFractionCallback moleFractionCallback(elemCtx); + + const auto& face = elemCtx.stencil(timeIdx).interiorFace(faceIdx); + const auto& normal = face.normal(); + const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx); + + const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx); + const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + moleFractionCallback.setPhaseIndex(phaseIdx); + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + moleFractionCallback.setComponentIndex(compIdx); + + DimEvalVector moleFractionGradient(0.0); + gradCalc.calculateGradient(moleFractionGradient, + elemCtx, + faceIdx, + moleFractionCallback); + + moleFractionGradientNormal_[phaseIdx][compIdx] = 0.0; + for (unsigned i = 0; i < normal.size(); ++i) + moleFractionGradientNormal_[phaseIdx][compIdx] += + normal[i]*moleFractionGradient[i]; + Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]); + + // use the arithmetic average for the effective + // diffusion coefficients. + effectiveDiffusionCoefficient_[phaseIdx][compIdx] = + (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) + + + intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx)) + / 2; + + Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]); + } + } + } + + template + void updateBoundary_(const Context& context, + unsigned bfIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + const auto& stencil = context.stencil(timeIdx); + const auto& face = stencil.boundaryFace(bfIdx); + + const auto& elemCtx = context.elementContext(); + unsigned insideScvIdx = face.interiorIndex(); + const auto& insideScv = stencil.subControlVolume(insideScvIdx); + + const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx); + const auto& fluidStateInside = intQuantsInside.fluidState(); + + // distance between the center of the SCV and center of the boundary face + DimVector distVec = face.integrationPos(); + distVec -= context.element().geometry().global(insideScv.localGeometry().center()); + + Scalar dist = distVec * face.normal(); + + // if the following assertation triggers, the center of the + // center of the interior SCV was not inside the element! + assert(dist > 0); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) { + // calculate mole fraction gradient using two-point + // gradients + moleFractionGradientNormal_[phaseIdx][compIdx] = + (fluidState.moleFraction(phaseIdx, compIdx) + - + fluidStateInside.moleFraction(phaseIdx, compIdx)) + / dist; + Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]); + + // use effective diffusion coefficients of the interior finite + // volume. + effectiveDiffusionCoefficient_[phaseIdx][compIdx] = + intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx); + + Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]); + } + } + } + +public: + /*! + * \brief The the gradient of the mole fraction times the face normal. + * + * \copydoc Doxygen::phaseIdxParam + * \copydoc Doxygen::compIdxParam + */ + const Evaluation& moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const + { return moleFractionGradientNormal_[phaseIdx][compIdx]; } + + /*! + * \brief The effective diffusion coeffcient of a component in a + * fluid phase at the face's integration point + * + * \copydoc Doxygen::phaseIdxParam + * \copydoc Doxygen::compIdxParam + */ + const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const + { return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; } + +private: + Evaluation moleFractionGradientNormal_[numPhases][numComponents]; + Evaluation effectiveDiffusionCoefficient_[numPhases][numComponents]; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/common/energymodule.hh b/opm/models/common/energymodule.hh new file mode 100644 index 000000000..84e3a142e --- /dev/null +++ b/opm/models/common/energymodule.hh @@ -0,0 +1,859 @@ +// -*- 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 + * + * \brief Contains the classes required to consider energy as a + * conservation quantity in a multi-phase module. + */ +#ifndef EWOMS_ENERGY_MODULE_HH +#define EWOMS_ENERGY_MODULE_HH + +#include +#include + +#include +#include +#include + +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(Indices); +NEW_PROP_TAG(EnableEnergy); +NEW_PROP_TAG(ThermalConductionLaw); +NEW_PROP_TAG(ThermalConductionLawParams); +NEW_PROP_TAG(SolidEnergyLaw); +NEW_PROP_TAG(SolidEnergyLawParams); + +END_PROPERTIES + +namespace Opm { +/*! + * \ingroup Energy + * \brief Provides the auxiliary methods required for consideration of + * the energy equation. + */ +template +class EnergyModule; + +/*! + * \copydoc Opm::EnergyModule + */ +template +class EnergyModule +{ + 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, RateVector) RateVector; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, Model) Model; + + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + + typedef Dune::FieldVector EvalEqVector; + +public: + /*! + * \brief Register all run-time parameters for the energy module. + */ + static void registerParameters() + {} + + /*! + * \brief Returns the name of a primary variable or an empty + * string if the specified primary variable index does not belong to + * the energy module. + */ + static std::string primaryVarName(unsigned pvIdx OPM_UNUSED) + { return ""; } + + /*! + * \brief Returns the name of an equation or an empty + * string if the specified equation index does not belong to + * the energy module. + */ + static std::string eqName(unsigned eqIdx OPM_UNUSED) + { return ""; } + + /*! + * \brief Returns the relative weight of a primary variable for + * calculating relative errors. + */ + static Scalar primaryVarWeight(const Model& model OPM_UNUSED, + unsigned globalDofIdx OPM_UNUSED, + unsigned pvIdx OPM_UNUSED) + { return -1; } + + /*! + * \brief Returns the relative weight of a equation of the residual. + */ + static Scalar eqWeight(const Model& model OPM_UNUSED, + unsigned globalDofIdx OPM_UNUSED, + unsigned eqIdx OPM_UNUSED) + { return -1; } + + /*! + * \brief Given a fluid state, set the temperature in the primary variables + */ + template + static void setPriVarTemperatures(PrimaryVariables& priVars OPM_UNUSED, + const FluidState& fs OPM_UNUSED) + {} + + /*! + * \brief Given a fluid state, set the enthalpy rate which emerges + * from a volumetric rate. + */ + template + static void setEnthalpyRate(RateVector& rateVec OPM_UNUSED, + const FluidState& fluidState OPM_UNUSED, + unsigned phaseIdx OPM_UNUSED, + const Evaluation& volume OPM_UNUSED) + {} + + /*! + * \brief Add the rate of the enthalpy flux to a rate vector. + */ + static void setEnthalpyRate(RateVector& rateVec OPM_UNUSED, + const Evaluation& rate OPM_UNUSED) + {} + + /*! + * \brief Add the rate of the enthalpy flux to a rate vector. + */ + static void addToEnthalpyRate(RateVector& rateVec OPM_UNUSED, + const Evaluation& rate OPM_UNUSED) + {} + + /*! + * \brief Add the rate of the conductive energy flux to a rate vector. + */ + static Scalar thermalConductionRate(const ExtensiveQuantities& extQuants OPM_UNUSED) + { return 0.0; } + + /*! + * \brief Add the energy storage term for a fluid phase to an equation + * vector + */ + template + static void addPhaseStorage(Dune::FieldVector& storage OPM_UNUSED, + const IntensiveQuantities& intQuants OPM_UNUSED, + unsigned phaseIdx OPM_UNUSED) + {} + + /*! + * \brief Add the energy storage term for a fluid phase to an equation + * vector + */ + template + static void addFracturePhaseStorage(Dune::FieldVector& storage OPM_UNUSED, + const IntensiveQuantities& intQuants OPM_UNUSED, + const Scv& scv OPM_UNUSED, + unsigned phaseIdx OPM_UNUSED) + {} + + /*! + * \brief Add the energy storage term for the fracture part a fluid phase to an + * equation vector + */ + template + static void addSolidEnergyStorage(Dune::FieldVector& storage OPM_UNUSED, + const IntensiveQuantities& intQuants OPM_UNUSED) + {} + + /*! + * \brief Evaluates the advective energy fluxes over a face of a + * subcontrol volume and adds the result in the flux vector. + * + * This method is called by compute flux (base class) + */ + template + static void addAdvectiveFlux(RateVector& flux OPM_UNUSED, + const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + {} + + /*! + * \brief Evaluates the advective energy fluxes over a fracture + * which should be attributed to a face of a subcontrol + * volume and adds the result in the flux vector. + */ + template + static void handleFractureFlux(RateVector& flux OPM_UNUSED, + const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + {} + + /*! + * \brief Adds the diffusive energy flux to the flux vector over the face of a + * sub-control volume. + * + * This method is called by compute flux (base class) + */ + template + static void addDiffusiveFlux(RateVector& flux OPM_UNUSED, + const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + {} +}; + +/*! + * \copydoc Opm::EnergyModule + */ +template +class EnergyModule +{ + 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, EqVector) EqVector; + typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector; + typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, Model) Model; + + enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) }; + enum { numPhases = FluidSystem::numPhases }; + enum { energyEqIdx = Indices::energyEqIdx }; + enum { temperatureIdx = Indices::temperatureIdx }; + + typedef Dune::FieldVector EvalEqVector; + typedef Opm::MathToolbox Toolbox; + +public: + /*! + * \brief Register all run-time parameters for the energy module. + */ + static void registerParameters() + {} + + /*! + * \brief Returns the name of a primary variable or an empty + * string if the specified primary variable index does not belong to + * the energy module. + */ + static std::string primaryVarName(unsigned pvIdx) + { + if (pvIdx == temperatureIdx) + return "temperature"; + return ""; + } + + /*! + * \brief Returns the name of an equation or an empty + * string if the specified equation index does not belong to + * the energy module. + */ + static std::string eqName(unsigned eqIdx) + { + if (eqIdx == energyEqIdx) + return "energy"; + return ""; + } + + /*! + * \brief Returns the relative weight of a primary variable for + * calculating relative errors. + */ + static Scalar primaryVarWeight(const Model& model, unsigned globalDofIdx, unsigned pvIdx) + { + if (pvIdx != temperatureIdx) + return -1; + + // make the weight of the temperature primary variable inversly proportional to its value + return std::max(1.0/1000, 1.0/model.solution(/*timeIdx=*/0)[globalDofIdx][temperatureIdx]); + } + + /*! + * \brief Returns the relative weight of a equation. + */ + static Scalar eqWeight(const Model& model OPM_UNUSED, + unsigned globalDofIdx OPM_UNUSED, + unsigned eqIdx) + { + if (eqIdx != energyEqIdx) + return -1; + + // approximate change of internal energy of 1kg of liquid water for a temperature + // change of 30K + return 1.0 / (4.184e3 * 30.0); + } + + /*! + * \brief Set the rate of energy flux of a rate vector. + */ + static void setEnthalpyRate(RateVector& rateVec, const Evaluation& rate) + { rateVec[energyEqIdx] = rate; } + + /*! + * \brief Add the rate of the enthalpy flux to a rate vector. + */ + static void addToEnthalpyRate(RateVector& rateVec, const Evaluation& rate) + { rateVec[energyEqIdx] += rate; } + + /*! + * \brief Returns the conductive energy flux for a given flux integration point. + */ + static Evaluation thermalConductionRate(const ExtensiveQuantities& extQuants) + { return -extQuants.temperatureGradNormal() * extQuants.thermalConductivity(); } + + /*! + * \brief Given a fluid state, set the enthalpy rate which emerges + * from a volumetric rate. + */ + template + static void setEnthalpyRate(RateVector& rateVec, + const FluidState& fluidState, + unsigned phaseIdx, + const Evaluation& volume) + { + rateVec[energyEqIdx] = + volume + * fluidState.density(phaseIdx) + * fluidState.enthalpy(phaseIdx); + } + + /*! + * \brief Given a fluid state, set the temperature in the primary variables + */ + template + static void setPriVarTemperatures(PrimaryVariables& priVars, + const FluidState& fs) + { + priVars[temperatureIdx] = Toolbox::value(fs.temperature(/*phaseIdx=*/0)); +#ifndef NDEBUG + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + assert(std::abs(Toolbox::value(fs.temperature(/*phaseIdx=*/0)) + - Toolbox::value(fs.temperature(phaseIdx))) < 1e-30); + } +#endif + } + + /*! + * \brief Add the energy storage term for a fluid phase to an equation + * vector + */ + template + static void addPhaseStorage(Dune::FieldVector& storage, + const IntensiveQuantities& intQuants, + unsigned phaseIdx) + { + const auto& fs = intQuants.fluidState(); + storage[energyEqIdx] += + Toolbox::template decay(fs.density(phaseIdx)) + * Toolbox::template decay(fs.internalEnergy(phaseIdx)) + * Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(intQuants.porosity()); + } + + /*! + * \brief Add the energy storage term for a fluid phase to an equation + * vector + */ + template + static void addFracturePhaseStorage(Dune::FieldVector& storage, + const IntensiveQuantities& intQuants, + const Scv& scv, + unsigned phaseIdx) + { + const auto& fs = intQuants.fractureFluidState(); + storage[energyEqIdx] += + Toolbox::template decay(fs.density(phaseIdx)) + * Toolbox::template decay(fs.internalEnergy(phaseIdx)) + * Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(intQuants.fracturePorosity()) + * Toolbox::template decay(intQuants.fractureVolume())/scv.volume(); + } + + /*! + * \brief Add the energy storage term for a fluid phase to an equation + * vector + */ + template + static void addSolidEnergyStorage(Dune::FieldVector& storage, + const IntensiveQuantities& intQuants) + { storage[energyEqIdx] += Opm::decay(intQuants.solidInternalEnergy()); } + + /*! + * \brief Evaluates the advective energy fluxes for a flux integration point and adds + * the result in the flux vector. + * + * This method is called by compute flux (base class) + */ + template + static void addAdvectiveFlux(RateVector& flux, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) + { + const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx); + + // advective energy flux in all phases + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!context.model().phaseIsConsidered(phaseIdx)) + continue; + + // intensive quantities of the upstream and the downstream DOFs + unsigned upIdx = static_cast(extQuants.upstreamIndex(phaseIdx)); + const IntensiveQuantities& up = context.intensiveQuantities(upIdx, timeIdx); + + flux[energyEqIdx] += + extQuants.volumeFlux(phaseIdx) + * up.fluidState().enthalpy(phaseIdx) + * up.fluidState().density(phaseIdx); + } + } + + /*! + * \brief Evaluates the advective energy fluxes over a fracture which should be + * attributed to a face of a subcontrol volume and adds the result in the flux + * vector. + */ + template + static void handleFractureFlux(RateVector& flux, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) + { + const auto& scvf = context.stencil(timeIdx).interiorFace(spaceIdx); + const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx); + + // reduce the energy flux in the matrix by the half the width occupied by the + // fracture + flux[energyEqIdx] *= + 1 - extQuants.fractureWidth()/(2*scvf.area()); + + // advective energy flux in all phases + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!context.model().phaseIsConsidered(phaseIdx)) + continue; + + // intensive quantities of the upstream and the downstream DOFs + unsigned upIdx = static_cast(extQuants.upstreamIndex(phaseIdx)); + const IntensiveQuantities& up = context.intensiveQuantities(upIdx, timeIdx); + + flux[energyEqIdx] += + extQuants.fractureVolumeFlux(phaseIdx) + * up.fluidState().enthalpy(phaseIdx) + * up.fluidState().density(phaseIdx); + } + } + + /*! + * \brief Adds the diffusive energy flux to the flux vector over the face of a + * sub-control volume. + * + * This method is called by compute flux (base class) + */ + template + static void addDiffusiveFlux(RateVector& flux, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) + { + const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx); + + // diffusive energy flux + flux[energyEqIdx] += + - extQuants.temperatureGradNormal() + * extQuants.thermalConductivity(); + } +}; + +/*! + * \ingroup Energy + * \class Opm::EnergyIndices + * + * \brief Provides the indices required for the energy equation. + */ +template +struct EnergyIndices; + +/*! + * \copydoc Opm::EnergyIndices + */ +template +struct EnergyIndices +{ + //! The index of the primary variable representing temperature + enum { temperatureIdx = -1000 }; + + //! The index of the equation representing the conservation of energy + enum { energyEqIdx = -1000 }; + +protected: + enum { numEq_ = 0 }; +}; + +/*! + * \copydoc Opm::EnergyIndices + */ +template +struct EnergyIndices +{ + //! The index of the primary variable representing temperature + enum { temperatureIdx = PVOffset }; + + //! The index of the equation representing the conservation of energy + enum { energyEqIdx = PVOffset }; + +protected: + enum { numEq_ = 1 }; +}; + +/*! + * \ingroup Energy + * \class Opm::EnergyIntensiveQuantities + * + * \brief Provides the volumetric quantities required for the energy equation. + */ +template +class EnergyIntensiveQuantities; + +/*! + * \copydoc Opm::EnergyIntensiveQuantities + */ +template +class EnergyIntensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + typedef Opm::MathToolbox Toolbox; + +public: + /*! + * \brief Returns the volumetric internal energy \f$\mathrm{[J/(m^3]}\f$ of the + * solid matrix in the sub-control volume. + */ + Evaluation solidInternalEnergy() const + { + throw std::logic_error("solidInternalEnergy() does not make sense for isothermal models"); + } + + /*! + * \brief Returns the total thermal conductivity \f$\mathrm{[W/m^2 / (K/m)]}\f$ of + * the solid matrix in the sub-control volume. + */ + Evaluation thermalConductivity() const + { + throw std::logic_error("thermalConductivity() does not make sense for isothermal models"); + } + +protected: + /*! + * \brief Update the temperatures of the fluids of a fluid state. + */ + template + static void updateTemperatures_(FluidState& fluidState, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) + { + Scalar T = context.problem().temperature(context, spaceIdx, timeIdx); + fluidState.setTemperature(Toolbox::createConstant(T)); + } + + /*! + * \brief Update the quantities required to calculate + * energy fluxes. + */ + template + void update_(FluidState& fs OPM_UNUSED, + typename FluidSystem::template ParameterCache& paramCache OPM_UNUSED, + const ElementContext& elemCtx OPM_UNUSED, + unsigned dofIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + { } +}; + +/*! + * \copydoc Opm::EnergyIntensiveQuantities + */ +template +class EnergyIntensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, ThermalConductionLaw) ThermalConductionLaw; + typedef typename GET_PROP_TYPE(TypeTag, SolidEnergyLaw) SolidEnergyLaw; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + + enum { numPhases = FluidSystem::numPhases }; + enum { energyEqIdx = Indices::energyEqIdx }; + enum { temperatureIdx = Indices::temperatureIdx }; + + typedef Opm::MathToolbox Toolbox; + +protected: + /*! + * \brief Update the temperatures of the fluids of a fluid state. + */ + template + static void updateTemperatures_(FluidState& fluidState, + const Context& context, + unsigned spaceIdx, + unsigned timeIdx) + { + const auto& priVars = context.primaryVars(spaceIdx, timeIdx); + Evaluation val; + if (std::is_same::value) // finite differences + val = Toolbox::createConstant(priVars[temperatureIdx]); + else { + // automatic differentiation + if (timeIdx == 0) + val = Toolbox::createVariable(priVars[temperatureIdx], temperatureIdx); + else + val = Toolbox::createConstant(priVars[temperatureIdx]); + } + fluidState.setTemperature(val); + } + + /*! + * \brief Update the quantities required to calculate + * energy fluxes. + */ + template + void update_(FluidState& fs, + typename FluidSystem::template ParameterCache& paramCache, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) + { + // set the specific enthalpies of the fluids + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + fs.setEnthalpy(phaseIdx, + FluidSystem::enthalpy(fs, paramCache, phaseIdx)); + } + + // compute and set the volumetric internal energy of the solid phase + const auto& problem = elemCtx.problem(); + const auto& solidEnergyParams = problem.solidEnergyLawParams(elemCtx, dofIdx, timeIdx); + const auto& thermalCondParams = problem.thermalConductionLawParams(elemCtx, dofIdx, timeIdx); + + solidInternalEnergy_ = SolidEnergyLaw::solidInternalEnergy(solidEnergyParams, fs); + thermalConductivity_ = ThermalConductionLaw::thermalConductivity(thermalCondParams, fs); + + Opm::Valgrind::CheckDefined(solidInternalEnergy_); + Opm::Valgrind::CheckDefined(thermalConductivity_); + } + +public: + /*! + * \brief Returns the volumetric internal energy \f$\mathrm{[J/m^3]}\f$ of the + * solid matrix in the sub-control volume. + */ + const Evaluation& solidInternalEnergy() const + { return solidInternalEnergy_; } + + /*! + * \brief Returns the total conductivity capacity \f$\mathrm{[W/m^2 / (K/m)]}\f$ of + * the solid matrix in the sub-control volume. + */ + const Evaluation& thermalConductivity() const + { return thermalConductivity_; } + +private: + Evaluation solidInternalEnergy_; + Evaluation thermalConductivity_; +}; + +/*! + * \ingroup Energy + * \class Opm::EnergyExtensiveQuantities + * + * \brief Provides the quantities required to calculate energy fluxes. + */ +template +class EnergyExtensiveQuantities; + +/*! + * \copydoc Opm::EnergyExtensiveQuantities + */ +template +class EnergyExtensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + +protected: + /*! + * \brief Update the quantities required to calculate + * energy fluxes. + */ + void update_(const ElementContext& elemCtx OPM_UNUSED, + unsigned faceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) + {} + + template + void updateBoundary_(const Context& context OPM_UNUSED, + unsigned bfIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED, + const FluidState& fs OPM_UNUSED) + {} + +public: + /*! + * \brief The temperature gradient times the face normal [K m^2 / m] + */ + Scalar temperatureGradNormal() const + { + throw std::logic_error("Calling temperatureGradNormal() does not make sense " + "for isothermal models"); + } + + /*! + * \brief The total thermal conductivity at the face \f$\mathrm{[W/m^2 / (K/m)]}\f$ + */ + Scalar thermalConductivity() const + { + throw std::logic_error("Calling thermalConductivity() does not make sense for " + "isothermal models"); + } +}; + +/*! + * \copydoc Opm::EnergyExtensiveQuantities + */ +template +class EnergyExtensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + enum { dimWorld = GridView::dimensionworld }; + typedef Dune::FieldVector EvalDimVector; + typedef Dune::FieldVector DimVector; + +protected: + /*! + * \brief Update the quantities required to calculate + * energy fluxes. + */ + void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx) + { + const auto& gradCalc = elemCtx.gradientCalculator(); + Opm::TemperatureCallback temperatureCallback(elemCtx); + + EvalDimVector temperatureGrad; + gradCalc.calculateGradient(temperatureGrad, + elemCtx, + faceIdx, + temperatureCallback); + + // scalar product of temperature gradient and scvf normal + const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(faceIdx); + + temperatureGradNormal_ = 0; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + temperatureGradNormal_ += (face.normal()[dimIdx]*temperatureGrad[dimIdx]); + + const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx); + const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx); + const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx); + + // arithmetic mean + thermalConductivity_ = + 0.5 * (intQuantsInside.thermalConductivity() + intQuantsOutside.thermalConductivity()); + Opm::Valgrind::CheckDefined(thermalConductivity_); + } + + template + void updateBoundary_(const Context& context, unsigned bfIdx, unsigned timeIdx, const FluidState& fs) + { + const auto& stencil = context.stencil(timeIdx); + const auto& face = stencil.boundaryFace(bfIdx); + + const auto& elemCtx = context.elementContext(); + unsigned insideScvIdx = face.interiorIndex(); + const auto& insideScv = elemCtx.stencil(timeIdx).subControlVolume(insideScvIdx); + + const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx); + const auto& fsInside = intQuantsInside.fluidState(); + + // distance between the center of the SCV and center of the boundary face + DimVector distVec = face.integrationPos(); + distVec -= insideScv.geometry().center(); + + Scalar tmp = 0; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + tmp += distVec[dimIdx] * face.normal()[dimIdx]; + Scalar dist = tmp; + + // if the following assertation triggers, the center of the + // center of the interior SCV was not inside the element! + assert(dist > 0); + + // calculate the temperature gradient using two-point gradient + // appoximation + temperatureGradNormal_ = + (fs.temperature(/*phaseIdx=*/0) - fsInside.temperature(/*phaseIdx=*/0)) / dist; + + // take the value for thermal conductivity from the interior finite volume + thermalConductivity_ = intQuantsInside.thermalConductivity(); + } + +public: + /*! + * \brief The temperature gradient times the face normal [K m^2 / m] + */ + const Evaluation& temperatureGradNormal() const + { return temperatureGradNormal_; } + + /*! + * \brief The total thermal conductivity at the face \f$\mathrm{[W/m^2 / + * (K/m)]}\f$ + */ + const Evaluation& thermalConductivity() const + { return thermalConductivity_; } + +private: + Evaluation temperatureGradNormal_; + Evaluation thermalConductivity_; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/common/flux.hh b/opm/models/common/flux.hh new file mode 100644 index 000000000..0310d6bc2 --- /dev/null +++ b/opm/models/common/flux.hh @@ -0,0 +1,35 @@ +// -*- 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 + * + * \brief This file contains the necessary classes to calculate the + * velocity out of a pressure potential gradient. + */ +#ifndef EWOMS_VELOCITY_MODULES_HH +#define EWOMS_VELOCITY_MODULES_HH + +#include +#include + +#endif diff --git a/opm/models/common/forchheimerfluxmodule.hh b/opm/models/common/forchheimerfluxmodule.hh new file mode 100644 index 000000000..683f6550c --- /dev/null +++ b/opm/models/common/forchheimerfluxmodule.hh @@ -0,0 +1,583 @@ +// -*- 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 + * + * \brief This file contains the necessary classes to calculate the + * volumetric fluxes out of a pressure potential gradient using the + * Forchhheimer approach. + */ +#ifndef EWOMS_FORCHHEIMER_FLUX_MODULE_HH +#define EWOMS_FORCHHEIMER_FLUX_MODULE_HH + +#include "darcyfluxmodule.hh" + +#include + +#include +#include +#include + +#include +#include + +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(MaterialLaw); + +END_PROPERTIES + +namespace Opm { +template +class ForchheimerIntensiveQuantities; + +template +class ForchheimerExtensiveQuantities; + +template +class ForchheimerBaseProblem; + +/*! + * \ingroup FluxModules + * \brief Specifies a flux module which uses the Forchheimer relation. + */ +template +struct ForchheimerFluxModule +{ + typedef ForchheimerIntensiveQuantities FluxIntensiveQuantities; + typedef ForchheimerExtensiveQuantities FluxExtensiveQuantities; + typedef ForchheimerBaseProblem FluxBaseProblem; + + /*! + * \brief Register all run-time parameters for the flux module. + */ + static void registerParameters() + {} +}; + +/*! + * \ingroup FluxModules + * \brief Provides the defaults for the parameters required by the + * Forchheimer velocity approach. + */ +template +class ForchheimerBaseProblem +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + +public: + /*! + * \brief Returns the Ergun coefficient. + * + * The Ergun coefficient is a measure how much the velocity is + * reduced by turbolence. It is a quantity that does not depend on + * the fluid phase but only on the porous medium in question. A + * value of 0 means that the velocity is not influenced by + * turbolence. + */ + template + Scalar ergunCoefficient(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + throw std::logic_error("Not implemented: Problem::ergunCoefficient()"); + } + + /*! + * \brief Returns the ratio between the phase mobility + * \f$k_{r,\alpha}\f$ and its passability + * \f$\eta_{r,\alpha}\f$ for a given fluid phase + * \f$\alpha\f$. + * + * The passability coefficient specifies the influence of the + * other fluid phases on the turbolent behaviour of a given fluid + * phase. By default it is equal to the relative permeability. The + * mobility to passability ratio is the inverse of phase' the viscosity. + */ + template + Evaluation mobilityPassabilityRatio(Context& context, + unsigned spaceIdx, + unsigned timeIdx, + unsigned phaseIdx) const + { + return 1.0 / context.intensiveQuantities(spaceIdx, timeIdx).fluidState().viscosity(phaseIdx); + } +}; + +/*! + * \ingroup FluxModules + * \brief Provides the intensive quantities for the Forchheimer module + */ +template +class ForchheimerIntensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + +public: + /*! + * \brief Returns the Ergun coefficient. + * + * The Ergun coefficient is a measure how much the velocity is + * reduced by turbolence. A value of 0 means that it is not + * influenced. + */ + const Evaluation& ergunCoefficient() const + { return ergunCoefficient_; } + + /*! + * \brief Returns the passability of a phase. + */ + const Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const + { return mobilityPassabilityRatio_[phaseIdx]; } + +protected: + void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx) + { + const auto& problem = elemCtx.problem(); + ergunCoefficient_ = problem.ergunCoefficient(elemCtx, dofIdx, timeIdx); + + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) + mobilityPassabilityRatio_[phaseIdx] = + problem.mobilityPassabilityRatio(elemCtx, + dofIdx, + timeIdx, + phaseIdx); + } + +private: + Evaluation ergunCoefficient_; + Evaluation mobilityPassabilityRatio_[numPhases]; +}; + +/*! + * \ingroup FluxModules + * \brief Provides the Forchheimer flux module + * + * The commonly used Darcy relation looses its validity for Reynolds numbers \f$ Re < + * 1\f$. If one encounters flow velocities in porous media above this threshold, the + * Forchheimer approach can be used. Like the Darcy approach, it is a relation of with + * the fluid velocity in terms of the gradient of pressure potential. However, this + * relation is not linear (as in the Darcy case) any more. + * + * Therefore, the Newton scheme is used to solve the Forchheimer equation. This velocity + * is then used like the Darcy velocity e.g. by the local residual. + * + * Note that for Reynolds numbers above \f$\approx 500\f$ the standard Forchheimer + * relation also looses it's validity. + * + * The Forchheimer equation is given by the following relation: + * + * \f[ + \nabla p_\alpha - \rho_\alpha \vec{g} = + - \frac{\mu_\alpha}{k_{r,\alpha}} K^{-1}\vec{v}_\alpha + - \frac{\rho_\alpha C_E}{\eta_{r,\alpha}} \sqrt{K}^{-1} + \left| \vec{v}_\alpha \right| \vec{v}_\alpha + \f] + * + * Where \f$C_E\f$ is the modified Ergun parameter and \f$\eta_{r,\alpha}\f$ is the + * passability which is given by a closure relation (usually it is assumed to be + * identical to the relative permeability). To avoid numerical problems, the relation + * implemented by this class multiplies both sides with \f$\frac{k_{r_alpha}}{mu} K\f$, + * so we get + * + * \f[ + \frac{k_{r_alpha}}{mu} K \left( \nabla p_\alpha - \rho_\alpha \vec{g}\right) = + - \vec{v}_\alpha + - \frac{\rho_\alpha C_E}{\eta_{r,\alpha}} \frac{k_{r_alpha}}{mu} \sqrt{K} + \left| \vec{v}_\alpha \right| \vec{v}_\alpha + \f] + + */ +template +class ForchheimerExtensiveQuantities + : public DarcyExtensiveQuantities +{ + typedef DarcyExtensiveQuantities DarcyExtQuants; + 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, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation; + + enum { dimWorld = GridView::dimensionworld }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + + typedef Opm::MathToolbox Toolbox; + + typedef Dune::FieldVector DimVector; + typedef Dune::FieldVector DimEvalVector; + typedef Dune::FieldMatrix DimMatrix; + typedef Dune::FieldMatrix DimEvalMatrix; + +public: + /*! + * \brief Return the Ergun coefficent at the face's integration point. + */ + const Evaluation& ergunCoefficient() const + { return ergunCoefficient_; } + + /*! + * \brief Return the ratio of the mobility divided by the passability at the face's + * integration point for a given fluid phase. + * + * Usually, that's the inverse of the viscosity. + */ + Evaluation& mobilityPassabilityRatio(unsigned phaseIdx) const + { return mobilityPassabilityRatio_[phaseIdx]; } + +protected: + void calculateGradients_(const ElementContext& elemCtx, + unsigned faceIdx, + unsigned timeIdx) + { + DarcyExtQuants::calculateGradients_(elemCtx, faceIdx, timeIdx); + + auto focusDofIdx = elemCtx.focusDofIndex(); + unsigned i = static_cast(this->interiorDofIdx_); + unsigned j = static_cast(this->exteriorDofIdx_); + const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx); + const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx); + + // calculate the square root of the intrinsic permeability + assert(isDiagonal_(this->K_)); + sqrtK_ = 0.0; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]); + + // obtain the Ergun coefficient. Lacking better ideas, we use its the arithmetic mean. + if (focusDofIdx == i) { + ergunCoefficient_ = + (intQuantsIn.ergunCoefficient() + + Opm::getValue(intQuantsEx.ergunCoefficient()))/2; + } + else if (focusDofIdx == j) + ergunCoefficient_ = + (Opm::getValue(intQuantsIn.ergunCoefficient()) + + intQuantsEx.ergunCoefficient())/2; + else + ergunCoefficient_ = + (Opm::getValue(intQuantsIn.ergunCoefficient()) + + Opm::getValue(intQuantsEx.ergunCoefficient()))/2; + + // obtain the mobility to passability ratio for each phase. + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + unsigned upIdx = static_cast(this->upstreamIndex_(phaseIdx)); + const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx); + + if (focusDofIdx == upIdx) { + density_[phaseIdx] = + up.fluidState().density(phaseIdx); + mobilityPassabilityRatio_[phaseIdx] = + up.mobilityPassabilityRatio(phaseIdx); + } + else { + density_[phaseIdx] = + Opm::getValue(up.fluidState().density(phaseIdx)); + mobilityPassabilityRatio_[phaseIdx] = + Opm::getValue(up.mobilityPassabilityRatio(phaseIdx)); + } + } + } + + template + void calculateBoundaryGradients_(const ElementContext& elemCtx, + unsigned boundaryFaceIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + DarcyExtQuants::calculateBoundaryGradients_(elemCtx, + boundaryFaceIdx, + timeIdx, + fluidState); + + auto focusDofIdx = elemCtx.focusDofIndex(); + unsigned i = static_cast(this->interiorDofIdx_); + const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx); + + // obtain the Ergun coefficient. Because we are on the boundary here, we will + // take the Ergun coefficient of the interior + if (focusDofIdx == i) + ergunCoefficient_ = intQuantsIn.ergunCoefficient(); + else + ergunCoefficient_ = Opm::getValue(intQuantsIn.ergunCoefficient()); + + // calculate the square root of the intrinsic permeability + assert(isDiagonal_(this->K_)); + sqrtK_ = 0.0; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + sqrtK_[dimIdx] = std::sqrt(this->K_[dimIdx][dimIdx]); + + for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) + continue; + + if (focusDofIdx == i) { + density_[phaseIdx] = intQuantsIn.fluidState().density(phaseIdx); + mobilityPassabilityRatio_[phaseIdx] = intQuantsIn.mobilityPassabilityRatio(phaseIdx); + } + else { + density_[phaseIdx] = + Opm::getValue(intQuantsIn.fluidState().density(phaseIdx)); + mobilityPassabilityRatio_[phaseIdx] = + Opm::getValue(intQuantsIn.mobilityPassabilityRatio(phaseIdx)); + } + } + } + + /*! + * \brief Calculate the volumetric fluxes of all phases + * + * The pressure potentials and upwind directions must already be + * determined before calling this method! + */ + void calculateFluxes_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) + { + auto focusDofIdx = elemCtx.focusDofIndex(); + auto i = asImp_().interiorIndex(); + auto j = asImp_().exteriorIndex(); + const auto& intQuantsI = elemCtx.intensiveQuantities(i, timeIdx); + const auto& intQuantsJ = elemCtx.intensiveQuantities(j, timeIdx); + + const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx); + const auto& normal = scvf.normal(); + Opm::Valgrind::CheckDefined(normal); + + // obtain the Ergun coefficient from the intensive quantity object. Until a + // better method comes along, we use arithmetic averaging. + if (focusDofIdx == i) + ergunCoefficient_ = + (intQuantsI.ergunCoefficient() + + Opm::getValue(intQuantsJ.ergunCoefficient())) / 2; + else if (focusDofIdx == j) + ergunCoefficient_ = + (Opm::getValue(intQuantsI.ergunCoefficient()) + + intQuantsJ.ergunCoefficient()) / 2; + else + ergunCoefficient_ = + (Opm::getValue(intQuantsI.ergunCoefficient()) + + Opm::getValue(intQuantsJ.ergunCoefficient())) / 2; + + /////////////// + // calculate the weights of the upstream and the downstream control volumes + /////////////// + for (unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) { + this->filterVelocity_[phaseIdx] = 0.0; + this->volumeFlux_[phaseIdx] = 0.0; + continue; + } + + calculateForchheimerFlux_(phaseIdx); + + this->volumeFlux_[phaseIdx] = 0.0; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx) + this->volumeFlux_[phaseIdx] += + this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx]; + } + } + + /*! + * \brief Calculate the volumetric flux rates of all phases at the domain boundary + */ + void calculateBoundaryFluxes_(const ElementContext& elemCtx, + unsigned bfIdx, + unsigned timeIdx) + { + const auto& boundaryFace = elemCtx.stencil(timeIdx).boundaryFace(bfIdx); + const auto& normal = boundaryFace.normal(); + + /////////////// + // calculate the weights of the upstream and the downstream degrees of freedom + /////////////// + for (unsigned phaseIdx = 0; phaseIdx < numPhases; phaseIdx++) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) { + this->filterVelocity_[phaseIdx] = 0.0; + this->volumeFlux_[phaseIdx] = 0.0; + continue; + } + + calculateForchheimerFlux_(phaseIdx); + + this->volumeFlux_[phaseIdx] = 0.0; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + this->volumeFlux_[phaseIdx] += + this->filterVelocity_[phaseIdx][dimIdx]*normal[dimIdx]; + } + } + + void calculateForchheimerFlux_(unsigned phaseIdx) + { + // initial guess: filter velocity is zero + DimEvalVector& velocity = this->filterVelocity_[phaseIdx]; + velocity = 0.0; + + // the change of velocity between two consecutive Newton iterations + DimEvalVector deltaV(1e5); + // the function value that is to be minimized of the equation that is to be + // fulfilled + DimEvalVector residual; + // derivative of equation that is to be solved + DimEvalMatrix gradResid; + + // search by means of the Newton method for a root of Forchheimer equation + unsigned newtonIter = 0; + while (deltaV.one_norm() > 1e-11) { + if (newtonIter >= 50) + throw Opm::NumericalIssue("Could not determine Forchheimer velocity within " + +std::to_string(newtonIter)+" iterations"); + ++newtonIter; + + // calculate the residual and its Jacobian matrix + gradForchheimerResid_(residual, gradResid, phaseIdx); + + // newton method + gradResid.solve(deltaV, residual); + velocity -= deltaV; + } + } + + void forchheimerResid_(DimEvalVector& residual, unsigned phaseIdx) const + { + const DimEvalVector& velocity = this->filterVelocity_[phaseIdx]; + + // Obtaining the upstreamed quantities + const auto& mobility = this->mobility_[phaseIdx]; + const auto& density = density_[phaseIdx]; + const auto& mobilityPassabilityRatio = mobilityPassabilityRatio_[phaseIdx]; + + // optain the quantites for the integration point + const auto& pGrad = this->potentialGrad_[phaseIdx]; + + // residual = v_\alpha + residual = velocity; + + // residual += mobility_\alpha K(\grad p_\alpha - \rho_\alpha g) + // -> this->K_.usmv(mobility, pGrad, residual); + assert(isDiagonal_(this->K_)); + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++ dimIdx) + residual[dimIdx] += mobility*pGrad[dimIdx]*this->K_[dimIdx][dimIdx]; + + // Forchheimer turbulence correction: + // + // residual += + // \rho_\alpha + // * mobility_\alpha + // * C_E / \eta_{r,\alpha} + // * abs(v_\alpha) * sqrt(K)*v_\alpha + // + // -> sqrtK_.usmv(density*mobilityPassabilityRatio*ergunCoefficient_*velocity.two_norm(), + // velocity, + // residual); + Evaluation absVel = 0.0; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + absVel += velocity[dimIdx]*velocity[dimIdx]; + // the derivatives of the square root of 0 are undefined, so we must guard + // against this case + if (absVel <= 0.0) + absVel = 0.0; + else + absVel = Toolbox::sqrt(absVel); + const auto& alpha = density*mobilityPassabilityRatio*ergunCoefficient_*absVel; + for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) + residual[dimIdx] += sqrtK_[dimIdx]*alpha*velocity[dimIdx]; + Opm::Valgrind::CheckDefined(residual); + } + + void gradForchheimerResid_(DimEvalVector& residual, + DimEvalMatrix& gradResid, + unsigned phaseIdx) + { + // TODO (?) use AD for this. + DimEvalVector& velocity = this->filterVelocity_[phaseIdx]; + forchheimerResid_(residual, phaseIdx); + + Scalar eps = 1e-11; + DimEvalVector tmp; + for (unsigned i = 0; i < dimWorld; ++i) { + Scalar coordEps = std::max(eps, Toolbox::scalarValue(velocity[i]) * (1 + eps)); + velocity[i] += coordEps; + forchheimerResid_(tmp, phaseIdx); + tmp -= residual; + tmp /= coordEps; + gradResid[i] = tmp; + velocity[i] -= coordEps; + } + } + + /*! + * \brief Check whether all off-diagonal entries of a tensor are zero. + * + * \param K the tensor that is to be checked. + * \return True iff all off-diagonals are zero. + * + */ + bool isDiagonal_(const DimMatrix& K) const + { + for (unsigned i = 0; i < dimWorld; i++) { + for (unsigned j = 0; j < dimWorld; j++) { + if (i == j) + continue; + + if (std::abs(K[i][j]) > 1e-25) + return false; + } + } + return true; + } + +private: + Implementation& asImp_() + { return *static_cast(this); } + + const Implementation& asImp_() const + { return *static_cast(this); } + +protected: + // intrinsic permeability tensor and its square root + DimVector sqrtK_; + + // Ergun coefficient of all phases at the integration point + Evaluation ergunCoefficient_; + + // Passability of all phases at the integration point + Evaluation mobilityPassabilityRatio_[numPhases]; + + // Density of all phases at the integration point + Evaluation density_[numPhases]; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/common/multiphasebaseextensivequantities.hh b/opm/models/common/multiphasebaseextensivequantities.hh new file mode 100644 index 000000000..4caf0c486 --- /dev/null +++ b/opm/models/common/multiphasebaseextensivequantities.hh @@ -0,0 +1,177 @@ +// -*- 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::MultiPhaseBaseExtensiveQuantities + */ +#ifndef EWOMS_MULTI_PHASE_BASE_EXTENSIVE_QUANTITIES_HH +#define EWOMS_MULTI_PHASE_BASE_EXTENSIVE_QUANTITIES_HH + +#include "multiphasebaseproperties.hh" + +#include +#include +#include + +#include +#include + +#include + +namespace Opm { +/*! + * \ingroup Discretization + * + * \brief This class calculates the pressure potential gradients and + * the filter velocities for multi-phase flow in porous media + */ +template +class MultiPhaseBaseExtensiveQuantities + : public GET_PROP_TYPE(TypeTag, DiscExtensiveQuantities) + , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxExtensiveQuantities +{ + typedef typename GET_PROP_TYPE(TypeTag, DiscExtensiveQuantities) ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + + typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule; + typedef typename FluxModule::FluxExtensiveQuantities FluxExtensiveQuantities; + +public: + /*! + * \brief Register all run-time parameters for the extensive quantities. + */ + static void registerParameters() + { + FluxModule::registerParameters(); + } + + /*! + * \brief Update the extensive quantities for a given sub-control-volume-face. + * + * \param elemCtx Reference to the current element context. + * \param scvfIdx The local index of the sub-control-volume face for + * which the extensive quantities should be calculated. + * \param timeIdx The index used by the time discretization. + */ + void update(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx) + { + ParentType::update(elemCtx, scvfIdx, timeIdx); + + // compute the pressure potential gradients + FluxExtensiveQuantities::calculateGradients_(elemCtx, scvfIdx, timeIdx); + + // Check whether the pressure potential is in the same direction as the face + // normal or in the opposite one + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + if (!elemCtx.model().phaseIsConsidered(phaseIdx)) { + Opm::Valgrind::SetUndefined(upstreamScvIdx_[phaseIdx]); + Opm::Valgrind::SetUndefined(downstreamScvIdx_[phaseIdx]); + continue; + } + + upstreamScvIdx_[phaseIdx] = FluxExtensiveQuantities::upstreamIndex_(phaseIdx); + downstreamScvIdx_[phaseIdx] = FluxExtensiveQuantities::downstreamIndex_(phaseIdx); + } + + FluxExtensiveQuantities::calculateFluxes_(elemCtx, scvfIdx, timeIdx); + } + + + /*! + * \brief Update the extensive quantities for a given boundary face. + * + * \param context Reference to the current execution context. + * \param bfIdx The local index of the boundary face for which + * the extensive quantities should be calculated. + * \param timeIdx The index used by the time discretization. + * \param fluidState The FluidState on the domain boundary. + * \param paramCache The FluidSystem's parameter cache. + */ + template + void updateBoundary(const Context& context, + unsigned bfIdx, + unsigned timeIdx, + const FluidState& fluidState) + { + ParentType::updateBoundary(context, bfIdx, timeIdx, fluidState); + + FluxExtensiveQuantities::calculateBoundaryGradients_(context.elementContext(), + bfIdx, + timeIdx, + fluidState); + FluxExtensiveQuantities::calculateBoundaryFluxes_(context.elementContext(), + bfIdx, + timeIdx); + } + + /*! + * \brief Return the local index of the upstream control volume for a given phase as + * a function of the normal flux. + * + * \param phaseIdx The index of the fluid phase for which the upstream + * direction is requested. + */ + short upstreamIndex(unsigned phaseIdx) const + { return upstreamScvIdx_[phaseIdx]; } + + /*! + * \brief Return the local index of the downstream control volume + * for a given phase as a function of the normal flux. + * + * \param phaseIdx The index of the fluid phase for which the downstream + * direction is requested. + */ + short downstreamIndex(unsigned phaseIdx) const + { return downstreamScvIdx_[phaseIdx]; } + + /*! + * \brief Return the weight of the upstream control volume + * for a given phase as a function of the normal flux. + * + * \param phaseIdx The index of the fluid phase + */ + Scalar upstreamWeight(unsigned phaseIdx OPM_UNUSED) const + { return 1.0; } + + /*! + * \brief Return the weight of the downstream control volume + * for a given phase as a function of the normal flux. + * + * \param phaseIdx The index of the fluid phase + */ + Scalar downstreamWeight(unsigned phaseIdx) const + { return 1.0 - upstreamWeight(phaseIdx); } + +private: + short upstreamScvIdx_[numPhases]; + short downstreamScvIdx_[numPhases]; +}; + +} // namespace Opm + +#endif diff --git a/opm/models/common/multiphasebasemodel.hh b/opm/models/common/multiphasebasemodel.hh new file mode 100644 index 000000000..147501d7d --- /dev/null +++ b/opm/models/common/multiphasebasemodel.hh @@ -0,0 +1,253 @@ +// -*- 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::MultiPhaseBaseModel + */ +#ifndef EWOMS_MULTI_PHASE_BASE_MODEL_HH +#define EWOMS_MULTI_PHASE_BASE_MODEL_HH + +#include + +#include "multiphasebaseproperties.hh" +#include "multiphasebaseproblem.hh" +#include "multiphasebaseextensivequantities.hh" + +#include +#include + +#include +#include +#include +#include +#include + +namespace Opm { +template +class MultiPhaseBaseModel; +} + +BEGIN_PROPERTIES + +//! The generic type tag for problems using the immiscible multi-phase model +NEW_TYPE_TAG(MultiPhaseBaseModel, INHERITS_FROM(VtkMultiPhase, VtkTemperature)); + +//! Specify the splices of the MultiPhaseBaseModel type tag +SET_SPLICES(MultiPhaseBaseModel, SpatialDiscretizationSplice); + +//! Set the default spatial discretization +//! +//! We use a vertex centered finite volume method by default +SET_TAG_PROP(MultiPhaseBaseModel, SpatialDiscretizationSplice, VcfvDiscretization); + +//! set the number of equations to the number of phases +SET_INT_PROP(MultiPhaseBaseModel, NumEq, GET_PROP_TYPE(TypeTag, Indices)::numEq); +//! The number of phases is determined by the fluid system +SET_INT_PROP(MultiPhaseBaseModel, NumPhases, GET_PROP_TYPE(TypeTag, FluidSystem)::numPhases); +//! Number of chemical species in the system +SET_INT_PROP(MultiPhaseBaseModel, NumComponents, GET_PROP_TYPE(TypeTag, FluidSystem)::numComponents); + +//! The type of the base base class for actual problems +SET_TYPE_PROP(MultiPhaseBaseModel, BaseProblem, Opm::MultiPhaseBaseProblem); + +//! By default, use the Darcy relation to determine the phase velocity +SET_TYPE_PROP(MultiPhaseBaseModel, FluxModule, Opm::DarcyFluxModule); + +/*! + * \brief Set the material law to the null law by default. + */ +SET_PROP(MultiPhaseBaseModel, MaterialLaw) +{ +private: + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef Opm::NullMaterialTraits Traits; + +public: + typedef Opm::NullMaterial type; +}; + +/*! + * \brief Set the property for the material parameters by extracting + * it from the material law. + */ +SET_TYPE_PROP(MultiPhaseBaseModel, + MaterialLawParams, + typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params); + +//! set the energy storage law for the solid to the one which assumes zero heat capacity +//! by default +SET_TYPE_PROP(MultiPhaseBaseModel, + SolidEnergyLaw, + Opm::NullSolidEnergyLaw); + +//! extract the type of the parameter objects for the solid energy storage law from the +//! law itself +SET_TYPE_PROP(MultiPhaseBaseModel, + SolidEnergyLawParams, + typename GET_PROP_TYPE(TypeTag, SolidEnergyLaw)::Params); + +//! set the thermal conduction law to a dummy one by default +SET_TYPE_PROP(MultiPhaseBaseModel, + ThermalConductionLaw, + Opm::NullThermalConductionLaw); + +//! extract the type of the parameter objects for the thermal conduction law from the law +//! itself +SET_TYPE_PROP(MultiPhaseBaseModel, + ThermalConductionLawParams, + typename GET_PROP_TYPE(TypeTag, ThermalConductionLaw)::Params); + +//! disable gravity by default +SET_BOOL_PROP(MultiPhaseBaseModel, EnableGravity, false); + + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup MultiPhaseBaseModel + * \brief A base class for fully-implicit multi-phase porous-media flow models + * which assume multiple fluid phases. + */ +template +class MultiPhaseBaseModel : public GET_PROP_TYPE(TypeTag, Discretization) +{ + typedef typename GET_PROP_TYPE(TypeTag, Discretization) ParentType; + typedef typename GET_PROP_TYPE(TypeTag, Model) Implementation; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, ThreadManager) ThreadManager; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices; + typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef typename GridView::template Codim<0>::Iterator ElementIterator; + typedef typename GridView::template Codim<0>::Entity Element; + + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + enum { numComponents = FluidSystem::numComponents }; + +public: + MultiPhaseBaseModel(Simulator& simulator) + : ParentType(simulator) + { } + + /*! + * \brief Register all run-time parameters for the immiscible model. + */ + static void registerParameters() + { + ParentType::registerParameters(); + + // register runtime parameters of the VTK output modules + Opm::VtkMultiPhaseModule::registerParameters(); + Opm::VtkTemperatureModule::registerParameters(); + } + + /*! + * \brief Returns true iff a fluid phase is used by the model. + * + * \param phaseIdx The index of the fluid phase in question + */ + bool phaseIsConsidered(unsigned phaseIdx OPM_UNUSED) const + { return true; } + + /*! + * \brief Compute the total storage inside one phase of all + * conservation quantities. + * + * \copydetails Doxygen::storageParam + * \copydetails Doxygen::phaseIdxParam + */ + void globalPhaseStorage(EqVector& storage, unsigned phaseIdx) + { + assert(0 <= phaseIdx && phaseIdx < numPhases); + + storage = 0; + + ThreadedEntityIterator threadedElemIt(this->gridView()); + std::mutex mutex; +#ifdef _OPENMP +#pragma omp parallel +#endif + { + // Attention: the variables below are thread specific and thus cannot be + // moved in front of the #pragma! + unsigned threadId = ThreadManager::threadId(); + ElementContext elemCtx(this->simulator_); + ElementIterator elemIt = threadedElemIt.beginParallel(); + EqVector tmp; + + for (; !threadedElemIt.isFinished(elemIt); elemIt = threadedElemIt.increment()) { + const Element& elem = *elemIt; + if (elem.partitionType() != Dune::InteriorEntity) + continue; // ignore ghost and overlap elements + + elemCtx.updateStencil(elem); + elemCtx.updateIntensiveQuantities(/*timeIdx=*/0); + + const auto& stencil = elemCtx.stencil(/*timeIdx=*/0); + + for (unsigned dofIdx = 0; dofIdx < elemCtx.numDof(/*timeIdx=*/0); ++dofIdx) { + const auto& scv = stencil.subControlVolume(dofIdx); + const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0); + + tmp = 0; + this->localResidual(threadId).addPhaseStorage(tmp, + elemCtx, + dofIdx, + /*timeIdx=*/0, + phaseIdx); + tmp *= scv.volume()*intQuants.extrusionFactor(); + + mutex.lock(); + storage += tmp; + mutex.unlock(); + } + } + } + + storage = this->gridView_.comm().sum(storage); + } + + void registerOutputModules_() + { + ParentType::registerOutputModules_(); + + // add the VTK output modules which make sense for all multi-phase models + this->addOutputModule(new Opm::VtkMultiPhaseModule(this->simulator_)); + this->addOutputModule(new Opm::VtkTemperatureModule(this->simulator_)); + } + +private: + const Implementation& asImp_() const + { return *static_cast(this); } +}; +} // namespace Opm + +#endif diff --git a/opm/models/common/multiphasebaseproblem.hh b/opm/models/common/multiphasebaseproblem.hh new file mode 100644 index 000000000..36564250b --- /dev/null +++ b/opm/models/common/multiphasebaseproblem.hh @@ -0,0 +1,409 @@ +// -*- 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::MultiPhaseBaseProblem + */ +#ifndef EWOMS_MULTI_PHASE_BASE_PROBLEM_HH +#define EWOMS_MULTI_PHASE_BASE_PROBLEM_HH + +#include "multiphasebaseproperties.hh" + +#include +#include + +#include +#include +#include +#include + +#include +#include + +BEGIN_PROPERTIES + +NEW_PROP_TAG(SolidEnergyLawParams); +NEW_PROP_TAG(ThermalConductionLawParams); +NEW_PROP_TAG(EnableGravity); +NEW_PROP_TAG(FluxModule); + +END_PROPERTIES + +namespace Opm { + +/*! + * \ingroup Discretization + * + * \brief The base class for the problems of ECFV discretizations which deal + * with a multi-phase flow through a porous medium. + */ +template +class MultiPhaseBaseProblem + : public FvBaseProblem + , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxBaseProblem +{ +//! \cond SKIP_THIS + typedef Opm::FvBaseProblem ParentType; + + typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation; + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator; + typedef typename GET_PROP_TYPE(TypeTag, SolidEnergyLawParams) SolidEnergyLawParams; + typedef typename GET_PROP_TYPE(TypeTag, ThermalConductionLawParams) ThermalConductionLawParams; + typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params MaterialLawParams; + + enum { dimWorld = GridView::dimensionworld }; + enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) }; + typedef Dune::FieldVector DimVector; + typedef Dune::FieldMatrix DimMatrix; +//! \endcond + +public: + /*! + * \copydoc Problem::FvBaseProblem(Simulator& ) + */ + MultiPhaseBaseProblem(Simulator& simulator) + : ParentType(simulator) + { init_(); } + + /*! + * \brief Register all run-time parameters for the problem and the model. + */ + static void registerParameters() + { + ParentType::registerParameters(); + + EWOMS_REGISTER_PARAM(TypeTag, bool, EnableGravity, + "Use the gravity correction for the pressure gradients."); + } + + /*! + * \brief Returns the intrinsic permeability of an intersection. + * + * This method is specific to the finite volume discretizations. If left unspecified, + * it calls the intrinsicPermeability() method for the intersection's interior and + * exterior finite volumes and averages them harmonically. Note that if this function + * is defined, the intrinsicPermeability() method does not need to be defined by the + * problem (if a finite-volume discretization is used). + */ + template + void intersectionIntrinsicPermeability(DimMatrix& result, + const Context& context, + unsigned intersectionIdx, + unsigned timeIdx) const + { + const auto& scvf = context.stencil(timeIdx).interiorFace(intersectionIdx); + + const DimMatrix& K1 = asImp_().intrinsicPermeability(context, scvf.interiorIndex(), timeIdx); + const DimMatrix& K2 = asImp_().intrinsicPermeability(context, scvf.exteriorIndex(), timeIdx); + + // entry-wise harmonic mean. this is almost certainly wrong if + // you have off-main diagonal entries in your permeabilities! + for (unsigned i = 0; i < dimWorld; ++i) + for (unsigned j = 0; j < dimWorld; ++j) + result[i][j] = Opm::harmonicMean(K1[i][j], K2[i][j]); + } + + /*! + * \name Problem parameters + */ + // \{ + + /*! + * \brief Returns the intrinsic permeability tensor \f$[m^2]\f$ at a given position + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + const DimMatrix& intrinsicPermeability(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + throw std::logic_error("Not implemented: Problem::intrinsicPermeability()"); + } + + /*! + * \brief Returns the porosity [] of the porous medium for a given + * control volume. + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + Scalar porosity(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + throw std::logic_error("Not implemented: Problem::porosity()"); + } + + /*! + * \brief Returns the parameter object for the energy storage law of the solid in a + * sub-control volume. + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + const SolidEnergyLawParams& + solidEnergyParams(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + throw std::logic_error("Not implemented: Problem::solidEnergyParams()"); + } + + /*! + * \brief Returns the parameter object for the thermal conductivity law in a + * sub-control volume. + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + const ThermalConductionLawParams& + thermalConductionParams(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + throw std::logic_error("Not implemented: Problem::thermalConductionParams()"); + } + + /*! + * \brief Define the tortuosity. + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + Scalar tortuosity(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + throw std::logic_error("Not implemented: Problem::tortuosity()"); + } + + /*! + * \brief Define the dispersivity. + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + Scalar dispersivity(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + throw std::logic_error("Not implemented: Problem::dispersivity()"); + } + + /*! + * \brief Returns the material law parameters \f$\mathrm{[K]}\f$ within a control volume. + * + * If you get a compiler error at this method, you set the + * MaterialLaw property to something different than + * Opm::NullMaterialLaw. In this case, you have to overload the + * matererialLaw() method in the derived class! + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + const MaterialLawParams & + materialLawParams(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { + static MaterialLawParams dummy; + return dummy; + } + + /*! + * \brief Returns the temperature \f$\mathrm{[K]}\f$ within a control volume. + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + Scalar temperature(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { return asImp_().temperature(); } + + /*! + * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem. + * + * This is not specific to the discretization. By default it just + * throws an exception so it must be overloaded by the problem if + * no energy equation is to be used. + */ + Scalar temperature() const + { throw std::logic_error("Not implemented:temperature() method not implemented by the actual problem"); } + + + /*! + * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$. + * + * \param context Reference to the object which represents the + * current execution context. + * \param spaceIdx The local index of spatial entity defined by the context + * \param timeIdx The index used by the time discretization. + */ + template + const DimVector& gravity(const Context& context OPM_UNUSED, + unsigned spaceIdx OPM_UNUSED, + unsigned timeIdx OPM_UNUSED) const + { return asImp_().gravity(); } + + /*! + * \brief Returns the acceleration due to gravity \f$\mathrm{[m/s^2]}\f$. + * + * This method is used for problems where the gravitational + * acceleration does not depend on the spatial position. The + * default behaviour is that if the EnableGravity + * property is true, \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$ holds, + * else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$. + */ + const DimVector& gravity() const + { return gravity_; } + + /*! + * \brief Mark grid cells for refinement or coarsening + * + * \return The number of elements marked for refinement or coarsening. + */ + unsigned markForGridAdaptation() + { + typedef Opm::MathToolbox Toolbox; + + unsigned numMarked = 0; + ElementContext elemCtx( this->simulator() ); + auto gridView = this->simulator().vanguard().gridView(); + auto& grid = this->simulator().vanguard().grid(); + auto elemIt = gridView.template begin(); + auto elemEndIt = gridView.template end(); + for (; elemIt != elemEndIt; ++elemIt) + { + const auto& element = *elemIt ; + elemCtx.updateAll( element ); + + // HACK: this should better be part of an AdaptionCriterion class + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + Scalar minSat = 1e100 ; + Scalar maxSat = -1e100; + size_t nDofs = elemCtx.numDof(/*timeIdx=*/0); + for (unsigned dofIdx = 0; dofIdx < nDofs; ++dofIdx) + { + const auto& intQuant = elemCtx.intensiveQuantities( dofIdx, /*timeIdx=*/0 ); + minSat = std::min(minSat, + Toolbox::value(intQuant.fluidState().saturation(phaseIdx))); + maxSat = std::max(maxSat, + Toolbox::value(intQuant.fluidState().saturation(phaseIdx))); + } + + const Scalar indicator = + (maxSat - minSat)/(std::max(0.01, maxSat+minSat)/2); + if( indicator > 0.2 && element.level() < 2 ) { + grid.mark( 1, element ); + ++ numMarked; + } + else if ( indicator < 0.025 ) { + grid.mark( -1, element ); + ++ numMarked; + } + else + { + grid.mark( 0, element ); + } + } + } + + // get global sum so that every proc is on the same page + numMarked = this->simulator().vanguard().grid().comm().sum( numMarked ); + + return numMarked; + } + + // \} + +protected: + /*! + * \brief Converts a Scalar value to an isotropic Tensor + * + * This is convenient e.g. for specifying intrinsic permebilities: + * \code{.cpp} + * auto permTensor = this->toDimMatrix_(1e-12); + * \endcode + * + * \param val The scalar value which should be expressed as a tensor + */ + DimMatrix toDimMatrix_(Scalar val) const + { + DimMatrix ret(0.0); + for (unsigned i = 0; i < DimMatrix::rows; ++i) + ret[i][i] = val; + return ret; + } + + DimVector gravity_; + +private: + //! Returns the implementation of the problem (i.e. static polymorphism) + Implementation& asImp_() + { return *static_cast(this); } + //! \copydoc asImp_() + const Implementation& asImp_() const + { return *static_cast(this); } + + void init_() + { + gravity_ = 0.0; + if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) + gravity_[dimWorld-1] = -9.81; + } +}; + +} // namespace Opm + +#endif diff --git a/opm/models/common/multiphasebaseproperties.hh b/opm/models/common/multiphasebaseproperties.hh new file mode 100644 index 000000000..2cd115831 --- /dev/null +++ b/opm/models/common/multiphasebaseproperties.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 + * \ingroup MultiPhaseBaseModel + * + * \brief Defines the common properties required by the porous medium + * multi-phase models. + */ +#ifndef EWOMS_MULTI_PHASE_BASE_PROPERTIES_HH +#define EWOMS_MULTI_PHASE_BASE_PROPERTIES_HH + +#include +#include +#include + +BEGIN_PROPERTIES + +//! The splice to be used for the spatial discretization +NEW_PROP_TAG(SpatialDiscretizationSplice); +//! Number of fluid phases in the system +NEW_PROP_TAG(NumPhases); +//! Number of chemical species in the system +NEW_PROP_TAG(NumComponents); +//! Enumerations used by the model +NEW_PROP_TAG(Indices); +//! The material law which ought to be used (extracted from the spatial parameters) +NEW_PROP_TAG(MaterialLaw); +//! The context material law (extracted from the spatial parameters) +NEW_PROP_TAG(MaterialLawParams); +//! The material law for the energy stored in the solid matrix +NEW_PROP_TAG(SolidEnergyLaw); +//! The parameters of the material law for energy storage of the solid +NEW_PROP_TAG(SolidEnergyLawParams); +//! The material law for thermal conduction +NEW_PROP_TAG(ThermalConductionLaw); +//! The parameters of the material law for thermal conduction +NEW_PROP_TAG(ThermalConductionLawParams); +//!The fluid systems including the information about the phases +NEW_PROP_TAG(FluidSystem); +//! Specifies the relation used for velocity +NEW_PROP_TAG(FluxModule); + +//! Returns whether gravity is considered in the problem +NEW_PROP_TAG(EnableGravity); + +END_PROPERTIES + +#endif diff --git a/opm/models/common/quantitycallbacks.hh b/opm/models/common/quantitycallbacks.hh new file mode 100644 index 000000000..a1f2361c2 --- /dev/null +++ b/opm/models/common/quantitycallbacks.hh @@ -0,0 +1,482 @@ +// -*- 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 + * + * \brief This method contains all callback classes for quantities + * that are required by some extensive quantities + */ +#ifndef EWOMS_QUANTITY_CALLBACKS_HH +#define EWOMS_QUANTITY_CALLBACKS_HH + +#include + +#include +#include + +#include +#include + +namespace Opm { +/*! + * \ingroup Discretization + * + * \brief Callback class for temperature. + */ +template +class TemperatureCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(std::declval().fluidState()) IQFluidState; + typedef decltype(std::declval().temperature(0)) ResultRawType; + +public: + typedef typename std::remove_const::type>::type ResultType; + typedef typename Opm::MathToolbox::ValueType ResultValueType; + + TemperatureCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + {} + + /*! + * \brief Return the temperature given the index of a degree of freedom within an + * element context. + * + * In this context, we assume that thermal equilibrium applies, i.e. that the + * temperature of all phases is equal. + */ + ResultType operator()(unsigned dofIdx) const + { return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().temperature(/*phaseIdx=*/0); } + +private: + const ElementContext& elemCtx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for a phase pressure. + */ +template +class PressureCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(std::declval().fluidState()) IQFluidState; + typedef decltype(std::declval().pressure(0)) ResultRawType; + +public: + typedef typename std::remove_const::type>::type ResultType; + typedef typename Opm::MathToolbox::ValueType ResultValueType; + + PressureCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + { Opm::Valgrind::SetUndefined(phaseIdx_); } + + PressureCallback(const ElementContext& elemCtx, unsigned phaseIdx) + : elemCtx_(elemCtx) + , phaseIdx_(static_cast(phaseIdx)) + {} + + /*! + * \brief Set the index of the fluid phase for which the pressure + * should be returned. + */ + void setPhaseIndex(unsigned phaseIdx) + { phaseIdx_ = static_cast(phaseIdx); } + + /*! + * \brief Return the pressure of the specified phase given the index of a degree of + * freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { + Opm::Valgrind::CheckDefined(phaseIdx_); + return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().pressure(phaseIdx_); + } + +private: + const ElementContext& elemCtx_; + unsigned short phaseIdx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for a phase pressure. + */ +template +class BoundaryPressureCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar; + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(std::declval().fluidState()) IQRawFluidState; + typedef typename std::remove_const::type>::type IQFluidState; + typedef typename IQFluidState::Scalar IQScalar; + typedef Opm::MathToolbox Toolbox; + +public: + typedef IQScalar ResultType; + + BoundaryPressureCallback(const ElementContext& elemCtx, const FluidState& boundaryFs) + : elemCtx_(elemCtx) + , boundaryFs_(boundaryFs) + { Opm::Valgrind::SetUndefined(phaseIdx_); } + + BoundaryPressureCallback(const ElementContext& elemCtx, + const FluidState& boundaryFs, + unsigned phaseIdx) + : elemCtx_(elemCtx) + , boundaryFs_(boundaryFs) + , phaseIdx_(static_cast(phaseIdx)) + {} + + /*! + * \brief Set the index of the fluid phase for which the pressure + * should be returned. + */ + void setPhaseIndex(unsigned phaseIdx) + { phaseIdx_ = static_cast(phaseIdx); } + + /*! + * \brief Return the pressure of a phase given the index of a + * degree of freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { + Opm::Valgrind::CheckDefined(phaseIdx_); + return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().pressure(phaseIdx_); + } + + IQScalar boundaryValue() const + { + Opm::Valgrind::CheckDefined(phaseIdx_); + return boundaryFs_.pressure(phaseIdx_); + } + +private: + const ElementContext& elemCtx_; + const FluidState& boundaryFs_; + unsigned short phaseIdx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for the density of a phase. + */ +template +class DensityCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(std::declval().fluidState()) IQFluidState; + typedef decltype(std::declval().density(0)) ResultRawType; + +public: + typedef typename std::remove_const::type>::type ResultType; + typedef typename Opm::MathToolbox::ValueType ResultValueType; + + DensityCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + { Opm::Valgrind::SetUndefined(phaseIdx_); } + + DensityCallback(const ElementContext& elemCtx, unsigned phaseIdx) + : elemCtx_(elemCtx) + , phaseIdx_(static_cast(phaseIdx)) + {} + + /*! + * \brief Set the index of the fluid phase for which the density + * should be returned. + */ + void setPhaseIndex(unsigned phaseIdx) + { phaseIdx_ = static_cast(phaseIdx); } + + /*! + * \brief Return the density of a phase given the index of a + * degree of freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { + Opm::Valgrind::CheckDefined(phaseIdx_); + return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().density(phaseIdx_); + } + +private: + const ElementContext& elemCtx_; + unsigned short phaseIdx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for the molar density of a phase. + */ +template +class MolarDensityCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(std::declval().fluidState()) IQFluidState; + +public: + typedef decltype(std::declval().molarDensity(0)) ResultType; + typedef typename Opm::MathToolbox::ValueType ResultValueType; + + MolarDensityCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + { Opm::Valgrind::SetUndefined(phaseIdx_); } + + MolarDensityCallback(const ElementContext& elemCtx, unsigned phaseIdx) + : elemCtx_(elemCtx) + , phaseIdx_(static_cast(phaseIdx)) + {} + + /*! + * \brief Set the index of the fluid phase for which the molar + * density should be returned. + */ + void setPhaseIndex(unsigned phaseIdx) + { phaseIdx_ = static_cast(phaseIdx); } + + /*! + * \brief Return the molar density of a phase given the index of a + * degree of freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { + Opm::Valgrind::CheckDefined(phaseIdx_); + return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().molarDensity(phaseIdx_); + } + +private: + const ElementContext& elemCtx_; + unsigned short phaseIdx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for the viscosity of a phase. + */ +template +class ViscosityCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(std::declval().fluidState()) IQFluidState; + typedef decltype(std::declval().viscosity(0)) ResultRawType; + +public: + typedef typename std::remove_const::type>::type ResultType; + typedef typename Opm::MathToolbox::ValueType ResultValueType; + + ViscosityCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + { Opm::Valgrind::SetUndefined(phaseIdx_); } + + ViscosityCallback(const ElementContext& elemCtx, unsigned phaseIdx) + : elemCtx_(elemCtx) + , phaseIdx_(static_cast(phaseIdx)) + {} + + /*! + * \brief Set the index of the fluid phase for which the viscosity + * should be returned. + */ + void setPhaseIndex(unsigned phaseIdx) + { phaseIdx_ = static_cast(phaseIdx); } + + /*! + * \brief Return the viscosity of a phase given the index of a + * degree of freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { + Opm::Valgrind::CheckDefined(phaseIdx_); + return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().viscosity(phaseIdx_); + } + +private: + const ElementContext& elemCtx_; + unsigned short phaseIdx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for the velocity of a phase at the center of a DOF. + */ +template +class VelocityCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView; + + typedef decltype(IntensiveQuantities().velocityCenter()) ResultRawType; + + enum { dim = GridView::dimensionworld }; + +public: + typedef typename std::remove_const::type>::type ResultType; + typedef typename ResultType::field_type ResultFieldType; + typedef typename Opm::MathToolbox::ValueType ResultFieldValueType; + + VelocityCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + {} + + /*! + * \brief Return the velocity of a phase given the index of a + * degree of freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).velocityCenter(); } + +private: + const ElementContext& elemCtx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for the velocity of a phase at the center of a DOF. + */ +template +class VelocityComponentCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(IntensiveQuantities().velocityCenter()[0]) ResultRawType; + +public: + typedef typename std::remove_const::type>::type ResultType; + typedef typename Opm::MathToolbox::ValueType ResultValueType; + + VelocityComponentCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + { Opm::Valgrind::SetUndefined(dimIdx_); } + + VelocityComponentCallback(const ElementContext& elemCtx, unsigned dimIdx) + : elemCtx_(elemCtx) + , dimIdx_(dimIdx) + {} + + /*! + * \brief Set the index of the component of the velocity + * which should be returned. + */ + void setDimIndex(unsigned dimIdx) + { dimIdx_ = dimIdx; } + + /*! + * \brief Return the velocity of a phase given the index of a + * degree of freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { + Opm::Valgrind::CheckDefined(dimIdx_); + return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).velocityCenter()[dimIdx_]; + } + +private: + const ElementContext& elemCtx_; + unsigned dimIdx_; +}; + +/*! + * \ingroup Discretization + * + * \brief Callback class for a mole fraction of a component in a phase. + */ +template +class MoleFractionCallback +{ + typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext; + typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities; + + typedef decltype(std::declval().fluidState()) IQFluidState; + typedef decltype(std::declval().moleFraction(0, 0)) ResultRawType; + +public: + typedef typename std::remove_const::type>::type ResultType; + typedef typename Opm::MathToolbox::ValueType ResultValueType; + + MoleFractionCallback(const ElementContext& elemCtx) + : elemCtx_(elemCtx) + { + Opm::Valgrind::SetUndefined(phaseIdx_); + Opm::Valgrind::SetUndefined(compIdx_); + } + + MoleFractionCallback(const ElementContext& elemCtx, unsigned phaseIdx, unsigned compIdx) + : elemCtx_(elemCtx) + , phaseIdx_(static_cast(phaseIdx)) + , compIdx_(static_cast(compIdx)) + {} + + /*! + * \brief Set the index of the fluid phase for which a mole fraction should be + * returned. + */ + void setPhaseIndex(unsigned phaseIdx) + { phaseIdx_ = static_cast(phaseIdx); } + + /*! + * \brief Set the index of the component for which the mole fraction should be + * returned. + */ + void setComponentIndex(unsigned compIdx) + { compIdx_ = static_cast(compIdx); } + + /*! + * \brief Return the mole fraction of a component in a phase given the index of a + * degree of freedom within an element context. + */ + ResultType operator()(unsigned dofIdx) const + { + Opm::Valgrind::CheckDefined(phaseIdx_); + Opm::Valgrind::CheckDefined(compIdx_); + return elemCtx_.intensiveQuantities(dofIdx, /*timeIdx=*/0).fluidState().moleFraction(phaseIdx_, compIdx_); + } + +private: + const ElementContext& elemCtx_; + unsigned short phaseIdx_; + unsigned short compIdx_; +}; + +} // namespace Opm + +#endif