From 59bd1391a82ce73cfc05498c6f7df0a6fb1113b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 21 May 2015 16:29:01 +0200 Subject: [PATCH] Update polymer solver with convergence features from opm-autodiff. Original patches by Tor Harald Sandve and Markus Blatt. --- .../fullyimplicit/BlackoilPolymerModel.hpp | 4 +- .../BlackoilPolymerModel_impl.hpp | 159 +++++++++++++----- 2 files changed, 117 insertions(+), 46 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 3981f6760..f36ac964a 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -437,7 +437,9 @@ namespace Opm { std::array& R_sum, std::array& maxCoeff, std::array& B_avg, - int nc) const; + std::vector& maxNormWell, + int nc, + int nw) const; double dpMaxRel() const { return param_.dp_max_rel_; } double dsMax() const { return param_.ds_max_; } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 40933ed75..35d68ad77 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -1,8 +1,9 @@ /* Copyright 2013, 2015 SINTEF ICT, Applied Mathematics. - Copyright 2014 STATOIL ASA. - Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2014, 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2014, 2015 Statoil ASA. Copyright 2015 NTNU + Copyright 2015 IRIS AS This file is part of the Open Porous Media project (OPM). @@ -52,21 +53,21 @@ //#include // A debugging utility. -#define DUMP(foo) \ +#define OPM_AD_DUMP(foo) \ do { \ std::cout << "==========================================\n" \ << #foo ":\n" \ << collapseJacs(foo) << std::endl; \ } while (0) -#define DUMPVAL(foo) \ +#define OPM_AD_DUMPVAL(foo) \ do { \ std::cout << "==========================================\n" \ << #foo ":\n" \ << foo.value() << std::endl; \ } while (0) -#define DISKVAL(foo) \ +#define OPM_AD_DISKVAL(foo) \ do { \ std::ofstream os(#foo); \ os.precision(16); \ @@ -653,8 +654,8 @@ namespace detail { const int pos = pu.phase_pos[ phase ]; rq_[pos].b = fluidReciprocFVF(phase, state.canonical_phase_pressures[phase], temp, rs, rv, cond, cells_); rq_[pos].accum[aix] = pv_mult * rq_[pos].b * sat[pos]; - // DUMP(rq_[pos].b); - // DUMP(rq_[pos].accum[aix]); + // OPM_AD_DUMP(rq_[pos].b); + // OPM_AD_DUMP(rq_[pos].accum[aix]); } } @@ -669,7 +670,7 @@ namespace detail { rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix]; rq_[po].accum[aix] += state.rv * accum_gas_copy; - //DUMP(rq_[pg].accum[aix]); + // OPM_AD_DUMP(rq_[pg].accum[aix]); } if (has_polymer_) { @@ -825,14 +826,14 @@ namespace detail { computeWellConnectionPressures(state0, well_state); } - // DISKVAL(state.pressure); - // DISKVAL(state.saturation[0]); - // DISKVAL(state.saturation[1]); - // DISKVAL(state.saturation[2]); - // DISKVAL(state.rs); - // DISKVAL(state.rv); - // DISKVAL(state.qs); - // DISKVAL(state.bhp); + // OPM_AD_DISKVAL(state.pressure); + // OPM_AD_DISKVAL(state.saturation[0]); + // OPM_AD_DISKVAL(state.saturation[1]); + // OPM_AD_DISKVAL(state.saturation[2]); + // OPM_AD_DISKVAL(state.rs); + // OPM_AD_DISKVAL(state.rv); + // OPM_AD_DISKVAL(state.qs); + // OPM_AD_DISKVAL(state.bhp); // -------- Mass balance equations -------- @@ -875,7 +876,7 @@ namespace detail { residual_.material_balance_eq[ pg ] += ops_.div * (rs_face * rq_[po].mflux); residual_.material_balance_eq[ po ] += ops_.div * (rv_face * rq_[pg].mflux); - // DUMP(residual_.material_balance_eq[ Gas ]); + // OPM_AD_DUMP(residual_.material_balance_eq[ Gas ]); } @@ -1288,7 +1289,7 @@ namespace detail { } Selector alive_selector(aliveWells, Selector::NotEqualZero); residual_.well_eq = alive_selector.select(residual_.well_eq, rate_summer * state.qs); - // DUMP(residual_.well_eq); + // OPM_AD_DUMP(residual_.well_eq); } @@ -1307,15 +1308,54 @@ namespace detail { namespace detail { - - double infinityNorm( const ADB& a ) + /// \brief Compute the L-infinity norm of a vector + /// \warn This function is not suitable to compute on the well equations. + /// \param a The container to compute the infinity norm on. + /// It has to have one entry for each cell. + /// \param info In a parallel this holds the information about the data distribution. + double infinityNorm( const ADB& a, const boost::any& pinfo = boost::any() ) { +#if HAVE_MPI + if ( pinfo.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& real_info = + boost::any_cast(pinfo); + double result=0; + real_info.computeReduction(a.value(), Reduction::makeGlobalMaxFunctor(), result); + return result; + } + else +#endif + { + if( a.value().size() > 0 ) { + return a.value().matrix().lpNorm (); + } + else { // this situation can occur when no wells are present + return 0.0; + } + } + (void)pinfo; // Suppress unused warning for non-MPI. + } + + /// \brief Compute the L-infinity norm of a vector representing a well equation. + /// \param a The container to compute the infinity norm on. + /// \param info In a parallel this holds the information about the data distribution. + double infinityNormWell( const ADB& a, const boost::any& pinfo ) + { + double result=0; if( a.value().size() > 0 ) { - return a.value().matrix().lpNorm (); + result = a.value().matrix().lpNorm (); } - else { // this situation can occur when no wells are present - return 0.0; +#if HAVE_MPI + if ( pinfo.type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& real_info = + boost::any_cast(pinfo); + result = real_info.communicator().max(result); } +#endif + return result; + (void)pinfo; // Suppress unused warning for non-MPI. } } // namespace detail @@ -1764,8 +1804,8 @@ namespace detail { const ADB& b = rq_[ actph ].b; const ADB& mob = rq_[ actph ].mob; rq_[ actph ].mflux = upwind.select(b * mob) * head; - // DUMP(rq_[ actph ].mob); - // DUMP(rq_[ actph ].mflux); + // OPM_AD_DUMP(rq_[ actph ].mob); + // OPM_AD_DUMP(rq_[ actph ].mflux); } @@ -1844,12 +1884,14 @@ namespace detail { template double BlackoilPolymerModel::convergenceReduction(const Eigen::Array& B, - const Eigen::Array& tempV, - const Eigen::Array& R, - std::array& R_sum, - std::array& maxCoeff, - std::array& B_avg, - int nc) const + const Eigen::Array& tempV, + const Eigen::Array& R, + std::array& R_sum, + std::array& maxCoeff, + std::array& B_avg, + std::vector& maxNormWell, + int nc, + int nw) const { // Do the global reductions #if HAVE_MPI @@ -1877,16 +1919,23 @@ namespace detail { Opm::Reduction::makeGlobalMaxFunctor(), Opm::Reduction::makeGlobalSumFunctor()); info.computeReduction(containers, operators, values); - B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv); - maxCoeff[idx] = std::get<1>(values); - R_sum[idx] = std::get<2>(values); + B_avg[idx] = std::get<0>(values)/std::get<0>(nc_and_pv); + maxCoeff[idx] = std::get<1>(values); + R_sum[idx] = std::get<2>(values); + maxNormWell[idx] = 0.0; + for ( int w=0; w(nc_and_pv); } else @@ -1903,6 +1952,11 @@ namespace detail { { R_sum[idx] = B_avg[idx] = maxCoeff[idx] =0.0; } + maxNormWell[idx] = 0.0; + for ( int w=0; w bool BlackoilPolymerModel::getConvergence(const double dt, const int iteration) @@ -1923,6 +1980,7 @@ namespace detail { const double tol_wells = param_.tolerance_wells_; const int nc = Opm::AutoDiffGrid::numCells(grid_); + const int nw = wellsActive() ? wells().number_of_wells : 0; const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const V pv = geo_.poreVolume(); @@ -1934,10 +1992,12 @@ namespace detail { std::array B_avg = {{0., 0., 0., 0.}}; std::array maxCoeff = {{0., 0., 0., 0.}}; std::array mass_balance_residual = {{0., 0., 0., 0.}}; + std::array well_flux_residual = {{0., 0., 0., 0.}}; std::size_t cols = MaxNumPhases+1; // needed to pass the correct type to Eigen Eigen::Array B(nc, cols); Eigen::Array R(nc, cols); Eigen::Array tempV(nc, cols); + std::vector maxNormWell(MaxNumPhases); for ( int idx=0; idx 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(mass_balance_residual[MaxNumPhases]) || mass_balance_residual[MaxNumPhases] > maxResidualAllowed() || - std::isnan(CNV[Water]) || CNV[Water] > maxResidualAllowed() || + std::isnan(mass_balance_residual[MaxNumPhases]) || mass_balance_residual[MaxNumPhases] > maxResidualAllowed() || std::isnan(CNV[Water]) || CNV[Water] > maxResidualAllowed() || std::isnan(CNV[Oil]) || CNV[Oil] > maxResidualAllowed() || std::isnan(CNV[Gas]) || CNV[Gas] > maxResidualAllowed() || std::isnan(CNV[MaxNumPhases]) || CNV[MaxNumPhases] > maxResidualAllowed() || - std::isnan(residualWellFlux) || residualWellFlux > maxResidualAllowed() || + std::isnan(well_flux_residual[Water]) || well_flux_residual[Water] > maxResidualAllowed() || + std::isnan(well_flux_residual[Oil]) || well_flux_residual[Oil] > maxResidualAllowed() || + std::isnan(well_flux_residual[Gas]) || well_flux_residual[Gas] > maxResidualAllowed() || + std::isnan(well_flux_residual[MaxNumPhases]) || well_flux_residual[MaxNumPhases] > maxResidualAllowed() || std::isnan(residualWell) || residualWell > maxResidualAllowed() ) { OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or to large!"); @@ -1993,7 +2060,7 @@ namespace detail { { // Only rank 0 does print to std::cout if (iteration == 0) { - std::cout << "\nIter MB(OIL) MB(WATER) MB(GAS) CNVW CNVO CNVG CNVP WELL-FLOW WELL-CNTRL\n"; + std::cout << "\nIter MB(WATER) MB(OIL) MB(GAS) MB(POLY) CNVW CNVO CNVG CNVP W-FLUX(W) W-FLUX(O) W-FLUX(G) W-FLUX(P)\n"; } const std::streamsize oprec = std::cout.precision(3); const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific); @@ -2006,8 +2073,10 @@ namespace detail { << std::setw(11) << CNV[Oil] << std::setw(11) << CNV[Gas] << std::setw(11) << CNV[MaxNumPhases] - << std::setw(11) << residualWellFlux - << std::setw(11) << residualWell + << std::setw(11) << well_flux_residual[Water] + << std::setw(11) << well_flux_residual[Oil] + << std::setw(11) << well_flux_residual[Gas] + << std::setw(11) << well_flux_residual[MaxNumPhases] << std::endl; std::cout.precision(oprec); std::cout.flags(oflags);