OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-27 18:36:28 -06:00
Code Issues Packages Projects Releases Wiki Activity

Make convergence methods work for multi-segment wells.

Browse Source 
The convergenceReduction() method no longer takes a number-of-wells argument,
instead it infers the number of well fluxes to check from the size of the
well_flux_eq member of the residual.

After this, a custom getConvergence() method is no longer required for the
multi-segment well model.
This commit is contained in:
Atgeirr Flø Rasmussen 2015-10-12 15:49:10 +02:00 committed by Kai Bao
parent 1c4f4c646d
commit e6a81fca83
4 changed files with 7 additions and 148 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
opm/autodiff/BlackoilModelBase.hpp
View File

@ -534,9 +534,8 @@ namespace Opm {
/// maximum of tempV for the phase i.
/// \param[out] B_avg An array of size MaxNumPhases where entry i contains the average
/// of B for the phase i.
/// \param[out] maxNormWell The maximum of the well equations for each phase.
/// \param[out] maxNormWell The maximum of the well flux equations for each phase.
/// \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.
double
convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
@ -546,8 +545,7 @@ namespace Opm {
std::vector<double>& maxCoeff,
std::vector<double>& B_avg,
std::vector<double>& maxNormWell,
int nc,
int nw) const;
int nc) const;
double dpMaxRel() const { return param_.dp_max_rel_; }
double dsMax() const { return param_.ds_max_; }

