Factors out the reduction part into a separate function.

Theses reductions are global and will require a different
behaviour in the parallel case.
This commit is contained in:
Markus Blatt 2015-01-20 14:21:45 +01:00
parent 0b81cd8f6f
commit ddc2b820a7
2 changed files with 43 additions and 16 deletions

View File

@ -28,6 +28,8 @@
#include <opm/autodiff/LinearisedBlackoilResidual.hpp> #include <opm/autodiff/LinearisedBlackoilResidual.hpp>
#include <opm/autodiff/NewtonIterationBlackoilInterface.hpp> #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
#include <array>
struct UnstructuredGrid; struct UnstructuredGrid;
struct Wells; struct Wells;
@ -351,6 +353,16 @@ namespace Opm {
/// residual mass balance (tol_cnv). /// residual mass balance (tol_cnv).
bool getConvergence(const double dt, const int iteration); bool getConvergence(const double dt, const int iteration);
/// Compute the reduction within the convergence check.
void
convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& tempV,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& R,
std::array<double,MaxNumPhases>& B_avg,
std::array<double,MaxNumPhases>& maxCoeff,
std::array<double,MaxNumPhases>& R_sum,
int nc) const;
void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history, void detectNewtonOscillations(const std::vector<std::vector<double>>& residual_history,
const int it, const double relaxRelTol, const int it, const double relaxRelTol,
bool& oscillate, bool& stagnate) const; bool& oscillate, bool& stagnate) const;

View File

@ -1885,6 +1885,30 @@ namespace {
return; return;
} }
template<class T>
void
FullyImplicitBlackoilSolver<T>::convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& tempV,
const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& R,
std::array<double,MaxNumPhases>& R_sum,
std::array<double,MaxNumPhases>& maxCoeff,
std::array<double,MaxNumPhases>& B_avg,
int nc) const
{
const Opm::PhaseUsage& pu = fluid_.phaseUsage();
// Do the global reductions
for(int idx=0; idx<MaxNumPhases; ++idx)
{
if (active_[idx]) {
B_avg[pu.phase_pos[idx]] = B.sum()/nc;
maxCoeff[pu.phase_pos[idx]]=tempV.col(pu.phase_pos[idx]).maxCoeff();
R_sum[pu.phase_pos[idx]] = R.col(pu.phase_pos[idx]).sum();
}else{
R_sum[pu.phase_pos[idx]] = B_avg[pu.phase_pos[idx]] = maxCoeff[pu.phase_pos[idx]] =0.;
}
}
}
template<class T> template<class T>
bool bool
FullyImplicitBlackoilSolver<T>::getConvergence(const double dt, const int iteration) FullyImplicitBlackoilSolver<T>::getConvergence(const double dt, const int iteration)
@ -1900,11 +1924,11 @@ namespace {
const std::vector<PhasePresence> cond = phaseCondition(); const std::vector<PhasePresence> cond = phaseCondition();
double CNV[MaxNumPhases] = {0., 0., 0.}; std::array<double,MaxNumPhases> CNV = {{0., 0., 0.}};
double R_sum[MaxNumPhases] = {0., 0., 0.}; std::array<double,MaxNumPhases> R_sum = {{0., 0., 0.}};
double B_avg[MaxNumPhases] = {0., 0., 0.}; std::array<double,MaxNumPhases> B_avg = {{0., 0., 0.}};
double maxCoeff[MaxNumPhases] = {0., 0., 0.}; std::array<double,MaxNumPhases> maxCoeff = {{0., 0., 0.}};
double mass_balance_residual[MaxNumPhases] = {0., 0., 0.}; std::array<double,MaxNumPhases> mass_balance_residual = {{0., 0., 0.}};
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B; Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B;
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R; Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R;
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV; Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV;
@ -1919,17 +1943,8 @@ namespace {
tempV.col(pu.phase_pos[idx]) = R.col(pu.phase_pos[idx]).abs()/pv; tempV.col(pu.phase_pos[idx]) = R.col(pu.phase_pos[idx]).abs()/pv;
} }
} }
// Do the global reductions
for(int idx=0; idx<MaxNumPhases; ++idx) convergenceReduction(B, tempV, R, B_avg, maxCoeff, R_sum, nc);
{
if (active_[idx]) {
B_avg[pu.phase_pos[idx]] = B.sum()/nc;
maxCoeff[pu.phase_pos[idx]]=tempV.col(pu.phase_pos[idx]).maxCoeff();
R_sum[pu.phase_pos[idx]] = R.col(pu.phase_pos[idx]).sum();
}else{
R_sum[pu.phase_pos[idx]] = B_avg[pu.phase_pos[idx]] = 0.;
}
}
bool converged_MB = true; bool converged_MB = true;
bool converged_CNV = true; bool converged_CNV = true;