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
41#include <algorithm>
42#include <cmath>
43#include <vector>
44
45namespace Opm::Properties {
46
47template <class TypeTag, class MyTypeTag>
48struct DiscNewtonMethod;
49
50} // namespace Opm::Properties
51
52namespace Opm {
53
59template <class TypeTag>
60class BlackOilNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMethod>
61{
73
74 static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
75 static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
76
77public:
78 explicit BlackOilNewtonMethod(Simulator& simulator) : ParentType(simulator)
79 {
80 bparams_.read();
81 }
82
87 {
88 ParentType::finishInit();
89
90 wasSwitched_.resize(this->model().numTotalDof(), false);
91 }
92
96 static void registerParameters()
97 {
98 ParentType::registerParameters();
100 }
101
106 unsigned numPriVarsSwitched() const
107 { return numPriVarsSwitched_; }
108
109protected:
112
117 {
118 numPriVarsSwitched_ = 0;
119 ParentType::beginIteration_();
120 }
121
125 void endIteration_(SolutionVector& uCurrentIter,
126 const SolutionVector& uLastIter)
127 {
128#if HAVE_MPI
129 // in the MPI enabled case we need to add up the number of DOF
130 // for which the interpretation changed over all processes.
131 const int localSwitched = numPriVarsSwitched_;
132 MPI_Allreduce(&localSwitched,
133 &numPriVarsSwitched_,
134 /*num=*/1,
135 MPI_INT,
136 MPI_SUM,
137 MPI_COMM_WORLD);
138#endif // HAVE_MPI
139
140 this->simulator_.model().newtonMethod().endIterMsg()
141 << ", num switched=" << numPriVarsSwitched_;
142
143 ParentType::endIteration_(uCurrentIter, uLastIter);
144 }
145
146public:
147 void update_(SolutionVector& nextSolution,
148 const SolutionVector& currentSolution,
149 const GlobalEqVector& solutionUpdate,
150 const GlobalEqVector& currentResidual)
151 {
152 const auto& comm = this->simulator_.gridView().comm();
153
154 int succeeded;
155 try {
156 ParentType::update_(nextSolution,
157 currentSolution,
158 solutionUpdate,
159 currentResidual);
160 succeeded = 1;
161 }
162 catch (...) {
163 succeeded = 0;
164 }
165 succeeded = comm.min(succeeded);
166
167 if (!succeeded) {
168 throw NumericalProblem("A process did not succeed in adapting the primary variables");
169 }
170
171 numPriVarsSwitched_ = comm.sum(numPriVarsSwitched_);
172 }
173
174 template <class DofIndices>
175 void update_(SolutionVector& nextSolution,
176 const SolutionVector& currentSolution,
177 const GlobalEqVector& solutionUpdate,
178 const GlobalEqVector& currentResidual,
179 const DofIndices& dofIndices)
180 {
181 const auto zero = 0.0 * solutionUpdate[0];
182 for (auto dofIdx : dofIndices) {
183 if (solutionUpdate[dofIdx] == zero) {
184 continue;
185 }
187 nextSolution[dofIdx],
188 currentSolution[dofIdx],
189 solutionUpdate[dofIdx],
190 currentResidual[dofIdx]);
191 }
192 }
193
194protected:
198 void updatePrimaryVariables_(unsigned globalDofIdx,
199 PrimaryVariables& nextValue,
200 const PrimaryVariables& currentValue,
201 const EqVector& update,
202 const EqVector& currentResidual)
203 {
204 static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
205 static constexpr bool enableExtbo = Indices::zFractionIdx >= 0;
206 static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
207 static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
208 static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
209 static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0;
210 static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
211 static constexpr bool enableMICP = Indices::microbialConcentrationIdx >= 0;
212
213 currentValue.checkDefined();
214 Valgrind::CheckDefined(update);
215 Valgrind::CheckDefined(currentResidual);
216
217 // saturation delta for each phase
218 Scalar deltaSw = 0.0;
219 Scalar deltaSo = 0.0;
220 Scalar deltaSg = 0.0;
221 Scalar deltaSs = 0.0;
222
223 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw)
224 {
225 deltaSw = update[Indices::waterSwitchIdx];
226 deltaSo -= deltaSw;
227 }
228 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg)
229 {
230 deltaSg = update[Indices::compositionSwitchIdx];
231 deltaSo -= deltaSg;
232 }
233 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
234 deltaSs = update[Indices::solventSaturationIdx];
235 deltaSo -= deltaSs;
236 }
237
238 // maximum saturation delta
239 Scalar maxSatDelta = std::max(std::abs(deltaSg), std::abs(deltaSo));
240 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSw));
241 maxSatDelta = std::max(maxSatDelta, std::abs(deltaSs));
242
243 // scaling factor for saturation deltas to make sure that none of them exceeds
244 // the specified threshold value.
245 Scalar satAlpha = 1.0;
246 if (maxSatDelta > bparams_.dsMax_) {
247 satAlpha = bparams_.dsMax_ / maxSatDelta;
248 }
249
250 for (int pvIdx = 0; pvIdx < int(numEq); ++pvIdx) {
251 // calculate the update of the current primary variable. For the black-oil
252 // model we limit the pressure delta relative to the pressure's current
253 // absolute value (Default: 30%) and saturation deltas to an absolute change
254 // (Default: 20%). Further, we ensure that the R factors, solvent
255 // "saturation" and polymer concentration do not become negative after the
256 // update.
257 Scalar delta = update[pvIdx];
258
259 // limit pressure delta
260 if (pvIdx == Indices::pressureSwitchIdx) {
261 if (std::abs(delta) > bparams_.dpMaxRel_ * currentValue[pvIdx]) {
262 delta = signum(delta) * bparams_.dpMaxRel_ * currentValue[pvIdx];
263 }
264 }
265 // water saturation delta
266 else if (pvIdx == Indices::waterSwitchIdx)
267 if (currentValue.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
268 delta *= satAlpha;
269 }
270 else {
271 //Ensure Rvw and Rsw factor does not become negative
272 if (delta > currentValue[ Indices::waterSwitchIdx]) {
273 delta = currentValue[ Indices::waterSwitchIdx];
274 }
275 }
276 else if (pvIdx == Indices::compositionSwitchIdx) {
277 // the switching primary variable for composition is tricky because the
278 // "reasonable" value ranges it exhibits vary widely depending on its
279 // interpretation since it can represent Sg, Rs or Rv. For now, we only
280 // limit saturation deltas and ensure that the R factors do not become
281 // negative.
282 if (currentValue.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Sg) {
283 delta *= satAlpha;
284 }
285 else {
286 // Ensure Rv and Rs factor does not become negative
287 if (delta > currentValue[Indices::compositionSwitchIdx]) {
288 delta = currentValue[Indices::compositionSwitchIdx];
289 }
290 }
291 }
292 else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
293 // solvent saturation updates are also subject to the Appleyard chop
294 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
295 delta *= satAlpha;
296 }
297 else {
298 // Ensure Rssolw factor does not become negative
299 if (delta > currentValue[Indices::solventSaturationIdx]) {
300 delta = currentValue[Indices::solventSaturationIdx];
301 }
302 }
303 }
304 else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
305 // z fraction updates are also subject to the Appleyard chop
306 const auto& curr = currentValue[Indices::zFractionIdx]; // or currentValue[pvIdx] given the block condition
307 delta = std::clamp(delta, curr - Scalar{1.0}, curr);
308 }
309 else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
310 const double sign = delta >= 0. ? 1. : -1.;
311 // maximum change of polymer molecular weight, the unit is MDa.
312 // applying this limit to stabilize the simulation. The value itself is still experimental.
313 const Scalar maxMolarWeightChange = 100.0;
314 delta = sign * std::min(std::abs(delta), maxMolarWeightChange);
315 delta *= satAlpha;
316 }
317 else if (enableEnergy && pvIdx == Indices::temperatureIdx) {
318 const double sign = delta >= 0. ? 1. : -1.;
319 delta = sign * std::min(std::abs(delta), bparams_.maxTempChange_);
320 }
321 else if (enableBrine && pvIdx == Indices::saltConcentrationIdx &&
322 enableSaltPrecipitation &&
323 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
324 {
325 const Scalar maxSaltSaturationChange = 0.1;
326 const Scalar sign = delta >= 0. ? 1. : -1.;
327 delta = sign * std::min(std::abs(delta), maxSaltSaturationChange);
328 }
329
330 // do the actual update
331 nextValue[pvIdx] = currentValue[pvIdx] - delta;
332
333 // keep the solvent saturation between 0 and 1
334 if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
335 if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
336 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
337 }
338 }
339
340 // keep the z fraction between 0 and 1
341 if (enableExtbo && pvIdx == Indices::zFractionIdx) {
342 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], Scalar{0.0}), Scalar{1.0});
343 }
344
345 // keep the polymer concentration above 0
346 if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx) {
347 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
348 }
349
350 if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
351 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
352 const double polymerConcentration = nextValue[Indices::polymerConcentrationIdx];
353 if (polymerConcentration < 1.e-10) {
354 nextValue[pvIdx] = 0.0;
355 }
356 }
357
358 // keep the foam concentration above 0
359 if (enableFoam && pvIdx == Indices::foamConcentrationIdx) {
360 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
361 }
362
363 if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
364 // keep the salt concentration above 0
365 if (!enableSaltPrecipitation ||
366 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Cs)
367 {
368 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
369 }
370 // keep the salt saturation below upperlimit
371 if (enableSaltPrecipitation &&
372 currentValue.primaryVarsMeaningBrine() == PrimaryVariables::BrineMeaning::Sp)
373 {
374 nextValue[pvIdx] = std::min(nextValue[pvIdx], Scalar{1.0-1.e-8});
375 }
376 }
377
378 // keep the temperature within given values
379 if (enableEnergy && pvIdx == Indices::temperatureIdx) {
380 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.tempMin_, bparams_.tempMax_);
381 }
382
383 if (pvIdx == Indices::pressureSwitchIdx) {
384 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_);
385 }
386
387 // keep the values above 0
388 // for the biofilm and calcite, we set an upper limit equal to the initial porosity
389 // minus 1e-8. This prevents singularities (e.g., one of the calcite source term is
390 // evaluated at 1/(iniPoro - calcite)). The value 1e-8 is taken from the salt precipitation
391 // clapping above.
392 if constexpr (enableMICP) {
393 if (pvIdx == Indices::microbialConcentrationIdx) {
394 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
395 }
396 if (pvIdx == Indices::oxygenConcentrationIdx) {
397 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
398 }
399 if (pvIdx == Indices::ureaConcentrationIdx) {
400 nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
401 }
402 if (pvIdx == Indices::biofilmConcentrationIdx) {
403 nextValue[pvIdx] = std::clamp(nextValue[pvIdx],
404 Scalar{0.0},
405 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
406 }
407 if (pvIdx == Indices::calciteConcentrationIdx) {
408 nextValue[pvIdx] = std::clamp(nextValue[pvIdx],
409 Scalar{0.0},
410 this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
411 }
412 }
413 }
414
415 // switch the new primary variables to something which is physically meaningful.
416 // use a threshold value after a switch to make it harder to switch back
417 // immediately.
418 if (wasSwitched_[globalDofIdx]) {
419 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
420 globalDofIdx,
421 bparams_.waterSaturationMax_,
422 bparams_.waterOnlyThreshold_,
423 bparams_.priVarOscilationThreshold_);
424 }
425 else {
426 wasSwitched_[globalDofIdx] = nextValue.adaptPrimaryVariables(this->problem(),
427 globalDofIdx,
428 bparams_.waterSaturationMax_,
429 bparams_.waterOnlyThreshold_);
430 }
431
432 if (wasSwitched_[globalDofIdx]) {
433 ++numPriVarsSwitched_;
434 }
435 if (bparams_.projectSaturations_) {
436 nextValue.chopAndNormalizeSaturations();
437 }
438
439 nextValue.checkDefined();
440 }
441
442private:
443 int numPriVarsSwitched_{};
444
445 BlackoilNewtonParams<Scalar> bparams_{};
446
447 // keep track of cells where the primary variable meaning has changed
448 // to detect and hinder oscillations
449 std::vector<bool> wasSwitched_{};
450};
451
452} // namespace Opm
453
454#endif // OPM_BLACK_OIL_NEWTHON_METHOD_HPP
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 &currentSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector &currentResidual)
Definition: blackoilnewtonmethod.hpp:147
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:125
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 &currentValue, const EqVector &update, const EqVector &currentResidual)
Update a single primary variables object.
Definition: blackoilnewtonmethod.hpp:198
friend ParentType
Definition: blackoilnewtonmethod.hpp:111
void update_(SolutionVector &nextSolution, const SolutionVector &currentSolution, const GlobalEqVector &solutionUpdate, const GlobalEqVector &currentResidual, const DofIndices &dofIndices)
Definition: blackoilnewtonmethod.hpp:175
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.