11
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -2513,11 +2513,12 @@ namespace detail {
std::vector<double>& maxCoeff,
std::vector<double>& B_avg,
std::vector<double>& maxNormWell,
int nc,
int nw) const
int nc) const
{
const int np = asImpl().numPhases();
const int nm = asImpl().numMaterials();
const int nw = residual_.well_flux_eq.size() / np;
assert(nw * np == int(residual_.well_flux_eq.size()));
// Do the global reductions
#if HAVE_MPI
@ -2598,7 +2599,6 @@ namespace detail {
const double tol_wells = param_.tolerance_wells_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int nw = localWellsActive() ? wells().number_of_wells : 0;
const int np = asImpl().numPhases();
const int nm = asImpl().numMaterials();
assert(int(rq_.size()) == nm);
@ -2623,7 +2623,7 @@ namespace detail {
const double pvSum = convergenceReduction(B, tempV, R,
R_sum, maxCoeff, B_avg, maxNormWell,
nc, nw);
nc);
std::vector<double> CNV(nm);
std::vector<double> mass_balance_residual(nm);
@ -2718,7 +2718,6 @@ namespace detail {
const double tol_wells = param_.tolerance_wells_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
const int nw = localWellsActive() ? wells().number_of_wells : 0;
const int np = asImpl().numPhases();
const int nm = asImpl().numMaterials();
@ -2738,7 +2737,7 @@ namespace detail {
tempV.col(idx) = R.col(idx).abs()/pv;
}
convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc, nw);
convergenceReduction(B, tempV, R, R_sum, maxCoeff, B_avg, maxNormWell, nc);
std::vector<double> well_flux_residual(np);
bool converged_Well = true;

7
opm/autodiff/BlackoilMultiSegmentModel.hpp
View File

@ -110,13 +110,6 @@ namespace Opm {
const bool initial_assembly);
/// Compute convergence based on total mass balance (tol_mb) and maximum
/// residual mass balance (tol_cnv).
/// \param[in] dt timestep length
/// \param[in] iteration current iteration number
bool getConvergence(const double dt, const int iteration);
/// Apply an update to the primary variables, chopped if appropriate.
/// \param[in] dx updates to apply to primary variables
/// \param[in, out] reservoir_state reservoir state variables

131
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
View File

@ -1682,137 +1682,6 @@ namespace Opm {
template <class Grid>
bool
BlackoilMultiSegmentModel<Grid>::getConvergence(const double dt, const int iteration)
{
const double tol_mb = param_.tolerance_mb_;
const double tol_cnv = param_.tolerance_cnv_;
const double tol_wells = param_.tolerance_wells_;
const int nc = Opm::AutoDiffGrid::numCells(grid_);
// const int nw = localWellsActive() ? wells().number_of_wells : 0;
const int nw = wellsMultiSegment().size();
// no good way to store nseg?
int nseg_total = 0;
for (int w = 0; w < nw; ++w) {
nseg_total += wellsMultiSegment()[w]->numberOfSegments();
}
const int np = numPhases();
const int nm = numMaterials();
assert(int(rq_.size()) == nm);
const V pv = geo_.poreVolume();
std::vector<double> R_sum(nm);
std::vector<double> B_avg(nm);
std::vector<double> maxCoeff(nm);
std::vector<double> maxNormWell(np);
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> B(nc, nm);
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> R(nc, nm);
Eigen::Array<V::Scalar, Eigen::Dynamic, Eigen::Dynamic> tempV(nc, nm);
for ( int idx = 0; idx < nm; ++idx )
{
const ADB& tempB = rq_[idx].b;
B.col(idx) = 1./tempB.value();
R.col(idx) = residual_.material_balance_eq[idx].value();
tempV.col(idx) = R.col(idx).abs()/pv;
}
const double pvSum = convergenceReduction(B, tempV, R,
R_sum, maxCoeff, B_avg, maxNormWell,
nc, nseg_total);
std::vector<double> CNV(nm);
std::vector<double> mass_balance_residual(nm);
std::vector<double> well_flux_residual(np);
bool converged_MB = true;
bool converged_CNV = true;
bool converged_Well = true;
// Finish computation
for ( int idx = 0; idx < nm; ++idx )
{
CNV[idx] = B_avg[idx] * dt * maxCoeff[idx];
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_CNV = converged_CNV && (CNV[idx] < tol_cnv);
// Well flux convergence is only for fluid phases, not other materials
// in our current implementation.
assert(nm >= np);
if (idx < np) {
well_flux_residual[idx] = B_avg[idx] * maxNormWell[idx];
converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells);
}
}
const double residualWell = detail::infinityNormWell(residual_.well_eq,
linsolver_.parallelInformation());
converged_Well = converged_Well && (residualWell < Opm::unit::barsa);
const bool converged = converged_MB && converged_CNV && converged_Well;
// Residual in Pascal can have high values and still be ok.
const double maxWellResidualAllowed = 1000.0 * maxResidualAllowed();
for (int idx = 0; idx < nm; ++idx) {
if (std::isnan(mass_balance_residual[idx])
|| std::isnan(CNV[idx])
|| (idx < np && std::isnan(well_flux_residual[idx]))) {
OPM_THROW(Opm::NumericalProblem, "NaN residual for phase " << materialName(idx));
}
if (mass_balance_residual[idx] > maxResidualAllowed()
|| CNV[idx] > maxResidualAllowed()
|| (idx < np && well_flux_residual[idx] > maxResidualAllowed())) {
OPM_THROW(Opm::NumericalProblem, "Too large residual for phase " << materialName(idx));
}
}
if (std::isnan(residualWell) || residualWell > maxWellResidualAllowed) {
OPM_THROW(Opm::NumericalProblem, "NaN or too large residual for well control equation");
}
if ( terminal_output_ )
{
// Only rank 0 does print to std::cout
if (iteration == 0) {
std::cout << "\nIter";
for (int idx = 0; idx < nm; ++idx) {
std::cout << " MB(" << materialName(idx).substr(0, 3) << ") ";
}
for (int idx = 0; idx < nm; ++idx) {
std::cout << " CNV(" << materialName(idx).substr(0, 1) << ") ";
}
for (int idx = 0; idx < np; ++idx) {
std::cout << " W-FLUX(" << materialName(idx).substr(0, 1) << ")";
}
std::cout << '\n';
}
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;
for (int idx = 0; idx < nm; ++idx) {
std::cout << std::setw(11) << mass_balance_residual[idx];
}
for (int idx = 0; idx < nm; ++idx) {
std::cout << std::setw(11) << CNV[idx];
}
for (int idx = 0; idx < np; ++idx) {
std::cout << std::setw(11) << well_flux_residual[idx];
}
std::cout << std::endl;
std::cout.precision(oprec);
std::cout.flags(oflags);
}
return converged;
}
template <class Grid>
void
BlackoilMultiSegmentModel<Grid>::computeSegmentDensities(const SolutionState& /* state */,