Renames residuals() to computeResidualNorms()

It took me quite some time to understand the computations done
e.g. during the detection of oscillations, where the stuff returned
by residuals() is used as a vector of doubles. It turns out that
residuals() actually returns the norm of the residuals. To clarify this
we rename residuals() to computeResidualNorms() and residuals to
residual_norms. Having my dare devil day today, I even try to document the
method. (This documented method might feel kind of lonely between the others,
now;). Hopefully this saves others some time.
This commit is contained in:
Markus Blatt 2015-01-27 11:57:36 +01:00
parent 7dfcf7b0c0
commit 5e94e2ab08
2 changed files with 17 additions and 11 deletions

View File

@ -272,7 +272,11 @@ namespace Opm {
double double
residualNorm() const; residualNorm() const;
std::vector<double> residuals() const; /// \brief Compute the residual norms of the mass balance for each phase,
/// the well flux, and the well equation.
/// \return a vector that contains for each phase the norm of the mass balance
/// and afterwards the norm of the residuum of the well flux and the well equation.
std::vector<double> computeResidualNorms() const;
ADB ADB
fluidViscosity(const int phase, fluidViscosity(const int phase,

View File

@ -281,7 +281,9 @@ namespace {
computeWellConnectionPressures(state, xw); computeWellConnectionPressures(state, xw);
} }
std::vector<std::vector<double>> residual_history; // For each iteration we store in a vector the norms of the residual of
// the mass balance for each active phase, the well flux and the well equations
std::vector<std::vector<double>> residual_norms_history;
assemble(pvdt, x, xw); assemble(pvdt, x, xw);
@ -289,7 +291,7 @@ namespace {
bool converged = false; bool converged = false;
double omega = 1.; double omega = 1.;
residual_history.push_back(residuals()); residual_norms_history.push_back(computeResidualNorms());
int it = 0; int it = 0;
converged = getConvergence(dt,it); converged = getConvergence(dt,it);
@ -308,7 +310,7 @@ namespace {
// store number of linear iterations used // store number of linear iterations used
linearIterations += linsolver_.iterations(); linearIterations += linsolver_.iterations();
detectNewtonOscillations(residual_history, it, relaxRelTol(), isOscillate, isStagnate); detectNewtonOscillations(residual_norms_history, it, relaxRelTol(), isOscillate, isStagnate);
if (isOscillate) { if (isOscillate) {
omega -= relaxIncrement(); omega -= relaxIncrement();
@ -322,7 +324,7 @@ namespace {
assemble(pvdt, x, xw); assemble(pvdt, x, xw);
residual_history.push_back(residuals()); residual_norms_history.push_back(computeResidualNorms());
// increase iteration counter // increase iteration counter
++it; ++it;
@ -1783,9 +1785,9 @@ namespace {
template<class T> template<class T>
std::vector<double> std::vector<double>
FullyImplicitBlackoilSolver<T>::residuals() const FullyImplicitBlackoilSolver<T>::computeResidualNorms() const
{ {
std::vector<double> residual; std::vector<double> residualNorms;
std::vector<ADB>::const_iterator massBalanceIt = residual_.material_balance_eq.begin(); std::vector<ADB>::const_iterator massBalanceIt = residual_.material_balance_eq.begin();
const std::vector<ADB>::const_iterator endMassBalanceIt = residual_.material_balance_eq.end(); const std::vector<ADB>::const_iterator endMassBalanceIt = residual_.material_balance_eq.end();
@ -1796,7 +1798,7 @@ namespace {
OPM_THROW(Opm::NumericalProblem, OPM_THROW(Opm::NumericalProblem,
"Encountered a non-finite residual"); "Encountered a non-finite residual");
} }
residual.push_back(massBalanceResid); residualNorms.push_back(massBalanceResid);
} }
// the following residuals are not used in the oscillation detection now // the following residuals are not used in the oscillation detection now
@ -1805,16 +1807,16 @@ namespace {
OPM_THROW(Opm::NumericalProblem, OPM_THROW(Opm::NumericalProblem,
"Encountered a non-finite residual"); "Encountered a non-finite residual");
} }
residual.push_back(wellFluxResid); residualNorms.push_back(wellFluxResid);
const double wellResid = infinityNorm( residual_.well_eq ); const double wellResid = infinityNorm( residual_.well_eq );
if (!std::isfinite(wellResid)) { if (!std::isfinite(wellResid)) {
OPM_THROW(Opm::NumericalProblem, OPM_THROW(Opm::NumericalProblem,
"Encountered a non-finite residual"); "Encountered a non-finite residual");
} }
residual.push_back(wellResid); residualNorms.push_back(wellResid);
return residual; return residualNorms;
} }
template<class T> template<class T>