From 5f6027ba0183177b5417c43827e16abd9d19dc18 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 25 May 2015 00:06:17 +0200 Subject: [PATCH] Bugfix: we compute no well flux residual for polymer, do not try to use. --- .../BlackoilPolymerModel_impl.hpp | 38 ++++++++++--------- 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 35d68ad77..d9be4b553 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -1908,7 +1908,8 @@ namespace detail { auto nc_and_pv_containers = std::make_tuple(v, geo_.poreVolume()); info.computeReduction(nc_and_pv_containers, nc_and_pv_operators, nc_and_pv); - for ( int idx=0; idx(0.0 ,0.0 ,0.0); @@ -1933,7 +1934,7 @@ namespace detail { maxNormWell[idx] = R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0; } } - info.communicator().max(&maxNormWell[0], MaxNumPhases); + info.communicator().max(&maxNormWell[0], MaxNumPhases+1); // Compute pore volume #warning Missing polymer code for MPI version return std::get<1>(nc_and_pv); @@ -1941,9 +1942,9 @@ namespace detail { else #endif { - for ( int idx=0; idx 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::array well_flux_residual = {{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); @@ -2029,9 +2031,10 @@ namespace detail { 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_residual[idx] = B_avg[idx] * dt * maxNormWell[idx]; - - converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); + if (idx != MaxNumPhases) { // No well flux residual for polymer. + well_flux_residual[idx] = B_avg[idx] * dt * maxNormWell[idx]; + converged_Well = converged_Well && (well_flux_residual[idx] < tol_wells); + } } const double residualWell = detail::infinityNormWell(residual_.well_eq, @@ -2040,27 +2043,27 @@ namespace detail { const 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() || + 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(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(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!"); + OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!"); } if ( terminal_output_ ) { // Only rank 0 does print to std::cout if (iteration == 0) { - 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"; + std::cout << "\nIter MB(WATER) MB(OIL) MB(GAS) MB(POLY) CNVW CNVO CNVG CNVP W-FLUX(W) W-FLUX(O) W-FLUX(G)\n"; } const std::streamsize oprec = std::cout.precision(3); const std::ios::fmtflags oflags = std::cout.setf(std::ios::scientific); @@ -2076,7 +2079,6 @@ namespace detail { << 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);