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 enableBioeffects = Indices::biofilmVolumeFractionIdx >= 0;
215 static constexpr bool enableMICP = Indices::enableMICP;
217 currentValue.checkDefined();
218 Valgrind::CheckDefined(update);
219 Valgrind::CheckDefined(currentResidual);
222 Scalar deltaSw = 0.0;
223 Scalar deltaSo = 0.0;
224 Scalar deltaSg = 0.0;
225 Scalar deltaSs = 0.0;
227 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
229 deltaSw = update[Indices::waterSwitchIdx];
232 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
234 deltaSg = update[Indices::compositionSwitchIdx];
237 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
238 deltaSs = update[Indices::solventSaturationIdx];
243 Scalar maxSatDelta = std::max(std::abs(deltaSg), std::abs(deltaSo));
244 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSw));
245 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSs));
249 Scalar satAlpha = 1.0;
250 if (maxSatDelta > bparams_.dsMax_) {
251 satAlpha = bparams_.dsMax_ / maxSatDelta;
254 for (
int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
261 Scalar delta = update[pvIdx];
264 if (pvIdx == Indices::pressureSwitchIdx) {
265 if (std::abs(delta) > bparams_.dpMaxRel_ * currentValue[pvIdx]) {
266 delta =
signum(delta) * bparams_.dpMaxRel_ * currentValue[pvIdx];
270 else if (pvIdx == Indices::waterSwitchIdx)
271 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
276 if (delta > currentValue[ Indices::waterSwitchIdx]) {
277 delta = currentValue[ Indices::waterSwitchIdx];
280 else if (pvIdx == Indices::compositionSwitchIdx) {
286 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
291 if (delta > currentValue[Indices::compositionSwitchIdx]) {
292 delta = currentValue[Indices::compositionSwitchIdx];
296 else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
298 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
303 if (delta > currentValue[Indices::solventSaturationIdx]) {
304 delta = currentValue[Indices::solventSaturationIdx];
308 else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
310 const auto& curr = currentValue[Indices::zFractionIdx];
311 delta = std::clamp(delta, curr - Scalar{1.0}, curr);
313 else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
314 const double sign = delta >= 0. ? 1. : -1.;
317 const Scalar maxMolarWeightChange = 100.0;
318 delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
321 else if (enableEnergy && pvIdx == Indices::temperatureIdx) {
322 const double sign = delta >= 0. ? 1. : -1.;
323 delta = sign * std::min(std::abs(delta), bparams_.maxTempChange_);
325 else if (enableBrine && pvIdx == Indices::saltConcentrationIdx &&
326 enableSaltPrecipitation &&
327 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
329 const Scalar maxSaltSaturationChange = 0.1;
330 const Scalar sign = delta >= 0. ? 1. : -1.;
331 delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
335 nextValue[pvIdx] = currentValue[pvIdx] - delta;
338 if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
339 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
340 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
345 if (enableExtbo && pvIdx == Indices::zFractionIdx) {
346 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
350 if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx) {
351 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
354 if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
355 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
356 const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
357 if (polymerConcentration < 1.e-10) {
358 nextValue[pvIdx] = 0.0;
363 if (enableFoam && pvIdx == Indices::foamConcentrationIdx) {
364 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
367 if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
369 if (!enableSaltPrecipitation ||
370 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs)
372 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
375 if (enableSaltPrecipitation &&
376 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
378 nextValue[pvIdx] = std::min(nextValue[pvIdx], Scalar{1.0-1.e-8});
383 if (enableEnergy && pvIdx == Indices::temperatureIdx) {
384 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.tempMin_, bparams_.tempMax_);
387 if (pvIdx == Indices::pressureSwitchIdx) {
388 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_);
396 if constexpr (enableBioeffects) {
397 if (pvIdx == Indices::microbialConcentrationIdx) {
398 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
400 if (pvIdx == Indices::biofilmVolumeFractionIdx) {
401 nextValue[pvIdx] = std::clamp(nextValue[pvIdx],
403 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
405 if constexpr (enableMICP) {
406 if (pvIdx == Indices::oxygenConcentrationIdx) {
407 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
409 if (pvIdx == Indices::ureaConcentrationIdx) {
410 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
412 if (pvIdx == Indices::calciteVolumeFractionIdx) {
413 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0},
414 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
423 if (wasSwitched_[globalDofIdx]) {
424 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
426 bparams_.waterSaturationMax_,
427 bparams_.waterOnlyThreshold_,
428 bparams_.priVarOscilationThreshold_);
431 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
433 bparams_.waterSaturationMax_,
434 bparams_.waterOnlyThreshold_);
437 if (wasSwitched_[globalDofIdx]) {
438 ++numPriVarsSwitched_;
440 if (bparams_.projectSaturations_) {
441 nextValue.chopAndNormalizeSaturations();
444 nextValue.checkDefined();
448 int numPriVarsSwitched_{};
450 BlackoilNewtonParams<Scalar> bparams_{};
454 std::vector<bool> wasSwitched_{};
Contains the classes required to extend the black-oil model by bioeffects.
Declares the properties required by the black oil model.
Contains the high level supplements required to extend the black oil model by bioeffects.
Definition: blackoilbioeffectsmodules.hh:93
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: blackoilbioeffectsmodules.hh:43
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.