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_();
129 const SolutionVector& uLastIter)
134 const int localSwitched = numPriVarsSwitched_;
135 MPI_Allreduce(&localSwitched,
136 &numPriVarsSwitched_,
143 this->simulator_.model().newtonMethod().endIterMsg()
144 <<
", num switched=" << numPriVarsSwitched_;
146 ParentType::endIteration_(uCurrentIter, uLastIter);
151 const SolutionVector& currentSolution,
152 const GlobalEqVector& solutionUpdate,
153 const GlobalEqVector& currentResidual)
155 const auto& comm = this->simulator_.gridView().comm();
159 ParentType::update_(nextSolution,
168 succeeded = comm.min(succeeded);
171 throw NumericalProblem(
"A process did not succeed in adapting the primary variables");
174 numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
177 template <
class DofIndices>
179 const SolutionVector& currentSolution,
180 const GlobalEqVector& solutionUpdate,
181 const GlobalEqVector& currentResidual,
182 const DofIndices& dofIndices)
184 const auto zero = 0.0 * solutionUpdate[0];
185 for (
auto dofIdx : dofIndices) {
186 if (solutionUpdate[dofIdx] == zero) {
190 nextSolution[dofIdx],
191 currentSolution[dofIdx],
192 solutionUpdate[dofIdx],
193 currentResidual[dofIdx]);
202 PrimaryVariables& nextValue,
203 const PrimaryVariables& currentValue,
204 const EqVector& update,
205 const EqVector& currentResidual)
207 static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
208 static constexpr bool enableExtbo = Indices::zFractionIdx >= 0;
209 static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
210 static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
211 static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
212 static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0;
213 static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
214 static constexpr bool enableMICP = Indices::microbialConcentrationIdx >= 0;
216 currentValue.checkDefined();
217 Valgrind::CheckDefined(update);
218 Valgrind::CheckDefined(currentResidual);
221 Scalar deltaSw = 0.0;
222 Scalar deltaSo = 0.0;
223 Scalar deltaSg = 0.0;
224 Scalar deltaSs = 0.0;
226 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
228 deltaSw = update[Indices::waterSwitchIdx];
231 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
233 deltaSg = update[Indices::compositionSwitchIdx];
236 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
237 deltaSs = update[Indices::solventSaturationIdx];
242 Scalar maxSatDelta = std::max(std::abs(deltaSg), std::abs(deltaSo));
243 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSw));
244 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSs));
248 Scalar satAlpha = 1.0;
249 if (maxSatDelta > bparams_.dsMax_) {
250 satAlpha = bparams_.dsMax_ / maxSatDelta;
253 for (
int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
260 Scalar delta = update[pvIdx];
263 if (pvIdx == Indices::pressureSwitchIdx) {
264 if (std::abs(delta) > bparams_.dpMaxRel_ * currentValue[pvIdx]) {
265 delta =
signum(delta) * bparams_.dpMaxRel_ * currentValue[pvIdx];
269 else if (pvIdx == Indices::waterSwitchIdx)
270 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
275 if (delta > currentValue[ Indices::waterSwitchIdx]) {
276 delta = currentValue[ Indices::waterSwitchIdx];
279 else if (pvIdx == Indices::compositionSwitchIdx) {
285 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
290 if (delta > currentValue[Indices::compositionSwitchIdx]) {
291 delta = currentValue[Indices::compositionSwitchIdx];
295 else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
297 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
302 if (delta > currentValue[Indices::solventSaturationIdx]) {
303 delta = currentValue[Indices::solventSaturationIdx];
307 else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
309 const auto& curr = currentValue[Indices::zFractionIdx];
310 delta = std::clamp(delta, curr - Scalar{1.0}, curr);
312 else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
313 const double sign = delta >= 0. ? 1. : -1.;
316 const Scalar maxMolarWeightChange = 100.0;
317 delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
320 else if (enableEnergy && pvIdx == Indices::temperatureIdx) {
321 const double sign = delta >= 0. ? 1. : -1.;
322 delta = sign * std::min(std::abs(delta), bparams_.maxTempChange_);
324 else if (enableBrine && pvIdx == Indices::saltConcentrationIdx &&
325 enableSaltPrecipitation &&
326 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
328 const Scalar maxSaltSaturationChange = 0.1;
329 const Scalar sign = delta >= 0. ? 1. : -1.;
330 delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
334 nextValue[pvIdx] = currentValue[pvIdx] - delta;
337 if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
338 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
339 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
344 if (enableExtbo && pvIdx == Indices::zFractionIdx) {
345 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
349 if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx) {
350 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
353 if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
354 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
355 const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
356 if (polymerConcentration < 1.e-10) {
357 nextValue[pvIdx] = 0.0;
362 if (enableFoam && pvIdx == Indices::foamConcentrationIdx) {
363 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
366 if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
368 if (!enableSaltPrecipitation ||
369 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs)
371 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
374 if (enableSaltPrecipitation &&
375 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
377 nextValue[pvIdx] = std::min(nextValue[pvIdx], Scalar{1.0-1.e-8});
382 if (enableEnergy && pvIdx == Indices::temperatureIdx) {
383 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.tempMin_, bparams_.tempMax_);
386 if (pvIdx == Indices::pressureSwitchIdx) {
387 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_);
395 if constexpr (enableMICP) {
396 if (pvIdx == Indices::microbialConcentrationIdx) {
397 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
399 if (pvIdx == Indices::oxygenConcentrationIdx) {
400 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
402 if (pvIdx == Indices::ureaConcentrationIdx) {
403 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
405 if (pvIdx == Indices::biofilmConcentrationIdx) {
406 nextValue[pvIdx] = std::clamp(nextValue[pvIdx],
408 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
410 if (pvIdx == Indices::calciteConcentrationIdx) {
411 nextValue[pvIdx] = std::clamp(nextValue[pvIdx],
413 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
421 if (wasSwitched_[globalDofIdx]) {
422 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
424 bparams_.waterSaturationMax_,
425 bparams_.waterOnlyThreshold_,
426 bparams_.priVarOscilationThreshold_);
429 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
431 bparams_.waterSaturationMax_,
432 bparams_.waterOnlyThreshold_);
435 if (wasSwitched_[globalDofIdx]) {
436 ++numPriVarsSwitched_;
438 if (bparams_.projectSaturations_) {
439 nextValue.chopAndNormalizeSaturations();
442 nextValue.checkDefined();
446 int numPriVarsSwitched_{};
448 BlackoilNewtonParams<Scalar> bparams_{};
452 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:150
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:128
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:201
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:178
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.