28#ifndef OPM_FLASH_INTENSIVE_QUANTITIES_HH
29#define OPM_FLASH_INTENSIVE_QUANTITIES_HH
37#include <opm/material/fluidstates/CompositionalFluidState.hpp>
38#include <opm/material/common/Valgrind.hpp>
40#include <dune/common/fvector.hh>
41#include <dune/common/fmatrix.hh>
51template <
class TypeTag>
52class FlashIntensiveQuantities
53 :
public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
54 ,
public DiffusionIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>() >
55 ,
public EnergyIntensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableEnergy>() >
56 ,
public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
58 using ParentType = GetPropType<TypeTag, Properties::DiscIntensiveQuantities>;
60 using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
61 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
62 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
63 using Indices = GetPropType<TypeTag, Properties::Indices>;
64 using FluxModule = GetPropType<TypeTag, Properties::FluxModule>;
65 using GridView = GetPropType<TypeTag, Properties::GridView>;
66 using ThreadManager = GetPropType<TypeTag, Properties::ThreadManager>;
69 enum { z0Idx = Indices::z0Idx };
70 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
71 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
72 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
73 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
74 enum { dimWorld = GridView::dimensionworld };
75 enum { pressure0Idx = Indices::pressure0Idx };
77 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
78 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
79 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
80 using FlashSolver = GetPropType<TypeTag, Properties::FlashSolver>;
82 using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
83 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
85 using FluxIntensiveQuantities =
typename FluxModule::FluxIntensiveQuantities;
91 using FluidState = Opm::CompositionalFluidState<Evaluation, FluidSystem, enableEnergy>;
102 void update(
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
104 ParentType::update(elemCtx, dofIdx, timeIdx);
105 EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
107 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
108 const auto& problem = elemCtx.problem();
110 const Scalar flashTolerance = Parameters::get<TypeTag, Properties::FlashTolerance>();
111 const int flashVerbosity = Parameters::get<TypeTag, Properties::FlashVerbosity>();
112 const std::string flashTwoPhaseMethod = Parameters::get<TypeTag, Properties::FlashTwoPhaseMethod>();
115 ComponentVector z(0.);
117 Evaluation lastZ = 1.0;
118 for (
unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
119 z[compIdx] = priVars.makeEvaluation(z0Idx + compIdx, timeIdx);
122 z[numComponents - 1] = lastZ;
124 Evaluation sumz = 0.0;
125 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
126 z[compIdx] = Opm::max(z[compIdx], 1e-8);
132 Evaluation p = priVars.makeEvaluation(pressure0Idx, timeIdx);
133 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
134 fluidState_.setPressure(phaseIdx, p);
137 const auto *hint = elemCtx.thermodynamicHint(dofIdx, timeIdx);
139 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
140 const Evaluation& Ktmp = hint->fluidState().K(compIdx);
141 fluidState_.setKvalue(compIdx, Ktmp);
143 const Evaluation& Ltmp = hint->fluidState().L();
144 fluidState_.setLvalue(Ltmp);
146 else if (timeIdx == 0 && elemCtx.thermodynamicHint(dofIdx, 1)) {
148 const auto& hint2 = elemCtx.thermodynamicHint(dofIdx, 1);
149 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
150 const Evaluation& Ktmp = hint2->fluidState().K(compIdx);
151 fluidState_.setKvalue(compIdx, Ktmp);
153 const Evaluation& Ltmp = hint2->fluidState().L();
154 fluidState_.setLvalue(Ltmp);
157 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
158 const Evaluation Ktmp = fluidState_.wilsonK_(compIdx);
159 fluidState_.setKvalue(compIdx, Ktmp);
161 const Evaluation& Ltmp = -1.0;
162 fluidState_.setLvalue(Ltmp);
168 if (flashVerbosity >= 1) {
169 const int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
170 std::cout <<
" updating the intensive quantities for Cell " << spatialIdx << std::endl;
172 FlashSolver::solve(fluidState_, z, flashTwoPhaseMethod, flashTolerance, flashVerbosity);
174 if (flashVerbosity >= 5) {
176 const int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
177 std::cout <<
" \n After flash solve for cell " << spatialIdx << std::endl;
178 ComponentVector x, y;
179 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
180 x[comp_idx] = fluidState_.moleFraction(FluidSystem::oilPhaseIdx, comp_idx);
181 y[comp_idx] = fluidState_.moleFraction(FluidSystem::gasPhaseIdx, comp_idx);
183 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
184 std::cout <<
" x for component: " << comp_idx <<
" is:" << std::endl;
185 std::cout << x[comp_idx] << std::endl;
187 std::cout <<
" y for component: " << comp_idx <<
"is:" << std::endl;
188 std::cout << y[comp_idx] << std::endl;
190 const Evaluation& L = fluidState_.L();
191 std::cout <<
" L is:" << std::endl;
192 std::cout << L << std::endl;
197 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
198 paramCache.updatePhase(fluidState_, FluidSystem::oilPhaseIdx);
200 const Scalar R = Opm::Constants<Scalar>::R;
201 Evaluation Z_L = (paramCache.molarVolume(FluidSystem::oilPhaseIdx) * fluidState_.pressure(FluidSystem::oilPhaseIdx) )/
202 (R * fluidState_.temperature(FluidSystem::oilPhaseIdx));
203 paramCache.updatePhase(fluidState_, FluidSystem::gasPhaseIdx);
204 Evaluation Z_V = (paramCache.molarVolume(FluidSystem::gasPhaseIdx) * fluidState_.pressure(FluidSystem::gasPhaseIdx) )/
205 (R * fluidState_.temperature(FluidSystem::gasPhaseIdx));
210 Evaluation L = fluidState_.L();
211 Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
212 Evaluation Sg = Opm::max(1 - So, 0.0);
213 Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg);
217 fluidState_.setSaturation(0, So);
218 fluidState_.setSaturation(1, Sg);
220 fluidState_.setCompressFactor(0, Z_L);
221 fluidState_.setCompressFactor(1, Z_V);
224 if (flashVerbosity >= 5) {
225 std::cout <<
"So = " << So <<std::endl;
226 std::cout <<
"Sg = " << Sg <<std::endl;
230 if (flashVerbosity >= 5) {
231 std::cout <<
"So = " << So <<std::endl;
232 std::cout <<
"Sg = " << Sg <<std::endl;
233 std::cout <<
"Z_L = " << Z_L <<std::endl;
234 std::cout <<
"Z_V = " << Z_V <<std::endl;
240 const MaterialLawParams& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
243 MaterialLaw::relativePermeabilities(relativePermeability_,
244 materialParams, fluidState_);
245 Opm::Valgrind::CheckDefined(relativePermeability_);
248 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
249 paramCache.updatePhase(fluidState_, phaseIdx);
251 const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
253 fluidState_.setViscosity(phaseIdx, mu);
255 mobility_[phaseIdx] = relativePermeability_[phaseIdx] / mu;
256 Opm::Valgrind::CheckDefined(mobility_[phaseIdx]);
258 const Evaluation& rho = FluidSystem::density(fluidState_, paramCache, phaseIdx);
259 fluidState_.setDensity(phaseIdx, rho);
267 porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
268 Opm::Valgrind::CheckDefined(porosity_);
271 intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
274 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
277 EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
280 DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
287 {
return fluidState_; }
293 {
return intrinsicPerm_; }
299 {
return relativePermeability_[phaseIdx]; }
304 const Evaluation&
mobility(
unsigned phaseIdx)
const
306 return mobility_[phaseIdx];
313 {
return porosity_; }
316 DimMatrix intrinsicPerm_;
318 Evaluation porosity_;
319 std::array<Evaluation,numPhases> relativePermeability_;
320 std::array<Evaluation,numPhases> mobility_;
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes.
Definition: diffusionmodule.hh:140
Provides the volumetric quantities required for the energy equation.
Definition: energymodule.hh:530
Contains the intensive quantities of the flash-based compositional multi-phase model.
Definition: flash/flashintensivequantities.hh:57
Opm::CompositionalFluidState< Evaluation, FluidSystem, enableEnergy > FluidState
The type of the object returned by the fluidState() method.
Definition: flash/flashintensivequantities.hh:90
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: ptflash/flashintensivequantities.hh:304
const DimMatrix & intrinsicPermeability() const
Returns the intrinsic permeability tensor a degree of freedom.
Definition: ptflash/flashintensivequantities.hh:292
FlashIntensiveQuantities(const FlashIntensiveQuantities &other)=default
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: ptflash/flashintensivequantities.hh:286
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition: ptflash/flashintensivequantities.hh:312
const Evaluation & relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition: ptflash/flashintensivequantities.hh:298
FlashIntensiveQuantities()=default
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: ptflash/flashintensivequantities.hh:102
FlashIntensiveQuantities & operator=(const FlashIntensiveQuantities &other)=default
Classes required for molecular diffusion.
Contains the classes required to consider energy as a conservation quantity in a multi-phase module.
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:242
Declares the properties required by the compositional multi-phase model based on flash calculations.