28#ifndef OPM_BLACK_OIL_NEWTON_METHOD_HPP
29#define OPM_BLACK_OIL_NEWTON_METHOD_HPP
31#include <opm/common/Exceptions.hpp>
43template <
class TypeTag,
class MyTypeTag>
44struct DiscNewtonMethod;
55template <
class TypeTag>
70 static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
71 static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
84 ParentType::finishInit();
86 wasSwitched_.resize(this->model().numTotalDof(),
false);
94 ParentType::registerParameters();
103 {
return numPriVarsSwitched_; }
114 numPriVarsSwitched_ = 0;
115 ParentType::beginIteration_();
122 const SolutionVector& uLastIter)
127 int localSwitched = numPriVarsSwitched_;
128 MPI_Allreduce(&localSwitched,
129 &numPriVarsSwitched_,
136 this->simulator_.model().newtonMethod().endIterMsg()
137 <<
", num switched=" << numPriVarsSwitched_;
139 ParentType::endIteration_(uCurrentIter, uLastIter);
144 const SolutionVector& currentSolution,
145 const GlobalEqVector& solutionUpdate,
146 const GlobalEqVector& currentResidual)
148 const auto& comm = this->simulator_.gridView().comm();
152 ParentType::update_(nextSolution,
161 succeeded = comm.min(succeeded);
164 throw NumericalProblem(
"A process did not succeed in adapting the primary variables");
167 numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
170 template <
class DofIndices>
172 const SolutionVector& currentSolution,
173 const GlobalEqVector& solutionUpdate,
174 const GlobalEqVector& currentResidual,
175 const DofIndices& dofIndices)
177 const auto zero = 0.0 * solutionUpdate[0];
178 for (
auto dofIdx : dofIndices) {
179 if (solutionUpdate[dofIdx] == zero) {
183 nextSolution[dofIdx],
184 currentSolution[dofIdx],
185 solutionUpdate[dofIdx],
186 currentResidual[dofIdx]);
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) {
384 if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx) {
387 if (enableMICP && pvIdx == Indices::ureaConcentrationIdx) {
390 if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx) {
393 if (enableMICP && pvIdx == Indices::calciteConcentrationIdx) {
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();
426 int numPriVarsSwitched_{};
428 BlackoilNewtonParams<Scalar> bparams_{};
432 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:49
static const Scalar densityBiofilm()
Definition: blackoilmicpmodules.hh:263
static const Scalar toleranceBeforeClogging()
Definition: blackoilmicpmodules.hh:338
static const Scalar maximumUreaConcentration()
Definition: blackoilmicpmodules.hh:308
static const std::vector< Scalar > phi()
Definition: blackoilmicpmodules.hh:348
static const Scalar maximumOxygenConcentration()
Definition: blackoilmicpmodules.hh:303
A newton solver which is specific to the black oil model.
Definition: blackoilnewtonmethod.hpp:57
void beginIteration_()
Indicates the beginning of a Newton iteration.
Definition: blackoilnewtonmethod.hpp:112
static void registerParameters()
Register all run-time parameters for the blackoil newton method.
Definition: blackoilnewtonmethod.hpp:92
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual)
Definition: blackoilnewtonmethod.hpp:143
void finishInit()
Finialize the construction of the object.
Definition: blackoilnewtonmethod.hpp:82
void endIteration_(SolutionVector &uCurrentIter, const SolutionVector &uLastIter)
Indicates that one Newton iteration was finished.
Definition: blackoilnewtonmethod.hpp:121
unsigned numPriVarsSwitched() const
Returns the number of degrees of freedom for which the interpretation has changed for the most recent...
Definition: blackoilnewtonmethod.hpp:102
BlackOilNewtonMethod(Simulator &simulator)
Definition: blackoilnewtonmethod.hpp:74
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector ¤tResidual)
Update a single primary variables object.
Definition: blackoilnewtonmethod.hpp:194
friend ParentType
Definition: blackoilnewtonmethod.hpp:107
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual, const DofIndices &dofIndices)
Definition: blackoilnewtonmethod.hpp:171
The multi-dimensional Newton method.
Definition: newtonmethod.hh:92
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
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:235
int signum(Scalar val)
Template function which returns the sign of a floating point value.
Definition: signum.hh:40
static void registerParameters()
Registers the parameters in parameter system.