From 09f8da0324717d80c701ab647849336c9ad8bbf4 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 6 Oct 2015 13:24:34 +0200 Subject: [PATCH] fixing compilation and running after rebasing. --- opm/autodiff/BlackoilMultiSegmentModel.hpp | 2 + .../BlackoilMultiSegmentModel_impl.hpp | 136 ++++++++++-------- 2 files changed, 80 insertions(+), 58 deletions(-) diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index bb622a085..757ca2a85 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -134,6 +134,8 @@ namespace Opm { ReservoirState& reservoir_state, WellState& well_state); using Base::numPhases; + using Base::numMaterials; + using Base::materialName; protected: /* diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index d63bc0521..f9e9cd162 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -363,8 +363,15 @@ namespace Opm { const V depth = cellCentroidsZToEigen(grid_); const V pdepth = subset(depth, well_cells); std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); + // Surface density. - std::vector surf_dens(fluid_.surfaceDensity(), fluid_.surfaceDensity() + pu.num_phases); + DataBlock surf_dens(nperf, pu.num_phases); + for (int phase = 0; phase < pu.num_phases; ++ phase) { + surf_dens.col(phase) = V::Constant(nperf, fluid_.surfaceDensity()[pu.phase_pos[phase]]); + } + + std::vector surf_dens_perf(surf_dens.data(), surf_dens.data() + nperf * pu.num_phases); + // Gravity double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); @@ -372,7 +379,7 @@ namespace Opm { std::vector cd = WellDensitySegmented::computeConnectionDensities( wells(), xw, fluid_.phaseUsage(), - b_perf, rsmax_perf, rvmax_perf, surf_dens); + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // 3. Compute pressure deltas std::vector cdp = @@ -1579,102 +1586,115 @@ namespace Opm { 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 Opm::PhaseUsage& pu = fluid_.phaseUsage(); + const int np = numPhases(); + const int nm = numMaterials(); + assert(int(rq_.size()) == nm); const V pv = geo_.poreVolume(); - const std::vector cond = phaseCondition(); + std::vector R_sum(nm); + std::vector B_avg(nm); + std::vector maxCoeff(nm); + std::vector maxNormWell(np); + Eigen::Array B(nc, nm); + Eigen::Array R(nc, nm); + Eigen::Array tempV(nc, nm); - std::array CNV = {{0., 0., 0.}}; - std::array R_sum = {{0., 0., 0.}}; - std::array B_avg = {{0., 0., 0.}}; - std::array maxCoeff = {{0., 0., 0.}}; - std::array mass_balance_residual = {{0., 0., 0.}}; - std::array well_flux_residual = {{0., 0., 0.}}; - std::size_t cols = MaxNumPhases; // 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 CNV(nm); + std::vector mass_balance_residual(nm); + std::vector well_flux_residual(np); bool converged_MB = true; bool converged_CNV = true; bool converged_Well = true; // Finish computation - for ( int idx=0; idx= 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; + 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(); - // 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(CNV[Water]) || CNV[Water] > maxResidualAllowed() || - std::isnan(CNV[Oil]) || CNV[Oil] > maxResidualAllowed() || - std::isnan(CNV[Gas]) || CNV[Gas] > 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(residualWell) || residualWell > maxWellResidualAllowed ) - { - OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or too large!"); + 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 MB(WATER) MB(OIL) MB(GAS) CNVW CNVO CNVG W-FLUX(W) W-FLUX(O) W-FLUX(G)\n"; + 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 - << 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) << well_flux_residual[Water] - << std::setw(11) << well_flux_residual[Oil] - << std::setw(11) << well_flux_residual[Gas] - << std::endl; + 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); }