195 PrimaryVariables& nextValue,
196 const PrimaryVariables& currentValue,
197 const EqVector& update,
198 const EqVector& currentResidual)
200 static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
201 static constexpr bool enableExtbo = Indices::zFractionIdx >= 0;
202 static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
203 static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
204 static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
205 static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0;
206 static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
207 static constexpr bool enableMICP = Indices::microbialConcentrationIdx >= 0;
209 currentValue.checkDefined();
210 Valgrind::CheckDefined(update);
211 Valgrind::CheckDefined(currentResidual);
214 Scalar deltaSw = 0.0;
215 Scalar deltaSo = 0.0;
216 Scalar deltaSg = 0.0;
217 Scalar deltaSs = 0.0;
219 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
221 deltaSw = update[Indices::waterSwitchIdx];
224 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
226 deltaSg = update[Indices::compositionSwitchIdx];
229 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
230 deltaSs = update[Indices::solventSaturationIdx];
235 Scalar maxSatDelta = std::max(std::abs(deltaSg), std::abs(deltaSo));
236 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSw));
237 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSs));
241 Scalar satAlpha = 1.0;
242 if (maxSatDelta > bparams_.dsMax_) {
243 satAlpha = bparams_.dsMax_ / maxSatDelta;
246 for (
int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
253 Scalar delta = update[pvIdx];
256 if (pvIdx == Indices::pressureSwitchIdx) {
257 if (std::abs(delta) > bparams_.dpMaxRel_ * currentValue[pvIdx]) {
258 delta =
signum(delta) * bparams_.dpMaxRel_ * currentValue[pvIdx];
262 else if (pvIdx == Indices::waterSwitchIdx)
263 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
268 if (delta > currentValue[ Indices::waterSwitchIdx]) {
269 delta = currentValue[ Indices::waterSwitchIdx];
272 else if (pvIdx == Indices::compositionSwitchIdx) {
278 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
283 if (delta > currentValue[Indices::compositionSwitchIdx])
284 delta = currentValue[Indices::compositionSwitchIdx];
287 else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
289 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
293 if (delta > currentValue[Indices::solventSaturationIdx])
294 delta = currentValue[Indices::solventSaturationIdx];
297 else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
299 const auto& curr = currentValue[Indices::zFractionIdx];
300 delta = std::clamp(delta, curr - Scalar{1.0}, curr);
302 else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
303 const double sign = delta >= 0. ? 1. : -1.;
306 const Scalar maxMolarWeightChange = 100.0;
307 delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
310 else if (enableEnergy && pvIdx == Indices::temperatureIdx) {
311 const double sign = delta >= 0. ? 1. : -1.;
312 delta = sign * std::min(std::abs(delta), bparams_.maxTempChange_);
314 else if (enableBrine && pvIdx == Indices::saltConcentrationIdx &&
315 enableSaltPrecipitation &&
316 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
317 const Scalar maxSaltSaturationChange = 0.1;
318 const Scalar sign = delta >= 0. ? 1. : -1.;
319 delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
323 nextValue[pvIdx] = currentValue[pvIdx] - delta;
326 if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
327 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
328 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
333 if (enableExtbo && pvIdx == Indices::zFractionIdx) {
334 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
338 if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx) {
339 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
342 if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
343 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
344 const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
345 if (polymerConcentration < 1.e-10) {
346 nextValue[pvIdx] = 0.0;
351 if (enableFoam && pvIdx == Indices::foamConcentrationIdx) {
352 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
355 if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
357 if (!enableSaltPrecipitation || (enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs)) {
358 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
361 if ((enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)) {
362 nextValue[pvIdx] = std::min(nextValue[pvIdx], Scalar{1.0-1.e-8});
367 if (enableEnergy && pvIdx == Indices::temperatureIdx) {
368 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.tempMin_, bparams_.tempMax_);
371 if (pvIdx == Indices::pressureSwitchIdx) {
372 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_);
381 if (enableMICP && pvIdx == Indices::microbialConcentrationIdx) {
382 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::densityBiofilm());
384 if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx) {
385 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumOxygenConcentration());
387 if (enableMICP && pvIdx == Indices::ureaConcentrationIdx) {
388 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumUreaConcentration());
390 if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx) {
391 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
393 if (enableMICP && pvIdx == Indices::calciteConcentrationIdx) {
394 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
401 if (wasSwitched_[globalDofIdx]) {
402 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
404 bparams_.waterSaturationMax_,
405 bparams_.waterOnlyThreshold_,
406 bparams_.priVarOscilationThreshold_);
409 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
411 bparams_.waterSaturationMax_,
412 bparams_.waterOnlyThreshold_);
415 if (wasSwitched_[globalDofIdx]) {
416 ++numPriVarsSwitched_;
418 if (bparams_.projectSaturations_) {
419 nextValue.chopAndNormalizeSaturations();
422 nextValue.checkDefined();