28#ifndef EWOMS_BLACK_OIL_NEWTON_METHOD_HH
29#define EWOMS_BLACK_OIL_NEWTON_METHOD_HH
33#include <opm/common/Exceptions.hpp>
41template <
class TypeTag,
class MyTypeTag>
42struct DiscNewtonMethod;
44template<
class TypeTag,
class MyTypeTag>
46template<
class TypeTag,
class MyTypeTag>
48template<
class TypeTag,
class MyTypeTag>
50template<
class TypeTag,
class MyTypeTag>
52template<
class TypeTag,
class MyTypeTag>
54template<
class TypeTag,
class MyTypeTag>
56template<
class TypeTag,
class MyTypeTag>
58template<
class TypeTag,
class MyTypeTag>
60template<
class TypeTag,
class MyTypeTag>
62template<
class TypeTag,
class MyTypeTag>
64template<
class TypeTag,
class MyTypeTag>
66template<
class TypeTag>
70 static constexpr type value = 0.3;
72template<
class TypeTag>
76 static constexpr type value = 0.2;
78template<
class TypeTag>
82 static constexpr type value = 1e-5;
84template<
class TypeTag>
86template<
class TypeTag>
90 static constexpr type value = 5;
92template<
class TypeTag>
96 static constexpr type value = 1e9;
98template<
class TypeTag>
102 static constexpr type value = 0.0;
104template<
class TypeTag>
108 static constexpr type value = 1e99;
110template<
class TypeTag>
114 static constexpr type value = -1e99;
116template<
class TypeTag>
120 static constexpr type value = 1.0;
122template<
class TypeTag>
126 static constexpr type value = 1.0;
137template <
class TypeTag>
152 static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
153 static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
158 priVarOscilationThreshold_ = Parameters::get<TypeTag, Properties::PriVarOscilationThreshold>();
159 dpMaxRel_ = Parameters::get<TypeTag, Properties::DpMaxRel>();
160 dsMax_ = Parameters::get<TypeTag, Properties::DsMax>();
161 projectSaturations_ = Parameters::get<TypeTag, Properties::ProjectSaturations>();
162 maxTempChange_ = Parameters::get<TypeTag, Properties::MaxTemperatureChange>();
163 tempMax_ = Parameters::get<TypeTag, Properties::TemperatureMax>();
164 tempMin_ = Parameters::get<TypeTag, Properties::TemperatureMin>();
165 pressMax_ = Parameters::get<TypeTag, Properties::PressureMax>();
166 pressMin_ = Parameters::get<TypeTag, Properties::PressureMin>();
167 waterSaturationMax_ = Parameters::get<TypeTag, Properties::MaximumWaterSaturation>();
168 waterOnlyThreshold_ = Parameters::get<TypeTag, Properties::WaterOnlyThreshold>();
176 ParentType::finishInit();
178 wasSwitched_.resize(this->model().numTotalDof());
179 std::fill(wasSwitched_.begin(), wasSwitched_.end(),
false);
187 ParentType::registerParameters();
189 Parameters::registerParam<TypeTag, Properties::DpMaxRel>
190 (
"Maximum relative change of pressure in a single iteration");
191 Parameters::registerParam<TypeTag, Properties::DsMax>
192 (
"Maximum absolute change of any saturation in a single iteration");
193 Parameters::registerParam<TypeTag, Properties::PriVarOscilationThreshold>
194 (
"The threshold value for the primary variable switching conditions "
195 "after its meaning has switched to hinder oscilations");
196 Parameters::registerParam<TypeTag, Properties::ProjectSaturations>
197 (
"Option for doing saturation projection");
198 Parameters::registerParam<TypeTag, Properties::MaxTemperatureChange>
199 (
"Maximum absolute change of temperature in a single iteration");
200 Parameters::registerParam<TypeTag, Properties::TemperatureMax>
201 (
"Maximum absolute temperature");
202 Parameters::registerParam<TypeTag, Properties::TemperatureMin>
203 (
"Minimum absolute temperature");
204 Parameters::registerParam<TypeTag, Properties::PressureMax>
205 (
"Maximum absolute pressure");
206 Parameters::registerParam<TypeTag, Properties::PressureMin>
207 (
"Minimum absolute pressure");
208 Parameters::registerParam<TypeTag, Properties::MaximumWaterSaturation>
209 (
"Maximum water saturation");
210 Parameters::registerParam<TypeTag, Properties::WaterOnlyThreshold>
211 (
"Cells with water saturation above or equal is considered one-phase water only");
219 {
return numPriVarsSwitched_; }
230 numPriVarsSwitched_ = 0;
231 ParentType::beginIteration_();
238 const SolutionVector& uLastIter)
243 int localSwitched = numPriVarsSwitched_;
244 MPI_Allreduce(&localSwitched,
245 &numPriVarsSwitched_,
252 this->simulator_.model().newtonMethod().endIterMsg()
253 <<
", num switched=" << numPriVarsSwitched_;
255 ParentType::endIteration_(uCurrentIter, uLastIter);
260 const SolutionVector& currentSolution,
261 const GlobalEqVector& solutionUpdate,
262 const GlobalEqVector& currentResidual)
264 const auto& comm = this->simulator_.gridView().comm();
268 ParentType::update_(nextSolution,
277 succeeded = comm.min(succeeded);
280 throw NumericalProblem(
"A process did not succeed in adapting the primary variables");
282 numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
285 template <
class DofIndices>
287 const SolutionVector& currentSolution,
288 const GlobalEqVector& solutionUpdate,
289 const GlobalEqVector& currentResidual,
290 const DofIndices& dofIndices)
292 const auto zero = 0.0 * solutionUpdate[0];
293 for (
auto dofIdx : dofIndices) {
294 if (solutionUpdate[dofIdx] == zero) {
298 nextSolution[dofIdx],
299 currentSolution[dofIdx],
300 solutionUpdate[dofIdx],
301 currentResidual[dofIdx]);
310 PrimaryVariables& nextValue,
311 const PrimaryVariables& currentValue,
312 const EqVector& update,
313 const EqVector& currentResidual)
315 static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
316 static constexpr bool enableExtbo = Indices::zFractionIdx >= 0;
317 static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
318 static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
319 static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
320 static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0;
321 static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
322 static constexpr bool enableMICP = Indices::microbialConcentrationIdx >= 0;
324 currentValue.checkDefined();
325 Valgrind::CheckDefined(update);
326 Valgrind::CheckDefined(currentResidual);
329 Scalar deltaSw = 0.0;
330 Scalar deltaSo = 0.0;
331 Scalar deltaSg = 0.0;
332 Scalar deltaSs = 0.0;
334 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
336 deltaSw = update[Indices::waterSwitchIdx];
339 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
341 deltaSg = update[Indices::compositionSwitchIdx];
344 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
345 deltaSs = update[Indices::solventSaturationIdx];
350 Scalar maxSatDelta = std::max(std::abs(deltaSg), std::abs(deltaSo));
351 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSw));
352 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSs));
356 Scalar satAlpha = 1.0;
357 if (maxSatDelta > dsMax_)
358 satAlpha = dsMax_/maxSatDelta;
360 for (
int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
367 Scalar delta = update[pvIdx];
370 if (pvIdx == Indices::pressureSwitchIdx) {
371 if (std::abs(delta) > dpMaxRel_*currentValue[pvIdx])
372 delta =
signum(delta)*dpMaxRel_*currentValue[pvIdx];
375 else if (pvIdx == Indices::waterSwitchIdx)
376 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
380 if (delta > currentValue[ Indices::waterSwitchIdx])
381 delta = currentValue[ Indices::waterSwitchIdx];
383 else if (pvIdx == Indices::compositionSwitchIdx) {
389 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
393 if (delta > currentValue[Indices::compositionSwitchIdx])
394 delta = currentValue[Indices::compositionSwitchIdx];
397 else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
399 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
403 if (delta > currentValue[Indices::solventSaturationIdx])
404 delta = currentValue[Indices::solventSaturationIdx];
407 else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
409 const auto& curr = currentValue[Indices::zFractionIdx];
410 delta = std::clamp(delta, curr - 1.0, curr);
412 else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
413 const double sign = delta >= 0. ? 1. : -1.;
416 const double maxMolarWeightChange = 100.0;
417 delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
420 else if (enableEnergy && pvIdx == Indices::temperatureIdx) {
421 const double sign = delta >= 0. ? 1. : -1.;
422 delta = sign * std::min(std::abs(delta), maxTempChange_);
424 else if (enableBrine && pvIdx == Indices::saltConcentrationIdx &&
425 enableSaltPrecipitation &&
426 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp) {
427 const double maxSaltSaturationChange = 0.1;
428 const double sign = delta >= 0. ? 1. : -1.;
429 delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
433 nextValue[pvIdx] = currentValue[pvIdx] - delta;
436 if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
437 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss )
438 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
442 if (enableExtbo && pvIdx == Indices::zFractionIdx)
443 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
446 if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx)
447 nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
449 if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
450 nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
451 const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
452 if (polymerConcentration < 1.e-10)
453 nextValue[pvIdx] = 0.0;
457 if (enableFoam && pvIdx == Indices::foamConcentrationIdx)
458 nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
460 if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
462 if (!enableSaltPrecipitation || (enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs))
463 nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
465 if ((enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp))
466 nextValue[pvIdx] = std::min(nextValue[pvIdx], 1.0-1.e-8);
470 if (enableEnergy && pvIdx == Indices::temperatureIdx)
471 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], tempMin_, tempMax_);
473 if (pvIdx == Indices::pressureSwitchIdx) {
474 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], pressMin_, pressMax_);
483 if (enableMICP && pvIdx == Indices::microbialConcentrationIdx)
485 if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx)
487 if (enableMICP && pvIdx == Indices::ureaConcentrationIdx)
489 if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx)
491 if (enableMICP && pvIdx == Indices::calciteConcentrationIdx)
498 if (wasSwitched_[globalDofIdx])
499 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx, waterSaturationMax_, waterOnlyThreshold_, priVarOscilationThreshold_);
501 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(), globalDofIdx, waterSaturationMax_, waterOnlyThreshold_);
503 if (wasSwitched_[globalDofIdx])
504 ++ numPriVarsSwitched_;
505 if(projectSaturations_){
506 nextValue.chopAndNormalizeSaturations();
509 nextValue.checkDefined();
513 int numPriVarsSwitched_;
515 Scalar priVarOscilationThreshold_;
516 Scalar waterSaturationMax_;
517 Scalar waterOnlyThreshold_;
521 bool projectSaturations_;
522 Scalar maxTempChange_;
530 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:56
static const Scalar densityBiofilm()
Definition: blackoilmicpmodules.hh:351
static const Scalar toleranceBeforeClogging()
Definition: blackoilmicpmodules.hh:426
static const Scalar maximumUreaConcentration()
Definition: blackoilmicpmodules.hh:396
static const std::vector< Scalar > phi()
Definition: blackoilmicpmodules.hh:436
static const Scalar maximumOxygenConcentration()
Definition: blackoilmicpmodules.hh:391
A newton solver which is specific to the black oil model.
Definition: blackoilnewtonmethod.hh:139
void beginIteration_()
Indicates the beginning of a Newton iteration.
Definition: blackoilnewtonmethod.hh:228
static void registerParameters()
Register all run-time parameters for the immiscible model.
Definition: blackoilnewtonmethod.hh:185
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual)
Definition: blackoilnewtonmethod.hh:259
void finishInit()
Finialize the construction of the object.
Definition: blackoilnewtonmethod.hh:174
void endIteration_(SolutionVector &uCurrentIter, const SolutionVector &uLastIter)
Indicates that one Newton iteration was finished.
Definition: blackoilnewtonmethod.hh:237
unsigned numPriVarsSwitched() const
Returns the number of degrees of freedom for which the interpretation has changed for the most recent...
Definition: blackoilnewtonmethod.hh:218
BlackOilNewtonMethod(Simulator &simulator)
Definition: blackoilnewtonmethod.hh:156
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables ¤tValue, const EqVector &update, const EqVector ¤tResidual)
Update a single primary variables object.
Definition: blackoilnewtonmethod.hh:309
friend ParentType
Definition: blackoilnewtonmethod.hh:223
void update_(SolutionVector &nextSolution, const SolutionVector ¤tSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector ¤tResidual, const DofIndices &dofIndices)
Definition: blackoilnewtonmethod.hh:286
The multi-dimensional Newton method.
Definition: newtonmethod.hh:113
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:242
int signum(Scalar val)
Template function which returns the sign of a floating point value.
Definition: signum.hh:40
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:69
Definition: blackoilnewtonmethod.hh:45
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:75
Definition: blackoilnewtonmethod.hh:47
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:89
Definition: blackoilnewtonmethod.hh:53
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:119
Definition: blackoilnewtonmethod.hh:63
Specifies the type of the actual Newton method.
Definition: newtonmethodproperties.hh:32
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:107
Definition: blackoilnewtonmethod.hh:59
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:113
Definition: blackoilnewtonmethod.hh:61
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:81
Definition: blackoilnewtonmethod.hh:49
Definition: blackoilnewtonmethod.hh:51
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:95
Definition: blackoilnewtonmethod.hh:55
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:101
Definition: blackoilnewtonmethod.hh:57
a tag to mark properties as undefined
Definition: propertysystem.hh:40
GetPropType< TypeTag, Scalar > type
Definition: blackoilnewtonmethod.hh:125
Definition: blackoilnewtonmethod.hh:65