Compute scaled well residuals correctly for parallel runs.

This commit adapts the PR #375 for parallel runs. That is, the norms are
calculated over all wells, not just the ones that perforate the local grid
cells.
As this is a reduction, too, we move the computation to convergenceReduction
method.
This commit is contained in:
Markus Blatt 2015-05-13 12:00:38 +02:00
parent 9e0b2fed4f
commit 48ce90fcc7
2 changed files with 27 additions and 14 deletions

View File

@ -389,7 +389,9 @@ namespace Opm {
/// maximum of tempV for the phase i. /// maximum of tempV for the phase i.
/// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average /// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
/// of B for the phase i. /// of B for the phase i.
/// \param[out] maxNormWell The maximum of the well equations for each phase.
/// \param[in] nc The number of cells of the local grid. /// \param[in] nc The number of cells of the local grid.
/// \param[in] nw The number of wells on the local grid.
/// \return The total pore volume over all cells. /// \return The total pore volume over all cells.
double double
convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B, convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, MaxNumPhases>& B,
@ -398,7 +400,9 @@ namespace Opm {
std::array<double,MaxNumPhases>& R_sum, std::array<double,MaxNumPhases>& R_sum,
std::array<double,MaxNumPhases>& maxCoeff, std::array<double,MaxNumPhases>& maxCoeff,
std::array<double,MaxNumPhases>& B_avg, std::array<double,MaxNumPhases>& B_avg,
int nc) const; std::vector<double>& maxNormWell,
int nc,
int nw) 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,

View File

@ -1834,7 +1834,9 @@ namespace detail {
std::array<double,MaxNumPhases>& R_sum, std::array<double,MaxNumPhases>& R_sum,
std::array<double,MaxNumPhases>& maxCoeff, std::array<double,MaxNumPhases>& maxCoeff,
std::array<double,MaxNumPhases>& B_avg, std::array<double,MaxNumPhases>& B_avg,
int nc) const std::vector<double>& maxNormWell,
int nc,
int nw) const
{ {
// Do the global reductions // Do the global reductions
#if HAVE_MPI #if HAVE_MPI
@ -1865,12 +1867,18 @@ namespace detail {
B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv); B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv);
maxCoeff[idx] = std::get<1>(values); maxCoeff[idx] = std::get<1>(values);
R_sum[idx] = std::get<2>(values); R_sum[idx] = std::get<2>(values);
maxNormWell[idx] = 0.0;
for ( int w=0; w<nw; ++w )
{
maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
}
} }
else else
{ {
R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0; maxNormWell[idx] = R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0;
} }
} }
info.communicator().max(&maxNormWell[0], MaxNumPhases);
// Compute pore volume // Compute pore volume
return std::get<1>(nc_and_pv); return std::get<1>(nc_and_pv);
} }
@ -1888,6 +1896,11 @@ namespace detail {
{ {
R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0; R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0;
} }
maxNormWell[idx] = 0.0;
for ( int w=0; w<nw; ++w )
{
maxNormWell[idx] = std::max(maxNormWell[idx], std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
}
} }
// Compute total pore volume // Compute total pore volume
return geo_.poreVolume().sum(); return geo_.poreVolume().sum();
@ -1920,6 +1933,7 @@ namespace detail {
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B(nc, cols); Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> B(nc, cols);
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R(nc, cols); Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> R(nc, cols);
Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV(nc, cols); Eigen::Array<V::Scalar, Eigen::Dynamic, MaxNumPhases> tempV(nc, cols);
std::vector<double> maxNormWell(MaxNumPhases);
for ( int idx=0; idx<MaxNumPhases; ++idx ) for ( int idx=0; idx<MaxNumPhases; ++idx )
{ {
@ -1932,7 +1946,8 @@ namespace detail {
} }
} }
const double pvSum = convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, nc); const double pvSum = convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg,
maxNormWell, nc, nw);
bool converged_MB = true; bool converged_MB = true;
bool converged_CNV = true; bool converged_CNV = true;
@ -1944,13 +1959,7 @@ namespace detail {
mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum; mass_balance_residual[idx] = std::abs(B_avg[idx]*R_sum[idx]) * dt / pvSum;
converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb); converged_MB = converged_MB && (mass_balance_residual[idx] < tol_mb);
converged_CNV = converged_CNV && (CNV[idx] < tol_cnv); converged_CNV = converged_CNV && (CNV[idx] < tol_cnv);
well_flux_residual[idx] = B_avg[idx] * dt * maxNormWell[idx];
double maxNormWell = 0.0;
for ( int w=0; w<nw; ++w )
{
maxNormWell = std::max(maxNormWell, std::abs(residual_.well_flux_eq.value()[nw*idx + w]));
}
well_flux_residual[idx] = B_avg[idx] * dt * maxNormWell;
converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
} }