105 void update(
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
107 ParentType::update(elemCtx, dofIdx, timeIdx);
108 EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
110 const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
111 const auto& problem = elemCtx.problem();
113 const Scalar flashTolerance = Parameters::Get<Parameters::FlashTolerance<Scalar>>();
114 const int flashVerbosity = Parameters::Get<Parameters::FlashVerbosity>();
115 const std::string flashTwoPhaseMethod = Parameters::Get<Parameters::FlashTwoPhaseMethod>();
118 ComponentVector z(0.);
120 Evaluation lastZ = 1.0;
121 for (
unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
122 z[compIdx] = priVars.makeEvaluation(z0Idx + compIdx, timeIdx);
125 z[numComponents - 1] = lastZ;
127 Evaluation sumz = 0.0;
128 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
129 z[compIdx] = Opm::max(z[compIdx], 1e-8);
135 Evaluation p = priVars.makeEvaluation(pressure0Idx, timeIdx);
136 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
137 fluidState_.setPressure(phaseIdx, p);
140 const auto *hint = elemCtx.thermodynamicHint(dofIdx, timeIdx);
142 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
143 const Evaluation& Ktmp = hint->fluidState().K(compIdx);
144 fluidState_.setKvalue(compIdx, Ktmp);
146 const Evaluation& Ltmp = hint->fluidState().L();
147 fluidState_.setLvalue(Ltmp);
149 else if (timeIdx == 0 && elemCtx.thermodynamicHint(dofIdx, 1)) {
151 const auto& hint2 = elemCtx.thermodynamicHint(dofIdx, 1);
152 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
153 const Evaluation& Ktmp = hint2->fluidState().K(compIdx);
154 fluidState_.setKvalue(compIdx, Ktmp);
156 const Evaluation& Ltmp = hint2->fluidState().L();
157 fluidState_.setLvalue(Ltmp);
160 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
161 const Evaluation Ktmp = fluidState_.wilsonK_(compIdx);
162 fluidState_.setKvalue(compIdx, Ktmp);
164 const Evaluation& Ltmp = -1.0;
165 fluidState_.setLvalue(Ltmp);
171 if (flashVerbosity >= 1) {
172 const int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
173 std::cout <<
" updating the intensive quantities for Cell " << spatialIdx << std::endl;
175 FlashSolver::solve(fluidState_, z, flashTwoPhaseMethod, flashTolerance, flashVerbosity);
177 if (flashVerbosity >= 5) {
179 const int spatialIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
180 std::cout <<
" \n After flash solve for cell " << spatialIdx << std::endl;
181 ComponentVector x, y;
182 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
183 x[comp_idx] = fluidState_.moleFraction(FluidSystem::oilPhaseIdx, comp_idx);
184 y[comp_idx] = fluidState_.moleFraction(FluidSystem::gasPhaseIdx, comp_idx);
186 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
187 std::cout <<
" x for component: " << comp_idx <<
" is:" << std::endl;
188 std::cout << x[comp_idx] << std::endl;
190 std::cout <<
" y for component: " << comp_idx <<
"is:" << std::endl;
191 std::cout << y[comp_idx] << std::endl;
193 const Evaluation& L = fluidState_.L();
194 std::cout <<
" L is:" << std::endl;
195 std::cout << L << std::endl;
200 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
201 paramCache.updatePhase(fluidState_, FluidSystem::oilPhaseIdx);
203 const Scalar R = Opm::Constants<Scalar>::R;
204 Evaluation Z_L = (paramCache.molarVolume(FluidSystem::oilPhaseIdx) * fluidState_.pressure(FluidSystem::oilPhaseIdx) )/
205 (R * fluidState_.temperature(FluidSystem::oilPhaseIdx));
206 paramCache.updatePhase(fluidState_, FluidSystem::gasPhaseIdx);
207 Evaluation Z_V = (paramCache.molarVolume(FluidSystem::gasPhaseIdx) * fluidState_.pressure(FluidSystem::gasPhaseIdx) )/
208 (R * fluidState_.temperature(FluidSystem::gasPhaseIdx));
213 Evaluation L = fluidState_.L();
214 Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
215 Evaluation Sg = Opm::max(1 - So, 0.0);
216 Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg);
220 fluidState_.setSaturation(0, So);
221 fluidState_.setSaturation(1, Sg);
223 fluidState_.setCompressFactor(0, Z_L);
224 fluidState_.setCompressFactor(1, Z_V);
227 if (flashVerbosity >= 5) {
228 std::cout <<
"So = " << So <<std::endl;
229 std::cout <<
"Sg = " << Sg <<std::endl;
233 if (flashVerbosity >= 5) {
234 std::cout <<
"So = " << So <<std::endl;
235 std::cout <<
"Sg = " << Sg <<std::endl;
236 std::cout <<
"Z_L = " << Z_L <<std::endl;
237 std::cout <<
"Z_V = " << Z_V <<std::endl;
243 const MaterialLawParams& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
246 MaterialLaw::relativePermeabilities(relativePermeability_,
247 materialParams, fluidState_);
248 Opm::Valgrind::CheckDefined(relativePermeability_);
251 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
252 paramCache.updatePhase(fluidState_, phaseIdx);
254 const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
256 fluidState_.setViscosity(phaseIdx, mu);
258 mobility_[phaseIdx] = relativePermeability_[phaseIdx] / mu;
259 Opm::Valgrind::CheckDefined(mobility_[phaseIdx]);
261 const Evaluation& rho = FluidSystem::density(fluidState_, paramCache, phaseIdx);
262 fluidState_.setDensity(phaseIdx, rho);
270 porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
271 Opm::Valgrind::CheckDefined(porosity_);
274 intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
277 FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
280 EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
283 DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);