My Project
Loading...
Searching...
No Matches
blackoilnewtonmethod.hpp
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_BLACK_OIL_NEWTON_METHOD_HPP
29#define OPM_BLACK_OIL_NEWTON_METHOD_HPP
30
31#include <opm/common/Exceptions.hpp>
32
36
38
40
41namespace Opm::Properties {
42
43template <class TypeTag, class MyTypeTag>
44struct DiscNewtonMethod;
45
46} // namespace Opm::Properties
47
48namespace Opm {
49
55template <class TypeTag>
56class BlackOilNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMethod>
57{
69
70 static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
71 static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
72
73public:
74 BlackOilNewtonMethod(Simulator& simulator) : ParentType(simulator)
75 {
76 bparams_.read();
77 }
78
83 {
84 ParentType::finishInit();
85
86 wasSwitched_.resize(this->model().numTotalDof(), false);
87 }
88
92 static void registerParameters()
93 {
94 ParentType::registerParameters();
96 }
97
102 unsigned numPriVarsSwitched() const
103 { return numPriVarsSwitched_; }
104
105protected:
107 friend ParentType;
108
113 {
114 numPriVarsSwitched_ = 0;
115 ParentType::beginIteration_();
116 }
117
121 void endIteration_(SolutionVector& uCurrentIter,
122 const SolutionVector& uLastIter)
123 {
124#if HAVE_MPI
125 // in the MPI enabled case we need to add up the number of DOF
126 // for which the interpretation changed over all processes.
127 int localSwitched = numPriVarsSwitched_;
129 &numPriVarsSwitched_,
130 /*num=*/1,
131 MPI_INT,
132 MPI_SUM,
134#endif // HAVE_MPI
135
136 this->simulator_.model().newtonMethod().endIterMsg()
137 << ", num switched=" << numPriVarsSwitched_;
138
139 ParentType::endIteration_(uCurrentIter, uLastIter);
140 }
141
142public:
143 void update_(SolutionVector& nextSolution,
144 const SolutionVector& currentSolution,
145 const GlobalEqVector& solutionUpdate,
146 const GlobalEqVector& currentResidual)
147 {
148 const auto& comm = this->simulator_.gridView().comm();
149
150 int succeeded;
151 try {
152 ParentType::update_(nextSolution,
156 succeeded = 1;
157 }
158 catch (...) {
159 succeeded = 0;
160 }
161 succeeded = comm.min(succeeded);
162
163 if (!succeeded) {
164 throw NumericalProblem("A process did not succeed in adapting the primary variables");
165 }
166
167 numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
168 }
169
170 template <class DofIndices>
171 void update_(SolutionVector& nextSolution,
172 const SolutionVector& currentSolution,
173 const GlobalEqVector& solutionUpdate,
174 const GlobalEqVector& currentResidual,
175 const DofIndices& dofIndices)
176 {
177 const auto zero = 0.0 * solutionUpdate[0];
178 for (auto dofIdx : dofIndices) {
179 if (solutionUpdate[dofIdx] == zero) {
180 continue;
181 }
187 }
188 }
189
190protected:
195 PrimaryVariables& nextValue,
196 const PrimaryVariables& currentValue,
197 const EqVector& update,
198 const EqVector& currentResidual)
199 {
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;
208
209 currentValue.checkDefined();
210 Valgrind::CheckDefined(update);
211 Valgrind::CheckDefined(currentResidual);
212
213 // saturation delta for each phase
214 Scalar deltaSw = 0.0;
215 Scalar deltaSo = 0.0;
216 Scalar deltaSg = 0.0;
217 Scalar deltaSs = 0.0;
218
219 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
220 {
221 deltaSw = update[Indices::waterSwitchIdx];
222 deltaSo -= deltaSw;
223 }
224 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
225 {
226 deltaSg = update[Indices::compositionSwitchIdx];
227 deltaSo -= deltaSg;
228 }
229 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
230 deltaSs = update[Indices::solventSaturationIdx];
231 deltaSo -= deltaSs;
232 }
233
234 // maximum saturation delta
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));
238
239 // scaling factor for saturation deltas to make sure that none of them exceeds
240 // the specified threshold value.
241 Scalar satAlpha = 1.0;
242 if (maxSatDelta > bparams_.dsMax_) {
243 satAlpha = bparams_.dsMax_ / maxSatDelta;
244 }
245
246 for (int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
247 // calculate the update of the current primary variable. For the black-oil
248 // model we limit the pressure delta relative to the pressure's current
249 // absolute value (Default: 30%) and saturation deltas to an absolute change
250 // (Default: 20%). Further, we ensure that the R factors, solvent
251 // "saturation" and polymer concentration do not become negative after the
252 // update.
253 Scalar delta = update[pvIdx];
254
255 // limit pressure delta
256 if (pvIdx == Indices::pressureSwitchIdx) {
257 if (std::abs(delta) > bparams_.dpMaxRel_ * currentValue[pvIdx]) {
258 delta = signum(delta) * bparams_.dpMaxRel_ * currentValue[pvIdx];
259 }
260 }
261 // water saturation delta
262 else if (pvIdx == Indices::waterSwitchIdx)
263 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
264 delta *= satAlpha;
265 }
266 else {
267 //Ensure Rvw and Rsw factor does not become negative
268 if (delta > currentValue[ Indices::waterSwitchIdx]) {
269 delta = currentValue[ Indices::waterSwitchIdx];
270 }
271 }
272 else if (pvIdx == Indices::compositionSwitchIdx) {
273 // the switching primary variable for composition is tricky because the
274 // "reasonable" value ranges it exhibits vary widely depending on its
275 // interpretation since it can represent Sg, Rs or Rv. For now, we only
276 // limit saturation deltas and ensure that the R factors do not become
277 // negative.
278 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
279 delta *= satAlpha;
280 }
281 else {
282 //Ensure Rv and Rs factor does not become negative
283 if (delta > currentValue[Indices::compositionSwitchIdx])
284 delta = currentValue[Indices::compositionSwitchIdx];
285 }
286 }
287 else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
288 // solvent saturation updates are also subject to the Appleyard chop
289 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
290 delta *= satAlpha;
291 } else {
292 //Ensure Rssolw factor does not become negative
293 if (delta > currentValue[Indices::solventSaturationIdx])
294 delta = currentValue[Indices::solventSaturationIdx];
295 }
296 }
297 else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
298 // z fraction updates are also subject to the Appleyard chop
299 const auto& curr = currentValue[Indices::zFractionIdx]; // or currentValue[pvIdx] given the block condition
300 delta = std::clamp(delta, curr - Scalar{1.0}, curr);
301 }
302 else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
303 const double sign = delta >= 0. ? 1. : -1.;
304 // maximum change of polymer molecular weight, the unit is MDa.
305 // applying this limit to stabilize the simulation. The value itself is still experimental.
306 const Scalar maxMolarWeightChange = 100.0;
307 delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
308 delta *= satAlpha;
309 }
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_);
313 }
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);
320 }
321
322 // do the actual update
324
325 // keep the solvent saturation between 0 and 1
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});
329 }
330 }
331
332 // keep the z fraction between 0 and 1
333 if (enableExtbo && pvIdx == Indices::zFractionIdx) {
334 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
335 }
336
337 // keep the polymer concentration above 0
338 if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx) {
339 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
340 }
341
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;
347 }
348 }
349
350 // keep the foam concentration above 0
351 if (enableFoam && pvIdx == Indices::foamConcentrationIdx) {
352 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
353 }
354
355 if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
356 // keep the salt concentration above 0
357 if (!enableSaltPrecipitation || (enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs)) {
358 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
359 }
360 // keep the salt saturation below upperlimit
361 if ((enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)) {
362 nextValue[pvIdx] = std::min(nextValue[pvIdx], Scalar{1.0-1.e-8});
363 }
364 }
365
366 // keep the temperature within given values
367 if (enableEnergy && pvIdx == Indices::temperatureIdx) {
368 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.tempMin_, bparams_.tempMax_);
369 }
370
371 if (pvIdx == Indices::pressureSwitchIdx) {
372 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_);
373 }
374
375
376 // Limit the variables to [0, cmax] values to improve the convergence.
377 // For the microorganisms we set this value equal to the biomass density value.
378 // For the oxygen and urea we set this value to the maximum injected
379 // concentration (the urea concentration has been scaled by 10). For
380 // the biofilm and calcite, we set this value equal to the porosity minus the clogging tolerance.
381 if (enableMICP && pvIdx == Indices::microbialConcentrationIdx) {
382 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::densityBiofilm());
383 }
384 if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx) {
385 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumOxygenConcentration());
386 }
387 if (enableMICP && pvIdx == Indices::ureaConcentrationIdx) {
388 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumUreaConcentration());
389 }
390 if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx) {
391 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
392 }
393 if (enableMICP && pvIdx == Indices::calciteConcentrationIdx) {
394 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
395 }
396 }
397
398 // switch the new primary variables to something which is physically meaningful.
399 // use a threshold value after a switch to make it harder to switch back
400 // immediately.
401 if (wasSwitched_[globalDofIdx]) {
402 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
404 bparams_.waterSaturationMax_,
405 bparams_.waterOnlyThreshold_,
406 bparams_.priVarOscilationThreshold_);
407 }
408 else {
409 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
411 bparams_.waterSaturationMax_,
412 bparams_.waterOnlyThreshold_);
413 }
414
415 if (wasSwitched_[globalDofIdx]) {
416 ++numPriVarsSwitched_;
417 }
418 if (bparams_.projectSaturations_) {
419 nextValue.chopAndNormalizeSaturations();
420 }
421
422 nextValue.checkDefined();
423 }
424
425private:
426 int numPriVarsSwitched_{};
427
429
430 // keep track of cells where the primary variable meaning has changed
431 // to detect and hinder oscillations
432 std::vector<bool> wasSwitched_{};
433};
434
435} // namespace Opm
436
437#endif // OPM_BLACK_OIL_NEWTHON_METHOD_HPP
Contains the classes required to extend the black-oil model by MICP.
A newton solver which is specific to the black oil model.
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
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 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
void updatePrimaryVariables_(unsigned globalDofIdx, PrimaryVariables &nextValue, const PrimaryVariables &currentValue, const EqVector &update, const EqVector &currentResidual)
Update a single primary variables object.
Definition blackoilnewtonmethod.hpp:194
The multi-dimensional Newton method.
Definition newtonmethod.hh:92
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
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
The multi-dimensional Newton method.
static void registerParameters()
Registers the parameters in parameter system.
Definition blackoilnewtonmethodparams.cpp:32