28#ifndef OPM_BLACK_OIL_NEWTON_METHOD_HPP
29#define OPM_BLACK_OIL_NEWTON_METHOD_HPP
31#include <opm/common/Exceptions.hpp>
47template <
class TypeTag,
class MyTypeTag>
48struct DiscNewtonMethod;
59template <
class TypeTag>
74 static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
75 static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
88 ParentType::finishInit();
90 wasSwitched_.resize(this->model().numTotalDof(),
false);
98 ParentType::registerParameters();
107 {
return numPriVarsSwitched_; }
118 numPriVarsSwitched_ = 0;
119 ParentType::beginIteration_();
126 const SolutionVector& uLastIter)
131 const int localSwitched = numPriVarsSwitched_;
132 MPI_Allreduce(&localSwitched,
133 &numPriVarsSwitched_,
140 this->simulator_.model().newtonMethod().endIterMsg()
141 <<
", num switched=" << numPriVarsSwitched_;
143 ParentType::endIteration_(uCurrentIter, uLastIter);
148 const SolutionVector& currentSolution,
149 const GlobalEqVector& solutionUpdate,
150 const GlobalEqVector& currentResidual)
152 const auto& comm = this->simulator_.gridView().comm();
156 ParentType::update_(nextSolution,
165 succeeded = comm.min(succeeded);
168 throw NumericalProblem(
"A process did not succeed in adapting the primary variables");
171 numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
174 template <
class DofIndices>
176 const SolutionVector& currentSolution,
177 const GlobalEqVector& solutionUpdate,
178 const GlobalEqVector& currentResidual,
179 const DofIndices& dofIndices)
181 const auto zero = 0.0 * solutionUpdate[0];
182 for (
auto dofIdx : dofIndices) {
183 if (solutionUpdate[dofIdx] == zero) {
187 nextSolution[dofIdx],
188 currentSolution[dofIdx],
189 solutionUpdate[dofIdx],
190 currentResidual[dofIdx]);
199 PrimaryVariables& nextValue,
200 const PrimaryVariables& currentValue,
201 const EqVector& update,
202 const EqVector& currentResidual)
204 static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
205 static constexpr bool enableExtbo = Indices::zFractionIdx >= 0;
206 static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
207 static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
208 static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
209 static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0;
210 static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
211 static constexpr bool enableMICP = Indices::microbialConcentrationIdx >= 0;
213 currentValue.checkDefined();
214 Valgrind::CheckDefined(update);
215 Valgrind::CheckDefined(currentResidual);
218 Scalar deltaSw = 0.0;
219 Scalar deltaSo = 0.0;
220 Scalar deltaSg = 0.0;
221 Scalar deltaSs = 0.0;
223 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
225 deltaSw = update[Indices::waterSwitchIdx];
228 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
230 deltaSg = update[Indices::compositionSwitchIdx];
233 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
234 deltaSs = update[Indices::solventSaturationIdx];
239 Scalar maxSatDelta = std::max(std::abs(deltaSg), std::abs(deltaSo));
240 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSw));
241 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSs));
245 Scalar satAlpha = 1.0;
246 if (maxSatDelta > bparams_.dsMax_) {
247 satAlpha = bparams_.dsMax_ / maxSatDelta;
250 for (
int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
257 Scalar delta = update[pvIdx];
260 if (pvIdx == Indices::pressureSwitchIdx) {
261 if (std::abs(delta) > bparams_.dpMaxRel_ * currentValue[pvIdx]) {
262 delta =
signum(delta) * bparams_.dpMaxRel_ * currentValue[pvIdx];
266 else if (pvIdx == Indices::waterSwitchIdx)
267 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
272 if (delta > currentValue[ Indices::waterSwitchIdx]) {
273 delta = currentValue[ Indices::waterSwitchIdx];
276 else if (pvIdx == Indices::compositionSwitchIdx) {
282 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
287 if (delta > currentValue[Indices::compositionSwitchIdx]) {
288 delta = currentValue[Indices::compositionSwitchIdx];
292 else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
294 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
299 if (delta > currentValue[Indices::solventSaturationIdx]) {
300 delta = currentValue[Indices::solventSaturationIdx];
304 else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
306 const auto& curr = currentValue[Indices::zFractionIdx];
307 delta = std::clamp(delta, curr - Scalar{1.0}, curr);
309 else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
310 const double sign = delta >= 0. ? 1. : -1.;
313 const Scalar maxMolarWeightChange = 100.0;
314 delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
317 else if (enableEnergy && pvIdx == Indices::temperatureIdx) {
318 const double sign = delta >= 0. ? 1. : -1.;
319 delta = sign * std::min(std::abs(delta), bparams_.maxTempChange_);
321 else if (enableBrine && pvIdx == Indices::saltConcentrationIdx &&
322 enableSaltPrecipitation &&
323 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
325 const Scalar maxSaltSaturationChange = 0.1;
326 const Scalar sign = delta >= 0. ? 1. : -1.;
327 delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
331 nextValue[pvIdx] = currentValue[pvIdx] - delta;
334 if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
335 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
336 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
341 if (enableExtbo && pvIdx == Indices::zFractionIdx) {
342 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
346 if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx) {
347 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
350 if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
351 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
352 const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
353 if (polymerConcentration < 1.e-10) {
354 nextValue[pvIdx] = 0.0;
359 if (enableFoam && pvIdx == Indices::foamConcentrationIdx) {
360 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
363 if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
365 if (!enableSaltPrecipitation ||
366 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs)
368 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
371 if (enableSaltPrecipitation &&
372 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
374 nextValue[pvIdx] = std::min(nextValue[pvIdx], Scalar{1.0-1.e-8});
379 if (enableEnergy && pvIdx == Indices::temperatureIdx) {
380 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.tempMin_, bparams_.tempMax_);
383 if (pvIdx == Indices::pressureSwitchIdx) {
384 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_);
392 if constexpr (enableMICP) {
393 if (pvIdx == Indices::microbialConcentrationIdx) {
394 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
396 if (pvIdx == Indices::oxygenConcentrationIdx) {
397 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
399 if (pvIdx == Indices::ureaConcentrationIdx) {
400 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
402 if (pvIdx == Indices::biofilmConcentrationIdx) {
403 nextValue[pvIdx] = std::clamp(nextValue[pvIdx],
405 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
407 if (pvIdx == Indices::calciteConcentrationIdx) {
408 nextValue[pvIdx] = std::clamp(nextValue[pvIdx],
410 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
418 if (wasSwitched_[globalDofIdx]) {
419 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
421 bparams_.waterSaturationMax_,
422 bparams_.waterOnlyThreshold_,
423 bparams_.priVarOscilationThreshold_);
426 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
428 bparams_.waterSaturationMax_,
429 bparams_.waterOnlyThreshold_);
432 if (wasSwitched_[globalDofIdx]) {
433 ++numPriVarsSwitched_;
435 if (bparams_.projectSaturations_) {
436 nextValue.chopAndNormalizeSaturations();
439 nextValue.checkDefined();
443 int numPriVarsSwitched_{};
445 BlackoilNewtonParams<Scalar> bparams_{};
449 std::vector<bool> wasSwitched_{};
Contains the classes required to extend the black-oil model by MICP.
Declares the properties required by the black oil model.
Contains the high level supplements required to extend the black oil model by MICP.
Definition: blackoilmicpmodules.hh:54
A newton solver which is specific to the black oil model.
Definition: blackoilnewtonmethod.hpp:61
void beginIteration_()
Indicates the beginning of a Newton iteration.
Definition: blackoilnewtonmethod.hpp:116
static void registerParameters()
Register all run-time parameters for the blackoil newton method.
Definition: blackoilnewtonmethod.hpp:96
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual)
Definition: blackoilnewtonmethod.hpp:147
void finishInit()
Finialize the construction of the object.
Definition: blackoilnewtonmethod.hpp:86
void endIteration_(SolutionVector &uCurrentIter, const SolutionVector &uLastIter)
Indicates that one Newton iteration was finished.
Definition: blackoilnewtonmethod.hpp:125
unsigned numPriVarsSwitched() const
Returns the number of degrees of freedom for which the interpretation has changed for the most recent...
Definition: blackoilnewtonmethod.hpp:106
BlackOilNewtonMethod(Simulator &simulator)
Definition: blackoilnewtonmethod.hpp:78
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector ¤tResidual)
Update a single primary variables object.
Definition: blackoilnewtonmethod.hpp:198
friend ParentType
Definition: blackoilnewtonmethod.hpp:111
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual, const DofIndices &dofIndices)
Definition: blackoilnewtonmethod.hpp:175
The multi-dimensional Newton method.
Definition: newtonmethod.hh:100
Definition: blackoilmodel.hh:79
Definition: blackoilboundaryratevector.hh:39
constexpr int signum(Scalar val)
Template function which returns the sign of a floating point value.
Definition: signum.hh:41
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:233
static void registerParameters()
Registers the parameters in parameter system.