OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2800 from atgeirr/hnil-standard_nonlinearsolve

Browse Source 
Replace root finder with opm-common implementation
This commit is contained in:
Atgeirr Flø Rasmussen 2020-09-26 17:11:24 +02:00 committed by GitHub
commit 1aa1c2d322
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
1 changed files with 12 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

76
ebos/equil/equilibrationhelpers.hh
View File

@ -35,7 +35,9 @@
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
#include <opm/common/utility/numeric/RootFinders.hpp>
#include <cmath>
#include <memory>
@ -825,45 +827,17 @@ double satFromPc(const MaterialLawManager& materialLawManager,
const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, targetPc);
double f0 = f(s0);
double f1 = f(s1);
if (f0 <= 0.0)
return s0;
else if (f1 >= 0.0)
return s1;
assert(f0 > 0 && f1 < 0);
const int maxIter = 1000;
const double tol = 1e-10;
// regula falsi with the "Pegasus" method to avoid stagnation
for (int iter = 0; iter < maxIter; ++ iter) {
// determine the pivot
double s = (s1*f0 - s0*f1)/(f0 - f1);
// adapt the interval
double fs = f(s);
if (fs == 0.0)
return s;
else if ((fs > 0.0) == (f0 > 0.0)) {
// update interval and reverse
s0 = s1;
f0 = f1;
}
else
// "Pegasus" method
f0 *= f1/(f1 + fs);
s1 = s;
f1 = fs;
// check for convergence
if (std::abs(s1 - s0) < tol)
return (s1 + s0)/2;
}
throw std::runtime_error("Could not find solution for PcEq = 0.0 after "+std::to_string(maxIter)
+" iterations.");
// should at least converge in 2 times bisection but some safety here:
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
int usedIterations = -1;
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
return root;
}
@ -940,38 +914,12 @@ double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
return s1;
assert(f0 > 0.0 && f1 < 0.0);
const int maxIter = 60;
const double tol = 1e-10;
// regula falsi with the "Pegasus" method to avoid stagnation
for (int iter = 0; iter < maxIter; ++ iter) {
// determine the pivot
double s = (s1*f0 - s0*f1)/(f0 - f1);
// adapt the interval
double fs = f(s);
if (fs == 0.0)
return s;
else if ((fs > 0.0) == (f0 > 0.0)) {
// update interval and reverse
s0 = s1;
f0 = f1;
}
else
// "Pegasus" method
f0 *= f1 / (f1 + fs);
s1 = s;
f1 = fs;
// check for convergence
if (std::abs(s1 - s0) < tol)
return (s1 + s0)/2;
}
throw std::runtime_error("Could not find solution for PcEqSum = 0.0 after "+std::to_string(maxIter)
+" iterations.");
// should at least converge in 2 times bisection but some safety here:
const int maxIter = -2*static_cast<int>(std::log2(tol)) + 10;
int usedIterations = -1;
const double root = RegulaFalsiBisection<ThrowOnError>::solve(f, s0, s1, maxIter, tol, usedIterations);
return root;
}
/// Compute saturation from depth. Used for constant capillary pressure function