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