Reduces code duplication and reorganizes code when checking convergence.

The computations made to check the convergence are the same for all existing
phases. Therefore this patch uses loops over the phase indices when cmputing
them,

In the convergence check there are several reductions (maxCoeff(), sum())
that will trigger communication in a parallel run. This patch seperates the
reductions from the other computations. The idea is to one reduction for the
reductions that need to done as global communication is expensive.
This commit is contained in:
Markus Blatt 2015-01-19 20:22:34 +01:00
parent 7dfcf7b0c0
commit 0b81cd8f6f

View File

@ -1900,79 +1900,60 @@ namespace {
const std::vector<PhasePresence> cond = phaseCondition();
double CNVW = 0.;
double CNVO = 0.;
double CNVG = 0.;
double CNV[MaxNumPhases] = {0., 0., 0.};
double R_sum[MaxNumPhases] = {0., 0., 0.};
double B_avg[MaxNumPhases] = {0., 0., 0.};
double maxCoeff[MaxNumPhases] = {0., 0., 0.};
double mass_balance_residual[MaxNumPhases] = {0., 0., 0.};
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B;
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R;
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV;
double RW_sum = 0.;
double RO_sum = 0.;
double RG_sum = 0.;
double BW_avg = 0.;
double BO_avg = 0.;
double BG_avg = 0.;
if (active_[Water]) {
const int pos = pu.phase_pos[Water];
const ADB& tempBW = rq_[pos].b;
V BW = 1./tempBW.value();
V RW = residual_.material_balance_eq[Water].value();
BW_avg = BW.sum()/nc;
const V tempV = RW.abs()/pv;
CNVW = BW_avg * dt * tempV.maxCoeff();
RW_sum = RW.sum();
for(int idx=0; idx<MaxNumPhases; ++idx)
{
if (active_[idx]) {
const int pos = pu.phase_pos[idx];
const ADB& tempB = rq_[pos].b;
B.col(pu.phase_pos[idx]) = 1./tempB.value();
R.col(pu.phase_pos[idx]) = residual_.material_balance_eq[pu.phase_pos[idx]].value();
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)
{
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.;
}
}
if (active_[Oil]) {
// Omit the disgas here. We should add it.
const int pos = pu.phase_pos[Oil];
const ADB& tempBO = rq_[pos].b;
V BO = 1./tempBO.value();
V RO = residual_.material_balance_eq[Oil].value();
BO_avg = BO.sum()/nc;
const V tempV = RO.abs()/pv;
CNVO = BO_avg * dt * tempV.maxCoeff();
RO_sum = RO.sum();
bool converged_MB = true;
bool converged_CNV = true;
// Finish computation
for(int idx=0; idx<MaxNumPhases; ++idx)
{
CNV[idx] = B_avg[pu.phase_pos[idx]] * dt * maxCoeff[pu.phase_pos[idx]];
mass_balance_residual[idx] = std::abs(B_avg[pu.phase_pos[idx]]*R_sum[pu.phase_pos[idx]]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
}
if (active_[Gas]) {
// Omit the vapoil here. We should add it.
const int pos = pu.phase_pos[Gas];
const ADB& tempBG = rq_[pos].b;
V BG = 1./tempBG.value();
V RG = residual_.material_balance_eq[Gas].value();
BG_avg = BG.sum()/nc;
const V tempV = RG.abs()/pv;
CNVG = BG_avg * dt * tempV.maxCoeff();
RG_sum = RG.sum();
}
const double mass_balance_residual_water = std::abs(BW_avg*RW_sum) * dt / pvSum;
const double mass_balance_residual_oil = std::abs(BO_avg*RO_sum) * dt / pvSum;
const double mass_balance_residual_gas = std::abs(BG_avg*RG_sum) * dt / pvSum;
bool converged_MB = (mass_balance_residual_water < tol_mb)
&& (mass_balance_residual_oil < tol_mb)
&& (mass_balance_residual_gas < tol_mb);
bool converged_CNV = (CNVW < tol_cnv) && (CNVO < tol_cnv) && (CNVG < tol_cnv);
const double residualWellFlux = infinityNorm( residual_.well_flux_eq );
const double residualWell = infinityNorm( residual_.well_eq );
double residualWellFlux = infinityNorm(residual_.well_flux_eq);
double residualWell = infinityNorm(residual_.well_eq);
bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa);
bool converged = converged_MB && converged_CNV && converged_Well;
// if one of the residuals is NaN, throw exception, so that the solver can be restarted
if( std::isnan(mass_balance_residual_water) || mass_balance_residual_water > maxResidualAllowed() ||
std::isnan(mass_balance_residual_oil) || mass_balance_residual_oil > maxResidualAllowed() ||
std::isnan(mass_balance_residual_gas) || mass_balance_residual_gas > maxResidualAllowed() ||
std::isnan(CNVW) || CNVW > maxResidualAllowed() ||
std::isnan(CNVO) || CNVO > maxResidualAllowed() ||
std::isnan(CNVG) || CNVG > maxResidualAllowed() ||
if( std::isnan(mass_balance_residual[Water]) || mass_balance_residual[Water] > maxResidualAllowed() ||
std::isnan(mass_balance_residual[Oil]) || mass_balance_residual[Oil] > maxResidualAllowed() ||
std::isnan(mass_balance_residual[Gas]) || mass_balance_residual[Gas] > maxResidualAllowed() ||
std::isnan(CNV[Water]) || CNV[Water] > maxResidualAllowed() ||
std::isnan(CNV[Oil]) || CNV[Oil] > maxResidualAllowed() ||
std::isnan(CNV[Gas]) || CNV[Gas] > maxResidualAllowed() ||
std::isnan(residualWellFlux) || residualWellFlux > maxResidualAllowed() ||
std::isnan(residualWell) || residualWell > maxResidualAllowed() )
{
@ -1985,12 +1966,12 @@ namespace {
const std::streamsize oprec = std::cout.precision(3);
const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific);
std::cout << std::setw(4) << iteration
<< std::setw(11) << mass_balance_residual_water
<< std::setw(11) << mass_balance_residual_oil
<< std::setw(11) << mass_balance_residual_gas
<< std::setw(11) << CNVW
<< std::setw(11) << CNVO
<< std::setw(11) << CNVG
<< std::setw(11) << mass_balance_residual[Water]
<< std::setw(11) << mass_balance_residual[Oil]
<< std::setw(11) << mass_balance_residual[Gas]
<< std::setw(11) << CNV[Water]
<< std::setw(11) << CNV[Oil]
<< std::setw(11) << CNV[Gas]
<< std::setw(11) << residualWellFlux
<< std::setw(11) << residualWell
<< std::endl;