diff --git a/examples/sim_poly_fibo_ad.cpp b/examples/sim_poly_fibo_ad.cpp deleted file mode 100644 index 9189fe0db..000000000 --- a/examples/sim_poly_fibo_ad.cpp +++ /dev/null @@ -1,252 +0,0 @@ -/* - Copyright 2014 SINTEF ICT, Applied Mathematics. - Copyright 2014 STATOIL ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ -#include "config.h" - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include - -#include -#include - -#include -#include -#include -#include -#include -#include - -#include -#include - -#include -#include - -#include -#include -#include -#include -#include - - -namespace -{ - void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param) - { - if (param.anyUnused()) { - std::cout << "-------------------- Unused parameters: --------------------\n"; - param.displayUsage(); - std::cout << "----------------------------------------------------------------" << std::endl; - } - } -} // anon namespace - - - -// ----------------- Main program ----------------- -int -main(int argc, char** argv) -try -{ - using namespace Opm; - - std::cout << "\n================ Test program for fully implicit three-phase black-oil flow ===============\n\n"; - parameter::ParameterGroup param(argc, argv, false); - std::cout << "--------------- Reading parameters ---------------" << std::endl; - - // If we have a "deck_filename", grid and props will be read from that. - bool use_deck = param.has("deck_filename"); - if (!use_deck) { - OPM_THROW(std::runtime_error, "This program must be run with an input deck. " - "Specify the deck with deck_filename=deckname.data (for example)."); - } - std::shared_ptr grid; - std::shared_ptr props; - std::shared_ptr new_props; - std::shared_ptr rock_comp; - PolymerBlackoilState state; - // bool check_well_controls = false; - // int max_well_control_iterations = 0; - double gravity[3] = { 0.0 }; - std::string deck_filename = param.get("deck_filename"); - - Opm::ParserPtr newParser(new Opm::Parser() ); - Opm::DeckConstPtr deck = newParser->parseFile(deck_filename); - std::shared_ptr eclipseState(new EclipseState(deck)); - - // Grid init - std::vector porv; - if (eclipseState->hasDoubleGridProperty("PORV")) { - porv = eclipseState->getDoubleGridProperty("PORV")->getData(); - } - grid.reset(new GridManager(eclipseState->getEclipseGrid(), porv)); - auto &cGrid = *grid->c_grid(); - const PhaseUsage pu = Opm::phaseUsageFromDeck(deck); - Opm::EclipseWriter outputWriter(param, - eclipseState, - pu, - cGrid.number_of_cells, - cGrid.global_cell); - - // Rock and fluid init - props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid(), param)); - new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, *grid->c_grid())); - const bool polymer = deck->hasKeyword("POLYMER"); - const bool use_wpolymer = deck->hasKeyword("WPOLYMER"); - PolymerProperties polymer_props(deck, eclipseState); - PolymerPropsAd polymer_props_ad(polymer_props); - // check_well_controls = param.getDefault("check_well_controls", false); - // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); - // Rock compressibility. - rock_comp.reset(new RockCompressibility(deck, eclipseState)); - - // Gravity. - gravity[2] = deck->hasKeyword("NOGRAV") ? 0.0 : unit::gravity; - - // Init state variables (saturation and pressure). - if (param.has("init_saturation")) { - initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); - initBlackoilSurfvol(*grid->c_grid(), *props, state); - enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; - if (pu.phase_used[Oil] && pu.phase_used[Gas]) { - const int np = props->numPhases(); - const int nc = grid->c_grid()->number_of_cells; - for (int c = 0; c < nc; ++c) { - state.gasoilratio()[c] = state.surfacevol()[c*np + pu.phase_pos[Gas]] - / state.surfacevol()[c*np + pu.phase_pos[Oil]]; - } - } - } else if (deck->hasKeyword("EQUIL") && props->numPhases() == 3) { - state.init(*grid->c_grid(), props->numPhases()); - const double grav = param.getDefault("gravity", unit::gravity); - initStateEquil(*grid->c_grid(), *props, deck, eclipseState, grav, state.blackoilState()); - state.faceflux().resize(grid->c_grid()->number_of_faces, 0.0); - } else { - state.init(*grid->c_grid(), props->numPhases()); - initBlackoilStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state.blackoilState()); - } - - bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); - const double *grav = use_gravity ? &gravity[0] : 0; - - // Solver for Newton iterations. - std::unique_ptr fis_solver; - if (param.getDefault("use_cpr", true)) { - fis_solver.reset(new NewtonIterationBlackoilCPR(param)); - } else { - fis_solver.reset(new NewtonIterationBlackoilSimple(param)); - } - - // Write parameters used for later reference. - bool output = param.getDefault("output", true); - std::string output_dir; - if (output) { - // Create output directory if needed. - output_dir = - param.getDefault("output_dir", std::string("output")); - boost::filesystem::path fpath(output_dir); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - // Write simulation parameters. - param.writeParam(output_dir + "/simulation.param"); - } - - Opm::TimeMapConstPtr timeMap(eclipseState->getSchedule()->getTimeMap()); - SimulatorTimer simtimer; - - // initialize variables - simtimer.init(timeMap); - if (polymer){ - if (!use_wpolymer) { - OPM_MESSAGE("Warning: simulate polymer injection without WPOLYMER."); - } else { - if (param.has("polymer_start_days")) { - OPM_MESSAGE("Warning: Using WPOLYMER to control injection since it was found in deck." - "You seem to be trying to control it via parameter poly_start_days (etc.) as well."); - } - } - } else { - if (use_wpolymer) { - OPM_MESSAGE("Warning: use WPOLYMER in a non-polymer scenario."); - } - } - std::cout << "\n\n================ Starting main simulation loop ===============\n" - << std::flush; - SimulatorReport fullReport; - Opm::DerivedGeology geology(*grid->c_grid(), *new_props, eclipseState, grav); - - std::vector threshold_pressures = thresholdPressures(eclipseState, *grid->c_grid()); - SimulatorFullyImplicitBlackoilPolymer simulator(param, - *grid->c_grid(), - geology, - *new_props, - polymer_props_ad, - rock_comp->isActive() ? rock_comp.get() : 0, - *fis_solver, - grav, - deck->hasKeyword("DISGAS"), - deck->hasKeyword("VAPOIL"), - polymer, - eclipseState, - outputWriter, - deck, - threshold_pressures); - - - fullReport = simulator.run(simtimer, state); - std::cout << "\n\n================ End of simulation ===============\n\n"; - fullReport.report(std::cout); - - if (output) { - std::string filename = output_dir + "/walltime.txt"; - std::fstream tot_os(filename.c_str(),std::fstream::trunc | std::fstream::out); - fullReport.reportParam(tot_os); - warnIfUnusedParams(param); - } -} -catch (const std::exception &e) { - std::cerr << "Program threw an exception: " << e.what() << "\n"; - throw; -} - diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp index 47f6742c3..541ab69c2 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver.hpp @@ -1,7 +1,6 @@ /* Copyright 2013 SINTEF ICT, Applied Mathematics. Copyright 2014 STATOIL ASA. - This file is part of the Open Porous Media project (OPM). @@ -22,6 +21,8 @@ #ifndef OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED #define OPM_FULLYIMPLICITBLACKOILPOLYMERSOLVER_HEADER_INCLUDED +#include + #include #include #include @@ -30,6 +31,8 @@ #include #include +#include + struct UnstructuredGrid; struct Wells; @@ -56,29 +59,55 @@ namespace Opm { class FullyImplicitBlackoilPolymerSolver { public: + // the Newton relaxation type + enum RelaxType { DAMPEN, SOR }; + + // class holding the solver parameters + struct SolverParameter + { + double dp_max_rel_; + double ds_max_; + double dr_max_rel_; + enum RelaxType relax_type_; + double relax_max_; + double relax_increment_; + double relax_rel_tol_; + double max_residual_allowed_; + double tolerance_mb_; + double tolerance_cnv_; + double tolerance_wells_; + int max_iter_; + + SolverParameter( const parameter::ParameterGroup& param ); + SolverParameter(); + + void reset(); + }; + /// \brief The type of the grid that we use. typedef T Grid; /// Construct a solver. It will retain references to the /// arguments of this functions, and they are expected to /// remain in scope for the lifetime of the solver. - /// \param[in] param parameters + /// \param[in] param parameters /// \param[in] grid grid data structure /// \param[in] fluid fluid properties /// \param[in] geo rock properties /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] wells well structure /// \param[in] linsolver linear solver - FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param, + FullyImplicitBlackoilPolymerSolver(const SolverParameter& param, const Grid& grid , const BlackoilPropsAdInterface& fluid, const DerivedGeology& geo , const RockCompressibility* rock_comp_props, const PolymerPropsAd& polymer_props_ad, - const Wells& wells, + const Wells* wells, const NewtonIterationBlackoilInterface& linsolver, const bool has_disgas, const bool has_vapoil, - const bool has_polymer); + const bool has_polymer, + const bool terminal_output); /// \brief Set threshold pressures that prevent or reduce flow. /// This prevents flow across faces if the potential @@ -99,12 +128,16 @@ namespace Opm { /// \param[in] dt time step size /// \param[in] state reservoir state /// \param[in] wstate well state - void + /// \return number of linear iterations used + int step(const double dt , PolymerBlackoilState& state , WellStateFullyImplicitBlackoil& wstate, const std::vector& polymer_inflow); + unsigned int newtonIterations () const { return newtonIterations_; } + unsigned int linearIterations () const { return linearIterations_; } + private: // Types and enums typedef AutoDiffBlock ADB; @@ -133,21 +166,23 @@ namespace Opm { ADB rv; ADB concentration; ADB qs; - ADB bhp; + ADB bhp; + // Below are quantities stored in the state for optimization purposes. + std::vector canonical_phase_pressures; // Always has 3 elements, even if only 2 phases active. }; struct WellOps { - WellOps(const Wells& wells); + WellOps(const Wells* wells); M w2p; // well -> perf (scatter) M p2w; // perf -> well (gather) }; - enum { Water = BlackoilPropsAdInterface::Water, - Oil = BlackoilPropsAdInterface::Oil , - Gas = BlackoilPropsAdInterface::Gas }; + enum { Water = BlackoilPropsAdInterface::Water, + Oil = BlackoilPropsAdInterface::Oil , + Gas = BlackoilPropsAdInterface::Gas , + MaxNumPhases = BlackoilPropsAdInterface::MaxNumPhases + }; - // the Newton relaxation type - enum RelaxType { DAMPEN, SOR }; enum PrimalVariables { Sg = 0, RS = 1, RV = 2 }; // Member data @@ -155,8 +190,8 @@ namespace Opm { const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_; const RockCompressibility* rock_comp_props_; - const PolymerPropsAd& polymer_props_ad_; - const Wells& wells_; + const PolymerPropsAd& polymer_props_ad_; + const Wells* wells_; const NewtonIterationBlackoilInterface& linsolver_; // For each canonical phase -> true if active const std::vector active_; @@ -170,14 +205,8 @@ namespace Opm { const bool has_vapoil_; const bool has_polymer_; const int poly_pos_; - double dp_max_rel_; - double ds_max_; - double drs_max_rel_; - enum RelaxType relax_type_; - double relax_max_; - double relax_increment_; - double relax_rel_tol_; - int max_iter_; + + SolverParameter param_; bool use_threshold_pressure_; V threshold_pressures_by_interior_face_; @@ -187,16 +216,30 @@ namespace Opm { LinearisedBlackoilResidual residual_; + /// \brief Whether we print something to std::cout + bool terminal_output_; + unsigned int newtonIterations_; + unsigned int linearIterations_; + std::vector primalVariable_; // Private methods. + + // return true if wells are available + bool wellsActive() const { return wells_ ? wells_->number_of_wells > 0 : false ; } + // return wells object + const Wells& wells () const { assert( bool(wells_ != 0) ); return *wells_; } + SolutionState constantState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw); + const WellStateFullyImplicitBlackoil& xw) const; + + void + makeConstantState(SolutionState& state) const; SolutionState variableState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw); + const WellStateFullyImplicitBlackoil& xw) const; void computeAccum(const SolutionState& state, @@ -223,6 +266,7 @@ namespace Opm { void assemble(const V& dtpv, const PolymerBlackoilState& x, + const bool initial_assembly, WellStateFullyImplicitBlackoil& xw, const std::vector& polymer_inflow); @@ -235,6 +279,18 @@ namespace Opm { std::vector computePressures(const SolutionState& state) const; + std::vector + computePressures(const ADB& po, + const ADB& sw, + const ADB& so, + const ADB& sg) const; + + V + computeGasPressure(const V& po, + const V& sw, + const V& so, + const V& sg) const; + std::vector computeRelPerm(const SolutionState& state) const; @@ -244,35 +300,26 @@ namespace Opm { const std::vector& well_cells) const; void - computeMassFlux(const int actph , - const V& transi, - const ADB& kr , - const ADB& p , - const SolutionState& state ); - - void - computeMassFlux(const V& trans, - const std::vector& kr, + computeMassFlux(const V& transi, + const std::vector& kr , const std::vector& phasePressure, - const SolutionState& state); - + const SolutionState& state ); void computeCmax(PolymerBlackoilState& state); ADB computeMc(const SolutionState& state) const; - ADB - rockPorosity(const ADB& p) const; - - ADB - rockPermeability(const ADB& p) const; void applyThresholdPressures(ADB& dp); double residualNorm() const; - std::vector residuals() const; + /// \brief Compute the residual norms of the mass balance for each phase, + /// the well flux, and the well equation. + /// \return a vector that contains for each phase the norm of the mass balance + /// and afterwards the norm of the residual of the well flux and the well equation. + std::vector computeResidualNorms() const; ADB fluidViscosity(const int phase, @@ -321,7 +368,6 @@ namespace Opm { const ADB& so, const std::vector& cells) const; - ADB poroMult(const ADB& p) const; @@ -350,7 +396,35 @@ namespace Opm { /// Compute convergence based on total mass balance (tol_mb) and maximum /// residual mass balance (tol_cnv). - bool getConvergence(const double dt, const int it); + bool getConvergence(const double dt, const int iteration); + + /// \brief Compute the reduction within the convergence check. + /// \param[in] B A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. B.col(i) contains the values + /// for phase i. + /// \param[in] tempV A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. tempV.col(i) contains the + /// values + /// for phase i. + /// \param[in] R A matrix with MaxNumPhases columns and the same number rows + /// as the number of cells of the grid. B.col(i) contains the values + /// for phase i. + /// \param[out] R_sum An array of size MaxNumPhases where entry i contains the sum + /// of R for the phase i. + /// \param[out] maxCoeff An array of size MaxNumPhases where entry i contains the + /// 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[in] nc The number of cells of the local grid. + /// \return The total pore volume over all cells. + double + 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; void detectNewtonOscillations(const std::vector>& residual_history, const int it, const double relaxRelTol, @@ -358,14 +432,15 @@ namespace Opm { void stablizeNewton(V& dx, V& dxOld, const double omega, const RelaxType relax_type) const; - double dpMaxRel() const { return dp_max_rel_; } - double dsMax() const { return ds_max_; } - double drsMaxRel() const { return drs_max_rel_; } - enum RelaxType relaxType() const { return relax_type_; } - double relaxMax() const { return relax_max_; }; - double relaxIncrement() const { return relax_increment_; }; - double relaxRelTol() const { return relax_rel_tol_; }; - double maxIter() const { return max_iter_; } + double dpMaxRel() const { return param_.dp_max_rel_; } + double dsMax() const { return param_.ds_max_; } + double drMaxRel() const { return param_.dr_max_rel_; } + enum RelaxType relaxType() const { return param_.relax_type_; } + double relaxMax() const { return param_.relax_max_; }; + double relaxIncrement() const { return param_.relax_increment_; }; + double relaxRelTol() const { return param_.relax_rel_tol_; }; + double maxIter() const { return param_.max_iter_; } + double maxResidualAllowed() const { return param_.max_residual_allowed_; } }; } // namespace Opm diff --git a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp index 28826f5b5..af74ffe21 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitBlackoilPolymerSolver_impl.hpp @@ -1,6 +1,8 @@ /* Copyright 2013 SINTEF ICT, Applied Mathematics. Copyright 2014 STATOIL ASA. + Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services + Copyright 2015 NTNU This file is part of the Open Porous Media project (OPM). @@ -31,6 +33,7 @@ #include #include +#include #include #include #include @@ -42,6 +45,7 @@ #include #include #include +#include //#include // A debugging utility. @@ -78,7 +82,7 @@ typedef Eigen::Array DataBlock; -namespace { +namespace detail { std::vector @@ -93,35 +97,6 @@ namespace { - template - V computePerfPress(const Grid& grid, const Wells& wells, const V& rho, const double grav) - { - using namespace Opm::AutoDiffGrid; - const int nw = wells.number_of_wells; - const int nperf = wells.well_connpos[nw]; - const int dim = dimensions(grid); - V wdp = V::Zero(nperf,1); - assert(wdp.size() == rho.size()); - - // Main loop, iterate over all perforations, - // using the following formula: - // wdp(perf) = g*(perf_z - well_ref_z)*rho(perf) - // where the total density rho(perf) is taken to be - // sum_p (rho_p*saturation_p) in the perforation cell. - // [although this is computed on the outside of this function]. - for (int w = 0; w < nw; ++w) { - const double ref_depth = wells.depth_ref[w]; - for (int j = wells.well_connpos[w]; j < wells.well_connpos[w + 1]; ++j) { - const int cell = wells.well_cells[j]; - const double cell_depth = cellCentroid(grid, cell)[dim - 1]; - wdp[j] = rho[j]*grav*(cell_depth - ref_depth); - } - } - return wdp; - } - - - template std::vector activePhases(const PU& pu) @@ -155,6 +130,7 @@ namespace { } + template int polymerPos(const PU& pu) { @@ -168,69 +144,54 @@ namespace { return pos; } -} // Anonymous namespace - +} // namespace detail template - FullyImplicitBlackoilPolymerSolver:: - FullyImplicitBlackoilPolymerSolver(const parameter::ParameterGroup& param, - const Grid& grid , - const BlackoilPropsAdInterface& fluid, - const DerivedGeology& geo , - const RockCompressibility* rock_comp_props, - const PolymerPropsAd& polymer_props_ad, - const Wells& wells, - const NewtonIterationBlackoilInterface& linsolver, - const bool has_disgas, - const bool has_vapoil, - const bool has_polymer) - : grid_ (grid) - , fluid_ (fluid) - , geo_ (geo) - , rock_comp_props_(rock_comp_props) - , polymer_props_ad_(polymer_props_ad) - , wells_ (wells) - , linsolver_ (linsolver) - , active_(activePhases(fluid.phaseUsage())) - , canph_ (active2Canonical(fluid.phaseUsage())) - , cells_ (buildAllCells(Opm::AutoDiffGrid::numCells(grid))) - , ops_ (grid) - , wops_ (wells) - , cmax_(V::Zero(Opm::AutoDiffGrid::numCells(grid))) - , has_disgas_(has_disgas) - , has_vapoil_(has_vapoil) - , has_polymer_(has_polymer) - , poly_pos_(polymerPos(fluid.phaseUsage())) - , dp_max_rel_ (1.0e9) - , ds_max_ (0.2) - , drs_max_rel_ (1.0e9) - , relax_type_ (DAMPEN) - , relax_max_ (0.5) - , relax_increment_ (0.1) - , relax_rel_tol_ (0.2) - , max_iter_ (15) - , use_threshold_pressure_(false) - , rq_ (fluid.numPhases()) - , phaseCondition_(AutoDiffGrid::numCells(grid)) - , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), - ADB::null(), - ADB::null() } ) + void FullyImplicitBlackoilPolymerSolver::SolverParameter:: + reset() { - if (has_polymer_) { - if (!active_[Water]) { - OPM_THROW(std::logic_error, "Polymer must solved in water!\n"); - } - // If deck has polymer, residual_ should contain polymer equation. - rq_.resize(fluid_.numPhases()+1); - residual_.material_balance_eq.resize(fluid_.numPhases()+1, ADB::null()); - assert(poly_pos_ == fluid_.numPhases()); - } - dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_); - ds_max_ = param.getDefault("ds_max", ds_max_); - drs_max_rel_ = param.getDefault("drs_max_rel", drs_max_rel_); - relax_max_ = param.getDefault("relax_max", relax_max_); - max_iter_ = param.getDefault("max_iter", max_iter_); + // default values for the solver parameters + dp_max_rel_ = 1.0e9; + ds_max_ = 0.2; + dr_max_rel_ = 1.0e9; + relax_type_ = DAMPEN; + relax_max_ = 0.5; + relax_increment_ = 0.1; + relax_rel_tol_ = 0.2; + max_iter_ = 15; + max_residual_allowed_ = std::numeric_limits< double >::max(); + tolerance_mb_ = 1.0e-7; + tolerance_cnv_ = 1.0e-3; + tolerance_wells_ = 1./Opm::unit::day; + } + + template + FullyImplicitBlackoilPolymerSolver::SolverParameter:: + SolverParameter() + { + // set default values + reset(); + } + + template + FullyImplicitBlackoilPolymerSolver::SolverParameter:: + SolverParameter( const parameter::ParameterGroup& param ) + { + // set default values + reset(); + + // overload with given parameters + dp_max_rel_ = param.getDefault("dp_max_rel", dp_max_rel_); + ds_max_ = param.getDefault("ds_max", ds_max_); + dr_max_rel_ = param.getDefault("dr_max_rel", dr_max_rel_); + relax_max_ = param.getDefault("relax_max", relax_max_); + max_iter_ = param.getDefault("max_iter", max_iter_); + max_residual_allowed_ = param.getDefault("max_residual_allowed", max_residual_allowed_); + + tolerance_mb_ = param.getDefault("tolerance_mb", tolerance_mb_); + tolerance_cnv_ = param.getDefault("tolerance_cnv", tolerance_cnv_); + tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ ); std::string relaxation_type = param.getDefault("relax_type", std::string("dampen")); if (relaxation_type == "dampen") { @@ -243,6 +204,71 @@ namespace { } + template + FullyImplicitBlackoilPolymerSolver:: + FullyImplicitBlackoilPolymerSolver(const SolverParameter& param, + const Grid& grid , + const BlackoilPropsAdInterface& fluid, + const DerivedGeology& geo , + const RockCompressibility* rock_comp_props, + const PolymerPropsAd& polymer_props_ad, + const Wells* wells, + const NewtonIterationBlackoilInterface& linsolver, + const bool has_disgas, + const bool has_vapoil, + const bool has_polymer, + const bool terminal_output) + : grid_ (grid) + , fluid_ (fluid) + , geo_ (geo) + , rock_comp_props_(rock_comp_props) + , polymer_props_ad_(polymer_props_ad) + , wells_ (wells) + , linsolver_ (linsolver) + , active_(detail::activePhases(fluid.phaseUsage())) + , canph_ (detail::active2Canonical(fluid.phaseUsage())) + , cells_ (detail::buildAllCells(Opm::AutoDiffGrid::numCells(grid))) + , ops_ (grid) + , wops_ (wells_) + , cmax_(V::Zero(Opm::AutoDiffGrid::numCells(grid))) + , has_disgas_(has_disgas) + , has_vapoil_(has_vapoil) + , has_polymer_(has_polymer) + , poly_pos_(detail::polymerPos(fluid.phaseUsage())) + , param_( param ) + , use_threshold_pressure_(false) + , rq_ (fluid.numPhases()) + , phaseCondition_(AutoDiffGrid::numCells(grid)) + , residual_ ( { std::vector(fluid.numPhases(), ADB::null()), + ADB::null(), + ADB::null() } ) + , terminal_output_ (terminal_output) + , newtonIterations_( 0 ) + , linearIterations_( 0 ) + { +#if HAVE_MPI + if ( terminal_output_ ) { + if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast(linsolver_.parallelInformation()); + // Only rank 0 does print to std::cout if terminal_output is enabled + terminal_output_ = (info.communicator().rank()==0); + } + } +#endif + if (has_polymer_) { + if (!active_[Water]) { + OPM_THROW(std::logic_error, "Polymer must solved in water!\n"); + } + // If deck has polymer, residual_ should contain polymer equation. + rq_.resize(fluid_.numPhases()+1); + residual_.material_balance_eq.resize(fluid_.numPhases()+1, ADB::null()); + assert(poly_pos_ == fluid_.numPhases()); + } + } + + template void @@ -266,7 +292,7 @@ namespace { template - void + int FullyImplicitBlackoilPolymerSolver:: step(const double dt, PolymerBlackoilState& x , @@ -277,27 +303,22 @@ namespace { // Initial max concentration of this time step from PolymerBlackoilState. cmax_ = Eigen::Map(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_)); - if (active_[Gas]) { updatePrimalVariableFromState(x); } - { - const SolutionState state = constantState(x, xw); - computeAccum(state, 0); - computeWellConnectionPressures(state, xw); - } - std::vector> residual_history; + // For each iteration we store in a vector the norms of the residual of + // the mass balance for each active phase, the well flux and the well equations + std::vector> residual_norms_history; - assemble(pvdt, x, xw, polymer_inflow); + assemble(pvdt, x, true, xw, polymer_inflow); bool converged = false; double omega = 1.; - residual_history.push_back(residuals()); + residual_norms_history.push_back(computeResidualNorms()); int it = 0; - converged = getConvergence(dt, it); - + converged = getConvergence(dt,it); const int sizeNonLinear = residual_.sizeNonLinear(); V dxOld = V::Zero(sizeNonLinear); @@ -305,37 +326,50 @@ namespace { bool isOscillate = false; bool isStagnate = false; const enum RelaxType relaxtype = relaxType(); + int linearIterations = 0; while ((!converged) && (it < maxIter())) { V dx = solveJacobianSystem(); - detectNewtonOscillations(residual_history, it, relaxRelTol(), isOscillate, isStagnate); + // store number of linear iterations used + linearIterations += linsolver_.iterations(); + + detectNewtonOscillations(residual_norms_history, it, relaxRelTol(), isOscillate, isStagnate); if (isOscillate) { omega -= relaxIncrement(); omega = std::max(omega, relaxMax()); - std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl; + if (terminal_output_) + { + std::cout << " Oscillating behavior detected: Relaxation set to " << omega << std::endl; + } } stablizeNewton(dx, dxOld, omega, relaxtype); updateState(dx, x, xw); - assemble(pvdt, x, xw, polymer_inflow); + assemble(pvdt, x, false, xw, polymer_inflow); - residual_history.push_back(residuals()); + residual_norms_history.push_back(computeResidualNorms()); + // increase iteration counter ++it; - converged = getConvergence(dt, it); + converged = getConvergence(dt,it); } if (!converged) { - std::cerr << "Failed to compute converged solution in " << it << " iterations. Ignoring!\n"; - // OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); + // the runtime_error is caught by the AdaptiveTimeStepping + OPM_THROW(std::runtime_error, "Failed to compute converged solution in " << it << " iterations."); + return -1; } + + linearIterations_ += linearIterations; + newtonIterations_ += it; // Update max concentration. computeCmax(x); + return linearIterations; } @@ -366,6 +400,7 @@ namespace { , concentration( ADB::null()) , qs ( ADB::null()) , bhp ( ADB::null()) + , canonical_phase_pressures(3, ADB::null()) { } @@ -375,30 +410,34 @@ namespace { template FullyImplicitBlackoilPolymerSolver:: - WellOps::WellOps(const Wells& wells) - : w2p(wells.well_connpos[ wells.number_of_wells ], - wells.number_of_wells) - , p2w(wells.number_of_wells, - wells.well_connpos[ wells.number_of_wells ]) + WellOps::WellOps(const Wells* wells) + : w2p(), + p2w() { - const int nw = wells.number_of_wells; - const int* const wpos = wells.well_connpos; + if( wells ) + { + w2p = M(wells->well_connpos[ wells->number_of_wells ], wells->number_of_wells); + p2w = M(wells->number_of_wells, wells->well_connpos[ wells->number_of_wells ]); - typedef Eigen::Triplet Tri; + const int nw = wells->number_of_wells; + const int* const wpos = wells->well_connpos; - std::vector scatter, gather; - scatter.reserve(wpos[nw]); - gather .reserve(wpos[nw]); + typedef Eigen::Triplet Tri; - for (int w = 0, i = 0; w < nw; ++w) { - for (; i < wpos[ w + 1 ]; ++i) { - scatter.push_back(Tri(i, w, 1.0)); - gather .push_back(Tri(w, i, 1.0)); + std::vector scatter, gather; + scatter.reserve(wpos[nw]); + gather .reserve(wpos[nw]); + + for (int w = 0, i = 0; w < nw; ++w) { + for (; i < wpos[ w + 1 ]; ++i) { + scatter.push_back(Tri(i, w, 1.0)); + gather .push_back(Tri(w, i, 1.0)); + } } - } - w2p.setFromTriplets(scatter.begin(), scatter.end()); - p2w.setFromTriplets(gather .begin(), gather .end()); + w2p.setFromTriplets(scatter.begin(), scatter.end()); + p2w.setFromTriplets(gather .begin(), gather .end()); + } } @@ -408,10 +447,21 @@ namespace { template typename FullyImplicitBlackoilPolymerSolver::SolutionState FullyImplicitBlackoilPolymerSolver::constantState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw) + const WellStateFullyImplicitBlackoil& xw) const { auto state = variableState(x, xw); + makeConstantState(state); + return state; + } + + + + + template + void + FullyImplicitBlackoilPolymerSolver::makeConstantState(SolutionState& state) const + { // HACK: throw away the derivatives. this may not be the most // performant way to do things, but it will make the state // automatically consistent with variableState() (and doing @@ -421,12 +471,17 @@ namespace { state.rs = ADB::constant(state.rs.value()); state.rv = ADB::constant(state.rv.value()); state.concentration = ADB::constant(state.concentration.value()); - for (int phaseIdx= 0; phaseIdx < x.numPhases(); ++ phaseIdx) + const int num_phases = state.saturation.size(); + for (int phaseIdx = 0; phaseIdx < num_phases; ++ phaseIdx) { state.saturation[phaseIdx] = ADB::constant(state.saturation[phaseIdx].value()); + } state.qs = ADB::constant(state.qs.value()); state.bhp = ADB::constant(state.bhp.value()); - - return state; + assert(state.canonical_phase_pressures.size() == Opm::BlackoilPhases::MaxNumPhases); + for (int canphase = 0; canphase < Opm::BlackoilPhases::MaxNumPhases; ++canphase) { + ADB& pp = state.canonical_phase_pressures[canphase]; + pp = ADB::constant(pp.value()); + } } @@ -436,14 +491,14 @@ namespace { template typename FullyImplicitBlackoilPolymerSolver::SolutionState FullyImplicitBlackoilPolymerSolver::variableState(const PolymerBlackoilState& x, - const WellStateFullyImplicitBlackoil& xw) + const WellStateFullyImplicitBlackoil& xw) const { using namespace Opm::AutoDiffGrid; const int nc = numCells(grid_); const int np = x.numPhases(); std::vector vars0; - // p, Sw and Rs, Rv or Sg, concentration are used as primary depending on solution conditions + // p, Sw and Rs, Rv or Sg, concentration is used as primary depending on solution conditions vars0.reserve(np + 2); // Initial pressure. assert (not x.pressure().empty()); @@ -492,7 +547,7 @@ namespace { xvar = isRs*rs + isRv*rv + isSg*sg; vars0.push_back(xvar); } - + // Initial polymer concentration. if (has_polymer_) { assert (not x.concentration().empty()); @@ -501,19 +556,25 @@ namespace { } // Initial well rates. - assert (not xw.wellRates().empty()); - // Need to reshuffle well rates, from phase running fastest - // to wells running fastest. - const int nw = wells_.number_of_wells; - // The transpose() below switches the ordering. - const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); - const V qs = Eigen::Map(wrates.data(), nw*np); - vars0.push_back(qs); + if( wellsActive() ) { + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + const int nw = wells().number_of_wells; - // Initial well bottom-hole pressure. - assert (not xw.bhp().empty()); - const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); - vars0.push_back(bhp); + // The transpose() below switches the ordering. + const DataBlock wrates = Eigen::Map(& xw.wellRates()[0], nw, np).transpose(); + const V qs = Eigen::Map(wrates.data(), nw*np); + vars0.push_back(qs); + + // Initial well bottom-hole pressure. + assert (not xw.bhp().empty()); + const V bhp = Eigen::Map(& xw.bhp()[0], xw.bhp().size()); + vars0.push_back(bhp); + } else { + // push null states for qs and bhp + vars0.push_back(V()); + vars0.push_back(V()); + } std::vector vars = ADB::variables(vars0); @@ -523,45 +584,54 @@ namespace { int nextvar = 0; state.pressure = vars[ nextvar++ ]; - // Temperature. - const V temp = Eigen::Map(& x.temperature()[0], nc, 1); + // temperature + const V temp = Eigen::Map(& x.temperature()[0], x.temperature().size()); state.temperature = ADB::constant(temp); // Saturations const std::vector& bpat = vars[0].blockPattern(); { - ADB so = ADB::constant(V::Ones(nc, 1), bpat); + + ADB so = ADB::constant(V::Ones(nc, 1), bpat); if (active_[ Water ]) { ADB& sw = vars[ nextvar++ ]; state.saturation[pu.phase_pos[ Water ]] = sw; - so = so - sw; + so -= sw; } - if (active_[ Gas]) { + if (active_[ Gas ]) { // Define Sg Rs and Rv in terms of xvar. + // Xvar is only defined if gas phase is active const ADB& xvar = vars[ nextvar++ ]; - const ADB& sg = isSg*xvar + isRv* so; - state.saturation[ pu.phase_pos[ Gas ] ] = sg; - so = so - sg; - const ADB rsSat = fluidRsSat(state.pressure, so, cells_); - const ADB rvSat = fluidRvSat(state.pressure, so, cells_); + ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ]; + sg = isSg*xvar + isRv* so; + so -= sg; - if (has_disgas_) { - state.rs = (1-isRs) * rsSat + isRs*xvar; - } else { - state.rs = rsSat; - } - if (has_vapoil_) { - state.rv = (1-isRv) * rvSat + isRv*xvar; - } else { - state.rv = rvSat; + if (active_[ Oil ]) { + // RS and RV is only defined if both oil and gas phase are active. + const ADB& sw = (active_[ Water ] + ? state.saturation[ pu.phase_pos[ Water ] ] + : ADB::constant(V::Zero(nc, 1), bpat)); + state.canonical_phase_pressures = computePressures(state.pressure, sw, so, sg); + const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_); + if (has_disgas_) { + state.rs = (1-isRs) * rsSat + isRs*xvar; + } else { + state.rs = rsSat; + } + const ADB rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_); + if (has_vapoil_) { + state.rv = (1-isRv) * rvSat + isRv*xvar; + } else { + state.rv = rvSat; + } } } if (active_[ Oil ]) { // Note that so is never a primary variable. - state.saturation[ pu.phase_pos[ Oil ] ] = so; + state.saturation[pu.phase_pos[ Oil ]] = so; } } @@ -569,6 +639,7 @@ namespace { if (has_polymer_) { state.concentration = vars[nextvar++]; } + // Qs. state.qs = vars[ nextvar++ ]; @@ -599,7 +670,6 @@ namespace { const ADB& c = state.concentration; const std::vector cond = phaseCondition(); - std::vector pressure = computePressures(state); const ADB pv_mult = poroMult(press); @@ -607,7 +677,7 @@ namespace { for (int phase = 0; phase < maxnp; ++phase) { if (active_[ phase ]) { const int pos = pu.phase_pos[ phase ]; - rq_[pos].b = fluidReciprocFVF(phase, pressure[phase], temp, rs, rv, cond, cells_); + 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]); @@ -619,10 +689,15 @@ namespace { const int po = pu.phase_pos[ Oil ]; const int pg = pu.phase_pos[ Gas ]; + // Temporary copy to avoid contribution of dissolved gas in the vaporized oil + // when both dissolved gas and vaporized oil are present. + const ADB accum_gas_copy =rq_[pg].accum[aix]; + rq_[pg].accum[aix] += state.rs * rq_[po].accum[aix]; - rq_[po].accum[aix] += state.rv * rq_[pg].accum[aix]; + rq_[po].accum[aix] += state.rv * accum_gas_copy; //DUMP(rq_[pg].accum[aix]); } + if (has_polymer_) { // compute polymer properties. const ADB cmax = ADB::constant(cmax_, state.concentration.blockPattern()); @@ -634,12 +709,14 @@ namespace { rq_[poly_pos_].accum[aix] = pv_mult * rq_[pu.phase_pos[Water]].b * sat[pu.phase_pos[Water]] * c * (1. - dead_pore_vol) + pv_mult * rho_rock * (1. - phi) / phi * ads; } + } + template void FullyImplicitBlackoilPolymerSolver::computeCmax(PolymerBlackoilState& state) { @@ -659,15 +736,34 @@ namespace { void FullyImplicitBlackoilPolymerSolver::computeWellConnectionPressures(const SolutionState& state, const WellStateFullyImplicitBlackoil& xw) { + if( ! wellsActive() ) return ; + using namespace Opm::AutoDiffGrid; // 1. Compute properties required by computeConnectionPressureDelta(). // Note that some of the complexity of this part is due to the function // taking std::vector arguments, and not Eigen objects. - const int nperf = wells_.well_connpos[wells_.number_of_wells]; - const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); - // Compute b, rsmax, rvmax values for perforations. - const ADB perf_press = subset(state.pressure, well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); + + // Compute the average pressure in each well block + const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); + V avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + } + + // Use cell values for the temperature as the wells don't knows its temperature yet. const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); std::vector perf_cond(nperf); const std::vector& pc = phaseCondition(); for (int perf = 0; perf < nperf; ++perf) { @@ -675,27 +771,27 @@ namespace { } const PhaseUsage& pu = fluid_.phaseUsage(); DataBlock b(nperf, pu.num_phases); - std::vector rssat_perf(nperf, 0.0); - std::vector rvsat_perf(nperf, 0.0); + std::vector rsmax_perf(nperf, 0.0); + std::vector rvmax_perf(nperf, 0.0); if (pu.phase_used[BlackoilPhases::Aqua]) { - const ADB bw = fluid_.bWat(perf_press, perf_temp, well_cells); - b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw.value(); + const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; } assert(active_[Oil]); - const ADB perf_so = subset(state.saturation[pu.phase_pos[Oil]], well_cells); + const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); if (pu.phase_used[BlackoilPhases::Liquid]) { const ADB perf_rs = subset(state.rs, well_cells); - const ADB bo = fluid_.bOil(perf_press, perf_temp, perf_rs, perf_cond, well_cells); - b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo.value(); - const V rssat = fluidRsSat(perf_press.value(), perf_so.value(), well_cells); - rssat_perf.assign(rssat.data(), rssat.data() + nperf); + const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + const V rssat = fluidRsSat(avg_press, perf_so, well_cells); + rsmax_perf.assign(rssat.data(), rssat.data() + nperf); } if (pu.phase_used[BlackoilPhases::Vapour]) { const ADB perf_rv = subset(state.rv, well_cells); - const ADB bg = fluid_.bGas(perf_press, perf_temp, perf_rv, perf_cond, well_cells); - b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg.value(); - const V rvsat = fluidRvSat(perf_press.value(), perf_so.value(), well_cells); - rvsat_perf.assign(rvsat.data(), rvsat.data() + nperf); + const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); + rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); } // b is row major, so can just copy data. std::vector b_perf(b.data(), b.data() + nperf * pu.num_phases); @@ -719,8 +815,8 @@ namespace { // 2. Compute pressure deltas, and store the results. std::vector cdp = WellDensitySegmented - ::computeConnectionPressureDelta(wells_, xw, fluid_.phaseUsage(), - b_perf, rssat_perf, rvsat_perf, perf_depth, + ::computeConnectionPressureDelta(wells(), xw, fluid_.phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, perf_depth, surf_dens, grav); well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); } @@ -734,6 +830,7 @@ namespace { FullyImplicitBlackoilPolymerSolver:: assemble(const V& pvdt, const PolymerBlackoilState& x , + const bool initial_assembly, WellStateFullyImplicitBlackoil& xw, const std::vector& polymer_inflow) { @@ -741,6 +838,16 @@ namespace { // Create the primary variables. SolutionState state = variableState(x, xw); + if (initial_assembly) { + // Create the (constant, derivativeless) initial state. + SolutionState state0 = state; + makeConstantState(state0); + // Compute initial accumulation contributions + // and well connection pressures. + computeAccum(state0, 0); + computeWellConnectionPressures(state0, xw); + } + // DISKVAL(state.pressure); // DISKVAL(state.saturation[0]); // DISKVAL(state.saturation[1]); @@ -757,27 +864,19 @@ namespace { // These quantities are stored in rq_[phase].accum[1]. // The corresponding accumulation terms from the start of // the timestep (b^0_p*s^0_p etc.) were already computed - // in step() and stored in rq_[phase].accum[0]. + // on the initial call to assemble() and stored in rq_[phase].accum[0]. computeAccum(state, 1); // Set up the common parts of the mass balance equations // for each active phase. const V transi = subset(geo_.transmissibility(), ops_.internal_faces); const std::vector kr = computeRelPerm(state); - const std::vector pressures = computePressures(state); - computeMassFlux(transi, kr, pressures, state); + computeMassFlux(transi, kr, state.canonical_phase_pressures, state); for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) { - // std::cout << "===== kr[" << phase << "] = \n" << std::endl; - // std::cout << kr[phase]; - // std::cout << "===== rq_[" << phase << "].mflux = \n" << std::endl; - // std::cout << rq_[phase].mflux; + // computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], state.canonical_phase_pressures[canph_[phaseIdx]], state); residual_.material_balance_eq[ phaseIdx ] = pvdt*(rq_[phaseIdx].accum[1] - rq_[phaseIdx].accum[0]) + ops_.div*rq_[phaseIdx].mflux; - - - // DUMP(ops_.div*rq_[phase].mflux); - // DUMP(residual_.material_balance_eq[phase]); } // -------- Extra (optional) rs and rv contributions to the mass balance equations -------- @@ -809,7 +908,6 @@ namespace { residual_.material_balance_eq[poly_pos_] = pvdt*(rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0]) + ops_.div*rq_[poly_pos_].mflux; } - // Note: updateWellControls() can change all its arguments if // a well control is switched. updateWellControls(state.bhp, state.qs, xw); @@ -828,13 +926,15 @@ namespace { V& aliveWells, const std::vector& polymer_inflow) { + if( ! wellsActive() ) return ; + const int nc = Opm::AutoDiffGrid::numCells(grid_); - const int np = wells_.number_of_phases; - const int nw = wells_.number_of_wells; - const int nperf = wells_.well_connpos[nw]; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - V Tw = Eigen::Map(wells_.WI, nperf); - const std::vector well_cells(wells_.well_cells, wells_.well_cells + nperf); + V Tw = Eigen::Map(wells().WI, nperf); + const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); // pressure diffs computed already (once per step, not changing per iteration) const V& cdp = well_perforation_pressure_diffs_; @@ -847,8 +947,13 @@ namespace { // DUMPVAL(state.bhp); // DUMPVAL(ADB::constant(cdp)); + // Perforation pressure + const ADB perfpressure = (wops_.w2p * state.bhp) + cdp; + std::vector perfpressure_d(perfpressure.value().data(), perfpressure.value().data() + nperf); + xw.perfPress() = perfpressure_d; + // Pressure drawdown (also used to determine direction of flow) - const ADB drawdown = p_perfcell - (wops_.w2p * state.bhp + cdp); + const ADB drawdown = p_perfcell - perfpressure; // current injecting connections auto connInjInx = drawdown.value() < 0; @@ -856,7 +961,7 @@ namespace { // injector == 1, producer == 0 V isInj = V::Zero(nw); for (int w = 0; w < nw; ++w) { - if (wells_.type[w] == INJECTOR) { + if (wells().type[w] == INJECTOR) { isInj[w] = 1; } } @@ -921,7 +1026,7 @@ namespace { } // compute avg. and total wellbore phase volumetric rates at std. conds - const DataBlock compi = Eigen::Map(wells_.comp_frac, nw, np); + const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); std::vector wbq(np, ADB::null()); ADB wbqt = ADB::constant(V::Zero(nw), state.pressure.blockPattern()); for (int phase = 0; phase < np; ++phase) { @@ -1008,6 +1113,7 @@ namespace { residual_.material_balance_eq[phase] -= superset(cq_s[phase],well_cells,nc); } + // Add well contributions to polymer mass balance equation if (has_polymer_) { const ADB mc = computeMc(state); @@ -1019,6 +1125,7 @@ namespace { well_cells, nc); } + // Add WELL EQUATIONS ADB qs = state.qs; for (int phase = 0; phase < np; ++phase) { @@ -1042,7 +1149,7 @@ namespace { - namespace + namespace detail { double rateToCompare(const ADB& well_phase_flow_rate, const int well, @@ -1115,25 +1222,27 @@ namespace { return broken; } - } // anonymous namespace + } // namespace detail template void FullyImplicitBlackoilPolymerSolver::updateWellControls(ADB& bhp, - ADB& well_phase_flow_rate, - WellStateFullyImplicitBlackoil& xw) const + ADB& well_phase_flow_rate, + WellStateFullyImplicitBlackoil& xw) const { + if( ! wellsActive() ) return ; + std::string modestring[3] = { "BHP", "RESERVOIR_RATE", "SURFACE_RATE" }; // Find, for each well, if any constraints are broken. If so, // switch control to first broken constraint. - const int np = wells_.number_of_phases; - const int nw = wells_.number_of_wells; + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; bool bhp_changed = false; bool rates_changed = false; for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells_.ctrls[w]; + const WellControls* wc = wells().ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. @@ -1150,16 +1259,19 @@ namespace { // inequality constraint, and therefore skipped. continue; } - if (constraintBroken(bhp, well_phase_flow_rate, w, np, wells_.type[w], wc, ctrl_index)) { + if (detail::constraintBroken(bhp, well_phase_flow_rate, w, np, wells().type[w], wc, ctrl_index)) { // ctrl_index will be the index of the broken constraint after the loop. break; } } if (ctrl_index != nwc) { // Constraint number ctrl_index was broken, switch to it. - std::cout << "Switching control mode for well " << wells_.name[w] - << " from " << modestring[well_controls_iget_type(wc, current)] - << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + if (terminal_output_) + { + std::cout << "Switching control mode for well " << wells().name[w] + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + } xw.currentControls()[w] = ctrl_index; // Also updating well state and primary variables. // We can only be switching to BHP and SURFACE_RATE @@ -1215,17 +1327,19 @@ namespace { template void FullyImplicitBlackoilPolymerSolver::addWellControlEq(const SolutionState& state, - const WellStateFullyImplicitBlackoil& xw, - const V& aliveWells) + const WellStateFullyImplicitBlackoil& xw, + const V& aliveWells) { - const int np = wells_.number_of_phases; - const int nw = wells_.number_of_wells; + if( ! wellsActive() ) return; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; V bhp_targets = V::Zero(nw); V rate_targets = V::Zero(nw); M rate_distr(nw, np*nw); for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells_.ctrls[w]; + const WellControls* wc = wells().ctrls[w]; // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. @@ -1292,11 +1406,20 @@ namespace { - namespace { - struct Chop01 { - double operator()(double x) const { return std::max(std::min(x, 1.0), 0.0); } - }; - } + namespace detail + { + + double infinityNorm( const ADB& a ) + { + if( a.value().size() > 0 ) { + return a.value().matrix().lpNorm (); + } + else { // this situation can occur when no wells are present + return 0.0; + } + } + + } // namespace detail @@ -1304,17 +1427,16 @@ namespace { template void FullyImplicitBlackoilPolymerSolver::updateState(const V& dx, - PolymerBlackoilState& state, - WellStateFullyImplicitBlackoil& well_state) + PolymerBlackoilState& state, + WellStateFullyImplicitBlackoil& well_state) { using namespace Opm::AutoDiffGrid; const int np = fluid_.numPhases(); const int nc = numCells(grid_); - const int nw = wells_.number_of_wells; + const int nw = wellsActive() ? wells().number_of_wells : 0; const V null; assert(null.size() == 0); const V zero = V::Zero(nc); - const V one = V::Constant(nc, 1.0); // store cell status in vectors V isRs = V::Zero(nc,1); @@ -1346,10 +1468,10 @@ namespace { const V dxvar = active_[Gas] ? subset(dx, Span(nc, 1, varstart)): null; varstart += dxvar.size(); - + const V dc = (has_polymer_) ? subset(dx, Span(nc, 1, varstart)) : null; varstart += dc.size(); - + const V dqs = subset(dx, Span(np*nw, 1, varstart)); varstart += dqs.size(); const V dbhp = subset(dx, Span(nw, 1, varstart)); @@ -1369,6 +1491,7 @@ namespace { const Opm::PhaseUsage& pu = fluid_.phaseUsage(); const DataBlock s_old = Eigen::Map(& state.saturation()[0], nc, np); const double dsmax = dsMax(); + V so; V sw; V sg; @@ -1461,104 +1584,110 @@ namespace { } // Update rs and rv - const double drsmaxrel = drsMaxRel(); - const double drvmax = 1e9;//% same as in Mrst + const double drmaxrel = drMaxRel(); V rs; if (has_disgas_) { const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); const V drs = isRs * dxvar; - const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drsmaxrel); + const V drs_limited = sign(drs) * drs.abs().min(rs_old.abs()*drmaxrel); rs = rs_old - drs_limited; } V rv; if (has_vapoil_) { const V rv_old = Eigen::Map(&state.rv()[0], nc); const V drv = isRv * dxvar; - const V drv_limited = sign(drv) * drv.abs().min(drvmax); + const V drv_limited = sign(drv) * drv.abs().min(rv_old.abs()*drmaxrel); rv = rv_old - drv_limited; } - // Update the state - if (has_disgas_) - std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin()); - - if (has_vapoil_) - std::copy(&rv[0], &rv[0] + nc, state.rv().begin()); // Sg is used as primal variable for water only cells. const double epsilon = std::sqrt(std::numeric_limits::epsilon()); auto watOnly = sw > (1 - epsilon); // phase translation sg <-> rs - const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_); - const V rsSat = fluidRsSat(p, so, cells_); - std::fill(primalVariable_.begin(), primalVariable_.end(), PrimalVariables::Sg); if (has_disgas_) { + const V rsSat0 = fluidRsSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_); + const V rsSat = fluidRsSat(p, so, cells_); // The obvious case auto hasGas = (sg > 0 && isRs == 0); - // keep oil saturated if previous sg is sufficient large: - const int pos = pu.phase_pos[ Gas ]; - auto hadGas = (sg <= 0 && s_old.col(pos) > epsilon); // Set oil saturated if previous rs is sufficiently large const V rs_old = Eigen::Map(&state.gasoilratio()[0], nc); auto gasVaporized = ( (rs > rsSat * (1+epsilon) && isRs == 1 ) && (rs_old > rsSat0 * (1-epsilon)) ); - - auto useSg = watOnly || hasGas || hadGas || gasVaporized; + auto useSg = watOnly || hasGas || gasVaporized; for (int c = 0; c < nc; ++c) { - if (useSg[c]) { rs[c] = rsSat[c];} - else { primalVariable_[c] = PrimalVariables::RS; } - + if (useSg[c]) { + rs[c] = rsSat[c]; + } else { + primalVariable_[c] = PrimalVariables::RS; + } } + } // phase transitions so <-> rv - const V rvSat0 = fluidRvSat(p_old, s_old.col(pu.phase_pos[Oil]), cells_); - const V rvSat = fluidRvSat(p, so, cells_); - if (has_vapoil_) { + + // The gas pressure is needed for the rvSat calculations + const V gaspress_old = computeGasPressure(p_old, s_old.col(Water), s_old.col(Oil), s_old.col(Gas)); + const V gaspress = computeGasPressure(p, sw, so, sg); + const V rvSat0 = fluidRvSat(gaspress_old, s_old.col(pu.phase_pos[Oil]), cells_); + const V rvSat = fluidRvSat(gaspress, so, cells_); + // The obvious case auto hasOil = (so > 0 && isRv == 0); - // keep oil saturated if previous so is sufficient large: - const int pos = pu.phase_pos[ Oil ]; - auto hadOil = (so <= 0 && s_old.col(pos) > epsilon ); // Set oil saturated if previous rv is sufficiently large const V rv_old = Eigen::Map(&state.rv()[0], nc); auto oilCondensed = ( (rv > rvSat * (1+epsilon) && isRv == 1) && (rv_old > rvSat0 * (1-epsilon)) ); - auto useSg = watOnly || hasOil || hadOil || oilCondensed; + auto useSg = watOnly || hasOil || oilCondensed; for (int c = 0; c < nc; ++c) { - if (useSg[c]) { rv[c] = rvSat[c]; } - else {primalVariable_[c] = PrimalVariables::RV; } - + if (useSg[c]) { + rv[c] = rvSat[c]; + } else { + primalVariable_[c] = PrimalVariables::RV; + } } } + // Update the state + if (has_disgas_) { + std::copy(&rs[0], &rs[0] + nc, state.gasoilratio().begin()); + } + + if (has_vapoil_) { + std::copy(&rv[0], &rv[0] + nc, state.rv().begin()); + } + //Polymer concentration updates. if (has_polymer_) { const V c_old = Eigen::Map(&state.concentration()[0], nc, 1); const V c = (c_old - dc).max(zero); std::copy(&c[0], &c[0] + nc, state.concentration().begin()); } + + if( wellsActive() ) + { + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const V dwr = Eigen::Map(wwr.data(), nw*np); + const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const V wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); - // Qs update. - // Since we need to update the wellrates, that are ordered by wells, - // from dqs which are ordered by phase, the simplest is to compute - // dwr, which is the data from dqs but ordered by wells. - const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); - const V dwr = Eigen::Map(wwr.data(), nw*np); - const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); - const V wr = wr_old - dwr; - std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); - - // Bhp update. - const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); - const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); - const V bhp = bhp_old - dbhp_limited; - std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + // Bhp update. + const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); + const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); + const V bhp = bhp_old - dbhp_limited; + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + } // Update phase conditions used for property calculations. updatePhaseCondFromPrimalVariable(); @@ -1595,6 +1724,9 @@ namespace { } + + + template std::vector FullyImplicitBlackoilPolymerSolver::computePressures(const SolutionState& state) const @@ -1606,18 +1738,34 @@ namespace { const ADB null = ADB::constant(V::Zero(nc, 1), bpat); const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - const ADB sw = (active_[ Water ] + const ADB& sw = (active_[ Water ] ? state.saturation[ pu.phase_pos[ Water ] ] : null); - const ADB so = (active_[ Oil ] + const ADB& so = (active_[ Oil ] ? state.saturation[ pu.phase_pos[ Oil ] ] : null); - const ADB sg = (active_[ Gas ] + const ADB& sg = (active_[ Gas ] ? state.saturation[ pu.phase_pos[ Gas ] ] : null); + return computePressures(state.pressure, sw, so, sg); + + } + + + + + + template + std::vector + FullyImplicitBlackoilPolymerSolver:: + computePressures(const ADB& po, + const ADB& sw, + const ADB& so, + const ADB& sg) const + { // convert the pressure offsets to the capillary pressures std::vector pressure = fluid_.capPress(sw, so, sg, cells_); for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) { @@ -1634,9 +1782,9 @@ namespace { // This is an unfortunate inconsistency, but a convention we must handle. for (int phaseIdx = 0; phaseIdx < BlackoilPhases::MaxNumPhases; ++phaseIdx) { if (phaseIdx == BlackoilPhases::Aqua) { - pressure[phaseIdx] = state.pressure - pressure[phaseIdx]; + pressure[phaseIdx] = po - pressure[phaseIdx]; } else { - pressure[phaseIdx] += state.pressure; + pressure[phaseIdx] += po; } } @@ -1646,14 +1794,34 @@ namespace { + + template + V + FullyImplicitBlackoilPolymerSolver::computeGasPressure(const V& po, + const V& sw, + const V& so, + const V& sg) const + { + assert (active_[Gas]); + std::vector cp = fluid_.capPress(ADB::constant(sw), + ADB::constant(so), + ADB::constant(sg), + cells_); + return cp[Gas].value() + po; + } + + + + + template std::vector FullyImplicitBlackoilPolymerSolver::computeRelPermWells(const SolutionState& state, - const DataBlock& well_s, - const std::vector& well_cells) const + const DataBlock& well_s, + const std::vector& well_cells) const { - const int nw = wells_.number_of_wells; - const int nperf = wells_.well_connpos[nw]; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; const std::vector& bpat = state.pressure.blockPattern(); const ADB null = ADB::constant(V::Zero(nperf), bpat); @@ -1680,20 +1848,21 @@ namespace { template void - FullyImplicitBlackoilPolymerSolver::computeMassFlux(const V& transi, + FullyImplicitBlackoilPolymerSolver::computeMassFlux(const V& transi, const std::vector& kr , const std::vector& phasePressure, const SolutionState& state) { - for (int phase = 0; phase < fluid_.numPhases(); ++phase) { + for (int phase = 0; phase < fluid_.numPhases(); ++phase) { const int canonicalPhaseIdx = canph_[phase]; const std::vector cond = phaseCondition(); const ADB tr_mult = transMult(state.pressure); - const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond, cells_); + const ADB mu = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv,cond, cells_); + rq_[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu; - const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond, cells_); + const ADB rho = fluidDensity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv,cond, cells_); ADB& head = rq_[phase].head; @@ -1731,11 +1900,13 @@ namespace { const ADB& b = rq_[phase].b; const ADB& mob = rq_[phase].mob; rq_[phase].mflux = upwind.select(b * mob) * head; - } + } } + + template void FullyImplicitBlackoilPolymerSolver::applyThresholdPressures(ADB& dp) @@ -1794,38 +1965,38 @@ namespace { template std::vector - FullyImplicitBlackoilPolymerSolver::residuals() const + FullyImplicitBlackoilPolymerSolver::computeResidualNorms() const { - std::vector residual; + std::vector residualNorms; std::vector::const_iterator massBalanceIt = residual_.material_balance_eq.begin(); const std::vector::const_iterator endMassBalanceIt = residual_.material_balance_eq.end(); for (; massBalanceIt != endMassBalanceIt; ++massBalanceIt) { - const double massBalanceResid = (*massBalanceIt).value().matrix().template lpNorm(); + const double massBalanceResid = detail::infinityNorm( (*massBalanceIt) ); if (!std::isfinite(massBalanceResid)) { OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual"); } - residual.push_back(massBalanceResid); + residualNorms.push_back(massBalanceResid); } // the following residuals are not used in the oscillation detection now - const double wellFluxResid = residual_.well_flux_eq.value().matrix().template lpNorm(); + const double wellFluxResid = detail::infinityNorm( residual_.well_flux_eq ); if (!std::isfinite(wellFluxResid)) { - OPM_THROW(Opm::NumericalProblem, + OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual"); } - residual.push_back(wellFluxResid); + residualNorms.push_back(wellFluxResid); - const double wellResid = residual_.well_eq.value().matrix().template lpNorm(); + const double wellResid = detail::infinityNorm( residual_.well_eq ); if (!std::isfinite(wellResid)) { OPM_THROW(Opm::NumericalProblem, "Encountered a non-finite residual"); } - residual.push_back(wellResid); + residualNorms.push_back(wellResid); - return residual; + return residualNorms; } template @@ -1845,6 +2016,7 @@ namespace { return; } + stagnate = true; int oscillatePhase = 0; const std::vector& F0 = residual_history[it]; const std::vector& F1 = residual_history[it - 1]; @@ -1855,7 +2027,9 @@ namespace { oscillatePhase += (d1 < relaxRelTol) && (relaxRelTol < d2); - stagnate = ! (std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3); + // Process is 'stagnate' unless at least one phase + // exhibits significant residual change. + stagnate = (stagnate && !(std::abs((F1[p] - F2[p]) / F2[p]) > 1.0e-3)); } oscillate = (oscillatePhase > 1); @@ -1865,7 +2039,7 @@ namespace { template void FullyImplicitBlackoilPolymerSolver::stablizeNewton(V& dx, V& dxOld, const double omega, - const RelaxType relax_type) const + const RelaxType relax_type) const { // The dxOld is updated with dx. // If omega is equal to 1., no relaxtion will be appiled. @@ -1893,122 +2067,177 @@ namespace { return; } + template + double + FullyImplicitBlackoilPolymerSolver::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 + { + // Do the global reductions +#if HAVE_MPI + if ( linsolver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast(linsolver_.parallelInformation()); + + // Compute the global number of cells and porevolume + std::vector v(nc, 1); + auto nc_and_pv = std::tuple(0, 0.0); + auto nc_and_pv_operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), + Opm::Reduction::makeGlobalSumFunctor()); + 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); + auto containers = std::make_tuple(B.col(idx), + tempV.col(idx), + R.col(idx)); + auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), + 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); + } + else + { + R_sum[idx] = B_avg[idx] = maxCoeff[idx] = 0.0; + } + } + // Compute pore volume + return std::get<1>(nc_and_pv); + } + else +#endif + { + for ( int idx=0; idx bool FullyImplicitBlackoilPolymerSolver::getConvergence(const double dt, const int iteration) { - const double tol_mb = 1.0e-7; - const double tol_cnv = 1.0e-3; + 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 Opm::PhaseUsage& pu = fluid_.phaseUsage(); const V pv = geo_.poreVolume(); - const double pvSum = pv.sum(); const std::vector cond = phaseCondition(); - double CNVW = 0.; - double CNVO = 0.; - double CNVG = 0.; - double CNVP = 0.; + std::array CNV = {{0., 0., 0., 0.}}; + std::array R_sum = {{0., 0., 0., 0.}}; + 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::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); - double RW_sum = 0.; - double RO_sum = 0.; - double RG_sum = 0.; - double RP_sum = 0.; - - double BW_avg = 0.; - double BO_avg = 0.; - double BG_avg = 0.; - - if (active_[Water]) { - const int pos = pu.phase_pos[Water]; - const ADB& tempBW = rq_[pos].b; - V BW = 1./tempBW.value(); - V RW = residual_.material_balance_eq[pos].value(); - BW_avg = BW.sum()/nc; - const V tempV = RW.abs()/pv; - - CNVW = BW_avg * dt * tempV.maxCoeff(); - RW_sum = RW.sum(); + for ( int idx=0; idx(); - double residualWell = residual_.well_eq.value().matrix().template lpNorm(); - - bool converged_Well = (residualWellFlux < 1./Opm::unit::day) && (residualWell < Opm::unit::barsa); - - bool converged = converged_MB && converged_CNV && converged_Well; - - if (iteration == 0) { - std::cout << "\nIter MB(OIL) MB(WATER) MB(GAS) MB(POLYMER) CNVW CNVO CNVG CNVP WELL-FLOW WELL-CNTRL\n"; + bool converged_MB = true; + bool converged_CNV = true; + // Finish computation + 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(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(residualWell) || residualWell > maxResidualAllowed() ) + { + OPM_THROW(Opm::NumericalProblem,"One of the residuals is NaN or to large!"); + } + + if ( terminal_output_ ) + { + // 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"; + } + 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) << mass_balance_residual[MaxNumPhases] + << std::setw(11) << CNV[Water] + << std::setw(11) << CNV[Oil] + << std::setw(11) << CNV[Gas] + << std::setw(11) << CNV[MaxNumPhases] + << std::setw(11) << residualWellFlux + << std::setw(11) << residualWell + << std::endl; + std::cout.precision(oprec); + std::cout.flags(oflags); } - 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) << mass_balance_residual_polymer - << std::setw(11) << CNVW - << std::setw(11) << CNVO - << std::setw(11) << CNVG - << std::setw(11) << CNVP - << std::setw(11) << residualWellFlux - << std::setw(11) << residualWell - << std::endl; - std::cout.precision(oprec); - std::cout.flags(oflags); return converged; } @@ -2117,6 +2346,7 @@ namespace { return fluid_.rsSat(p, satOil, cells); } + template V FullyImplicitBlackoilPolymerSolver::fluidRvSat(const V& p, @@ -2139,10 +2369,6 @@ namespace { return fluid_.rvSat(p, satOil, cells); } - - - - template ADB FullyImplicitBlackoilPolymerSolver::computeMc(const SolutionState& state) const @@ -2152,7 +2378,6 @@ namespace { - template ADB FullyImplicitBlackoilPolymerSolver::poroMult(const ADB& p) const @@ -2169,10 +2394,10 @@ namespace { const int num_blocks = p.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { - jacs[block] = dpm_diag * p.derivative()[block]; + fastSparseProduct(dpm_diag, p.derivative()[block], jacs[block]); } return ADB::function(pm, jacs); - } else { + } else { return ADB::constant(V::Constant(n, 1.0), p.blockPattern()); } } @@ -2197,7 +2422,7 @@ namespace { const int num_blocks = p.numBlocks(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { - jacs[block] = dtm_diag * p.derivative()[block]; + fastSparseProduct(dtm_diag, p.derivative()[block], jacs[block]); } return ADB::function(tm, jacs); } else { @@ -2209,7 +2434,7 @@ namespace { /* template void - FullyImplicitBlackoilPolymerSolver:: + FullyImplicitBlackoilSolver:: classifyCondition(const SolutionState& state, std::vector& cond ) const { @@ -2316,7 +2541,7 @@ namespace { // For oil only cells Rs is used as primal variable. For cells almost full of water // the default primal variable (Sg) is used. - if (has_disgas_) { + if (has_disgas_) { for (V::Index c = 0, e = sg.size(); c != e; ++c) { if ( !watOnly[c] && hasOil[c] && !hasGas[c] ) {primalVariable_[c] = PrimalVariables::RS; } } diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp index 8b165d8a2..78caf9c26 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp @@ -1,5 +1,5 @@ /* - Copyright 2014 SINTEF ICT, Applied Mathematics. + Copyright 2013 SINTEF ICT, Applied Mathematics. Copyright 2014 STATOIL ASA. This file is part of the Open Porous Media project (OPM). @@ -26,6 +26,7 @@ struct UnstructuredGrid; struct Wells; +struct FlowBoundaryConditions; namespace Opm { @@ -37,9 +38,8 @@ namespace Opm class SimulatorTimer; class PolymerBlackoilState; class WellStateFullyImplicitBlackoil; - class WellsManager; class EclipseState; - class EclipseWriter; + class BlackoilOutputWriter; class PolymerPropsAd; class PolymerInflowInterface; struct SimulatorReport; @@ -70,7 +70,6 @@ namespace Opm /// \param[in] grid grid data structure /// \param[in] geo derived geological properties /// \param[in] props fluid and rock properties - /// \param[in] polymer_props polymer properties /// \param[in] rock_comp_props if non-null, rock compressibility properties /// \param[in] linsolver linear solver /// \param[in] gravity if non-null, gravity vector @@ -91,7 +90,7 @@ namespace Opm const bool vapoil, const bool polymer, std::shared_ptr eclipse_state, - EclipseWriter& output_writer, + BlackoilOutputWriter& output_writer, Opm::DeckConstPtr& deck, const std::vector& threshold_pressures_by_face); diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp deleted file mode 100644 index cb3d710db..000000000 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.cpp +++ /dev/null @@ -1,211 +0,0 @@ -#include "config.h" - -#include "SimulatorFullyImplicitBlackoilPolymerOutput.hpp" - -#include -#include -#include -#include - -#include - -#include -#include -#include - -#include - -#ifdef HAVE_DUNE_CORNERPOINT -#include -#include -#include -#include -#endif -namespace Opm -{ - - - void outputStateVtk(const UnstructuredGrid& grid, - const PolymerBlackoilState& state, - const int step, - const std::string& output_dir) - { - // Write data in VTK format. - std::ostringstream vtkfilename; - vtkfilename << output_dir << "/vtk_files"; - boost::filesystem::path fpath(vtkfilename.str()); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - vtkfilename << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; - std::ofstream vtkfile(vtkfilename.str().c_str()); - if (!vtkfile) { - OPM_THROW(std::runtime_error, "Failed to open " << vtkfilename.str()); - } - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - dm["concentration"] = &state.concentration(); - dm["maxconcentration"] = &state.maxconcentration(); - std::vector cell_velocity; - Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), - AutoDiffGrid::numFaces(grid), - AutoDiffGrid::beginFaceCentroids(grid), - AutoDiffGrid::faceCells(grid), - AutoDiffGrid::beginCellCentroids(grid), - AutoDiffGrid::beginCellVolumes(grid), - AutoDiffGrid::dimensions(grid), - state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - Opm::writeVtkData(grid, dm, vtkfile); - } - - - void outputStateMatlab(const UnstructuredGrid& grid, - const PolymerBlackoilState& state, - const int step, - const std::string& output_dir) - { - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - dm["surfvolume"] = &state.surfacevol(); - dm["rs"] = &state.gasoilratio(); - dm["rv"] = &state.rv(); - dm["concentration"] = &state.concentration(); - dm["maxconcentration"] = &state.maxconcentration(); - std::vector cell_velocity; - Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), - AutoDiffGrid::numFaces(grid), - AutoDiffGrid::beginFaceCentroids(grid), - UgGridHelpers::faceCells(grid), - AutoDiffGrid::beginCellCentroids(grid), - AutoDiffGrid::beginCellVolumes(grid), - AutoDiffGrid::dimensions(grid), - state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); - } - file.precision(15); - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } - } - void outputWellStateMatlab(const Opm::WellState& well_state, - const int step, - const std::string& output_dir) - { - Opm::DataMap dm; - dm["bhp"] = &well_state.bhp(); - dm["wellrates"] = &well_state.wellRates(); - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - OPM_THROW(std::runtime_error,"Failed to open " << fname.str()); - } - file.precision(15); - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } - } - -#if 0 - void outputWaterCut(const Opm::Watercut& watercut, - const std::string& output_dir) - { - // Write water cut curve. - std::string fname = output_dir + "/watercut.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - OPM_THROW(std::runtime_error, "Failed to open " << fname); - } - watercut.write(os); - } - - void outputWellReport(const Opm::WellReport& wellreport, - const std::string& output_dir) - { - // Write well report. - std::string fname = output_dir + "/wellreport.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - OPM_THROW(std::runtime_error, "Failed to open " << fname); - } - wellreport.write(os); - } -#endif - -#ifdef HAVE_DUNE_CORNERPOINT - void outputStateVtk(const Dune::CpGrid& grid, - const PolymerBlackoilState& state, - const int step, - const std::string& output_dir) - { - // Write data in VTK format. - std::ostringstream vtkfilename; - std::ostringstream vtkpath; - vtkpath << output_dir << "/vtk_files"; - vtkpath << "/output-" << std::setw(3) << std::setfill('0') << step; - boost::filesystem::path fpath(vtkpath.str()); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - vtkfilename << "output-" << std::setw(3) << std::setfill('0') << step; -#if DUNE_VERSION_NEWER(DUNE_GRID, 2, 3) - Dune::VTKWriter writer(grid.leafGridView(), Dune::VTK::nonconforming); -#else - Dune::VTKWriter writer(grid.leafView(), Dune::VTK::nonconforming); -#endif - writer.addCellData(state.saturation(), "saturation", state.numPhases()); - writer.addCellData(state.pressure(), "pressure", 1); - writer.addCellData(state.concentration(), "concentration", 1); - writer.addCellData(state.maxconcentration(), "maxconcentration", 1); - - std::vector cell_velocity; - Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), - AutoDiffGrid::numFaces(grid), - AutoDiffGrid::beginFaceCentroids(grid), - AutoDiffGrid::faceCells(grid), - AutoDiffGrid::beginCellCentroids(grid), - AutoDiffGrid::beginCellVolumes(grid), - AutoDiffGrid::dimensions(grid), - state.faceflux(), cell_velocity); - writer.addCellData(cell_velocity, "velocity", Dune::CpGrid::dimension); - writer.pwrite(vtkfilename.str(), vtkpath.str(), std::string("."), Dune::VTK::ascii); - } -#endif - -} diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp deleted file mode 100644 index 55a364324..000000000 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymerOutput.hpp +++ /dev/null @@ -1,96 +0,0 @@ -#ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED -#define OPM_SIMULATORFULLYIMPLICITBLACKOILPOLYMEROUTPUT_HEADER_INCLUDED -#include -#include -#include -#include -#include -#include - -#include - -#include -#include -#include -#include - -#include - -#ifdef HAVE_DUNE_CORNERPOINT -#include -#endif -namespace Opm -{ - - void outputStateVtk(const UnstructuredGrid& grid, - const PolymerBlackoilState& state, - const int step, - const std::string& output_dir); - - - void outputStateMatlab(const UnstructuredGrid& grid, - const PolymerBlackoilState& state, - const int step, - const std::string& output_dir); - - void outputWellStateMatlab(const Opm::WellState& well_state, - const int step, - const std::string& output_dir); -#ifdef HAVE_DUNE_CORNERPOINT - void outputStateVtk(const Dune::CpGrid& grid, - const PolymerBlackoilState& state, - const int step, - const std::string& output_dir); -#endif - - template - void outputStateMatlab(const Grid& grid, - const PolymerBlackoilState& state, - const int step, - const std::string& output_dir) - { - Opm::DataMap dm; - dm["saturation"] = &state.saturation(); - dm["pressure"] = &state.pressure(); - dm["surfvolume"] = &state.surfacevol(); - dm["rs"] = &state.gasoilratio(); - dm["rv"] = &state.rv(); - dm["concentration"] = &state.concentration(); - dm["maxconcentration"] = &state.maxconcentration(); - - std::vector cell_velocity; - Opm::estimateCellVelocity(AutoDiffGrid::numCells(grid), - AutoDiffGrid::numFaces(grid), - AutoDiffGrid::beginFaceCentroids(grid), - UgGridHelpers::faceCells(grid), - AutoDiffGrid::beginCellCentroids(grid), - AutoDiffGrid::beginCellVolumes(grid), - AutoDiffGrid::dimensions(grid), - state.faceflux(), cell_velocity); - dm["velocity"] = &cell_velocity; - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - OPM_THROW(std::runtime_error, "Failed to open " << fname.str()); - } - file.precision(15); - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } - } - - -} -#endif diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp index 0f9af3844..33cf621af 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer_impl.hpp @@ -1,5 +1,6 @@ /* - Copyright 2014 SINTEF ICT, Applied Mathematics. + Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2014 IRIS AS Copyright 2014 STATOIL ASA. This file is part of the Open Porous Media project (OPM). @@ -18,11 +19,13 @@ along with OPM. If not, see . */ -#include +//#include +#include #include #include #include #include + #include #include @@ -36,9 +39,9 @@ #include #include -#include #include #include +//#include #include #include #include @@ -46,6 +49,7 @@ #include +//#include #include #include @@ -80,7 +84,7 @@ namespace Opm const Grid& grid, const DerivedGeology& geo, BlackoilPropsAdInterface& props, - const PolymerPropsAd& polymer_props, + const PolymerPropsAd& polymer_props, const RockCompressibility* rock_comp_props, NewtonIterationBlackoilInterface& linsolver, const double* gravity, @@ -88,7 +92,7 @@ namespace Opm bool has_vapoil, bool has_polymer, std::shared_ptr eclipse_state, - EclipseWriter& output_writer, + BlackoilOutputWriter& output_writer, Opm::DeckConstPtr& deck, const std::vector& threshold_pressures_by_face); @@ -103,11 +107,6 @@ namespace Opm const parameter::ParameterGroup param_; - // Parameters for output. - bool output_; - bool output_vtk_; - std::string output_dir_; - int output_interval_; // Observed objects. const Grid& grid_; BlackoilPropsAdInterface& props_; @@ -122,10 +121,11 @@ namespace Opm const bool has_disgas_; const bool has_vapoil_; const bool has_polymer_; + bool terminal_output_; // eclipse_state std::shared_ptr eclipse_state_; // output_writer - EclipseWriter& output_writer_; + BlackoilOutputWriter& output_writer_; Opm::DeckConstPtr& deck_; RateConverterType rateConverter_; // Threshold pressures. @@ -134,7 +134,7 @@ namespace Opm void computeRESV(const std::size_t step, const Wells* wells, - const PolymerBlackoilState& x, + const BlackoilState& x, WellStateFullyImplicitBlackoil& xw); }; @@ -154,13 +154,13 @@ namespace Opm const bool has_vapoil, const bool has_polymer, std::shared_ptr eclipse_state, - EclipseWriter& output_writer, + BlackoilOutputWriter& output_writer, Opm::DeckConstPtr& deck, const std::vector& threshold_pressures_by_face) { - pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, - has_vapoil, has_polymer, eclipse_state, output_writer, deck, threshold_pressures_by_face)); + pimpl_.reset(new Impl(param, grid, geo, props, polymer_props, rock_comp_props, linsolver, gravity, has_disgas, has_vapoil, has_polymer, + eclipse_state, output_writer, deck, threshold_pressures_by_face)); } @@ -175,64 +175,6 @@ namespace Opm } - - static void outputWellStateMatlab(const Opm::WellStateFullyImplicitBlackoil& well_state, - const int step, - const std::string& output_dir) - { - Opm::DataMap dm; - dm["bhp"] = &well_state.bhp(); - dm["wellrates"] = &well_state.wellRates(); - - // Write data (not grid) in Matlab format - for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { - std::ostringstream fname; - fname << output_dir << "/" << it->first; - boost::filesystem::path fpath = fname.str(); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error,"Creating directories failed: " << fpath); - } - fname << "/" << std::setw(3) << std::setfill('0') << step << ".txt"; - std::ofstream file(fname.str().c_str()); - if (!file) { - OPM_THROW(std::runtime_error,"Failed to open " << fname.str()); - } - file.precision(15); - const std::vector& d = *(it->second); - std::copy(d.begin(), d.end(), std::ostream_iterator(file, "\n")); - } - } - -#if 0 - static void outputWaterCut(const Opm::Watercut& watercut, - const std::string& output_dir) - { - // Write water cut curve. - std::string fname = output_dir + "/watercut.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - OPM_THROW(std::runtime_error, "Failed to open " << fname); - } - watercut.write(os); - } - - static void outputWellReport(const Opm::WellReport& wellreport, - const std::string& output_dir) - { - // Write well report. - std::string fname = output_dir + "/wellreport.txt"; - std::ofstream os(fname.c_str()); - if (!os) { - OPM_THROW(std::runtime_error, "Failed to open " << fname); - } - wellreport.write(os); - } -#endif - - // \TODO: Treat bcs. template SimulatorFullyImplicitBlackoilPolymer::Impl::Impl(const parameter::ParameterGroup& param, @@ -247,7 +189,7 @@ namespace Opm const bool has_vapoil, const bool has_polymer, std::shared_ptr eclipse_state, - EclipseWriter& output_writer, + BlackoilOutputWriter& output_writer, Opm::DeckConstPtr& deck, const std::vector& threshold_pressures_by_face) : param_(param), @@ -261,34 +203,30 @@ namespace Opm has_disgas_(has_disgas), has_vapoil_(has_vapoil), has_polymer_(has_polymer), + terminal_output_(param.getDefault("output_terminal", true)), eclipse_state_(eclipse_state), output_writer_(output_writer), deck_(deck), rateConverter_(props_, std::vector(AutoDiffGrid::numCells(grid_), 0)), threshold_pressures_by_face_(threshold_pressures_by_face) { - // For output. - output_ = param.getDefault("output", true); - if (output_) { - output_vtk_ = param.getDefault("output_vtk", true); - output_dir_ = param.getDefault("output_dir", std::string("output")); - // Ensure that output dir exists - boost::filesystem::path fpath(output_dir_); - try { - create_directories(fpath); - } - catch (...) { - OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath); - } - output_interval_ = param.getDefault("output_interval", 1); - } - // Misc init. const int num_cells = AutoDiffGrid::numCells(grid); allcells_.resize(num_cells); for (int cell = 0; cell < num_cells; ++cell) { allcells_[cell] = cell; } +#if HAVE_MPI + if ( terminal_output_ ) { + if ( solver_.parallelInformation().type() == typeid(ParallelISTLInformation) ) + { + const ParallelISTLInformation& info = + boost::any_cast(solver_.parallelInformation()); + // Only rank 0 does print to std::cout + terminal_output_= (info.communicator().rank()==0); + } + } +#endif } @@ -306,14 +244,40 @@ namespace Opm Opm::time::StopWatch step_timer; Opm::time::StopWatch total_timer; total_timer.start(); - std::string tstep_filename = output_dir_ + "/step_timing.txt"; + std::string tstep_filename = output_writer_.outputDirectory() + "/step_timing.txt"; std::ofstream tstep_os(tstep_filename.c_str()); + typename FullyImplicitBlackoilPolymerSolver::SolverParameter solverParam( param_ ); + + //adaptive time stepping + // std::unique_ptr< AdaptiveTimeStepping > adaptiveTimeStepping; + // if( param_.getDefault("timestep.adaptive", bool(false) ) ) + // { + // adaptiveTimeStepping.reset( new AdaptiveTimeStepping( param_ ) ); + // } + + // init output writer + output_writer_.writeInit( timer ); + + std::string restorefilename = param_.getDefault("restorefile", std::string("") ); + if( ! restorefilename.empty() ) + { + // -1 means that we'll take the last report step that was written + const int desiredRestoreStep = param_.getDefault("restorestep", int(-1) ); + output_writer_.restore( timer, state.blackoilState(), prev_well_state, restorefilename, desiredRestoreStep ); + } + + unsigned int totalNewtonIterations = 0; + unsigned int totalLinearIterations = 0; + // Main simulation loop. while (!timer.done()) { // Report timestep. step_timer.start(); - timer.report(std::cout); + if ( terminal_output_ ) + { + timer.report(std::cout); + } // Create wells and well state. WellsManager wells_manager(eclipse_state_, @@ -328,6 +292,7 @@ namespace Opm const Wells* wells = wells_manager.c_wells(); WellStateFullyImplicitBlackoil well_state; well_state.init(wells, state.blackoilState(), prev_well_state); + // compute polymer inflow std::unique_ptr polymer_inflow_ptr; if (deck_->hasKeyword("WPOLYMER")) { @@ -344,42 +309,54 @@ namespace Opm polymer_inflow_ptr->getInflowValues(timer.simulationTimeElapsed(), timer.simulationTimeElapsed() + timer.currentStepLength(), polymer_inflow_c); - // Output state at start of time step. - if (output_ && (timer.currentStepNum() % output_interval_ == 0)) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputWellStateMatlab(well_state,timer.currentStepNum(), output_dir_); - } - if (output_) { - if (timer.currentStepNum() == 0) { - output_writer_.writeInit(timer); - } - output_writer_.writeTimeStep(timer, state.blackoilState(), well_state); - } + + // write simulation state at the report stage + output_writer_.writeTimeStep( timer, state.blackoilState(), well_state ); // Max oil saturation (for VPPARS), hysteresis update. props_.updateSatOilMax(state.saturation()); props_.updateSatHyst(state.saturation(), allcells_); // Compute reservoir volumes for RESV controls. - computeRESV(timer.currentStepNum(), wells, state, well_state); + computeRESV(timer.currentStepNum(), wells, state.blackoilState(), well_state); - // Run a single step of the solver. + // Run a multiple steps of the solver depending on the time step control. solver_timer.start(); - FullyImplicitBlackoilPolymerSolver solver(param_, grid_, props_, geo_, rock_comp_props_, polymer_props_, *wells, solver_, has_disgas_, has_vapoil_, has_polymer_); + + FullyImplicitBlackoilPolymerSolver solver(solverParam, grid_, props_, geo_, rock_comp_props_, polymer_props_, wells, solver_, has_disgas_, has_vapoil_, has_polymer_, terminal_output_); if (!threshold_pressures_by_face_.empty()) { solver.setThresholdPressures(threshold_pressures_by_face_); } + + // If sub stepping is enabled allow the solver to sub cycle + // in case the report steps are to large for the solver to converge + // + // \Note: The report steps are met in any case + // \Note: The sub stepping will require a copy of the state variables + // if( adaptiveTimeStepping ) { + // adaptiveTimeStepping->step( timer, solver, state, well_state, output_writer_ ); + // } else { + // solve for complete report step solver.step(timer.currentStepLength(), state, well_state, polymer_inflow_c); + // } + + // take time that was used to solve system for this reportStep solver_timer.stop(); + // accumulate the number of Newton and Linear Iterations + totalNewtonIterations += solver.newtonIterations(); + totalLinearIterations += solver.linearIterations(); + // Report timing. const double st = solver_timer.secsSinceStart(); - std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; + + if ( terminal_output_ ) + { + std::cout << "Fully implicit solver took: " << st << " seconds." << std::endl; + } + stime += st; - if (output_) { + if ( output_writer_.output() ) { SimulatorReport step_report; step_report.pressure_time = st; step_report.total_time = step_timer.secsSinceStart(); @@ -392,14 +369,7 @@ namespace Opm } // Write final simulation state. - if (output_) { - if (output_vtk_) { - outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_); - } - outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_); - outputWellStateMatlab(prev_well_state, timer.currentStepNum(), output_dir_); - output_writer_.writeTimeStep(timer, state.blackoilState(), prev_well_state); - } + output_writer_.writeTimeStep( timer, state.blackoilState(), prev_well_state ); // Stop timer and create timing report total_timer.stop(); @@ -407,6 +377,8 @@ namespace Opm report.pressure_time = stime; report.transport_time = 0.0; report.total_time = total_timer.secsSinceStart(); + report.total_newton_iterations = totalNewtonIterations; + report.total_linear_iterations = totalLinearIterations; return report; } @@ -444,17 +416,16 @@ namespace Opm } inline bool - is_resv_prod(const Wells& wells, - const int w) + is_resv(const Wells& wells, + const int w) { - return ((wells.type[w] == PRODUCER) && - (0 <= resv_control(wells.ctrls[w]))); + return (0 <= resv_control(wells.ctrls[w])); } inline bool - is_resv_prod(const WellMap& wmap, - const std::string& name, - const std::size_t step) + is_resv(const WellMap& wmap, + const std::string& name, + const std::size_t step) { bool match = false; @@ -465,29 +436,34 @@ namespace Opm match = (wp->isProducer(step) && wp->getProductionProperties(step) - .hasProductionControl(WellProducer::RESV)); + .hasProductionControl(WellProducer::RESV)) + || (wp->isInjector(step) && + wp->getInjectionProperties(step) + .hasInjectionControl(WellInjector::RESV)); } return match; } inline std::vector - resvProducers(const Wells& wells, - const std::size_t step, - const WellMap& wmap) + resvWells(const Wells* wells, + const std::size_t step, + const WellMap& wmap) { - std::vector resv_prod; - - for (int w = 0, nw = wells.number_of_wells; w < nw; ++w) { - if (is_resv_prod(wells, w) || - ((wells.name[w] != 0) && - is_resv_prod(wmap, wells.name[w], step))) - { - resv_prod.push_back(w); + std::vector resv_wells; + if( wells ) + { + for (int w = 0, nw = wells->number_of_wells; w < nw; ++w) { + if (is_resv(*wells, w) || + ((wells->name[w] != 0) && + is_resv(wmap, wells->name[w], step))) + { + resv_wells.push_back(w); + } } } - return resv_prod; + return resv_wells; } inline void @@ -527,7 +503,7 @@ namespace Opm SimulatorFullyImplicitBlackoilPolymer:: Impl::computeRESV(const std::size_t step, const Wells* wells, - const PolymerBlackoilState& x, + const BlackoilState& x, WellStateFullyImplicitBlackoil& xw) { typedef SimFIBODetails::WellMap WellMap; @@ -535,24 +511,24 @@ namespace Opm const std::vector& w_ecl = eclipse_state_->getSchedule()->getWells(step); const WellMap& wmap = SimFIBODetails::mapWells(w_ecl); - const std::vector& resv_prod = - SimFIBODetails::resvProducers(*wells, step, wmap); + const std::vector& resv_wells = SimFIBODetails::resvWells(wells, step, wmap); - if (! resv_prod.empty()) { + if (! resv_wells.empty()) { const PhaseUsage& pu = props_.phaseUsage(); const std::vector::size_type np = props_.numPhases(); - rateConverter_.defineState(x.blackoilState()); + rateConverter_.defineState(x); std::vector distr (np); std::vector hrates(np); std::vector prates(np); for (std::vector::const_iterator - rp = resv_prod.begin(), e = resv_prod.end(); + rp = resv_wells.begin(), e = resv_wells.end(); rp != e; ++rp) { WellControls* ctrl = wells->ctrls[*rp]; + const bool is_producer = wells->type[*rp] == PRODUCER; // RESV control mode, all wells { @@ -561,11 +537,17 @@ namespace Opm if (0 <= rctrl) { const std::vector::size_type off = (*rp) * np; - // Convert to positive rates to avoid issues - // in coefficient calculations. - std::transform(xw.wellRates().begin() + (off + 0*np), - xw.wellRates().begin() + (off + 1*np), - prates.begin(), std::negate()); + if (is_producer) { + // Convert to positive rates to avoid issues + // in coefficient calculations. + std::transform(xw.wellRates().begin() + (off + 0*np), + xw.wellRates().begin() + (off + 1*np), + prates.begin(), std::negate()); + } else { + std::copy(xw.wellRates().begin() + (off + 0*np), + xw.wellRates().begin() + (off + 1*np), + prates.begin()); + } const int fipreg = 0; // Hack. Ignore FIP regions. rateConverter_.calcCoeff(prates, fipreg, distr); @@ -576,7 +558,7 @@ namespace Opm // RESV control, WCONHIST wells. A bit of duplicate // work, regrettably. - if (wells->name[*rp] != 0) { + if (is_producer && wells->name[*rp] != 0) { WellMap::const_iterator i = wmap.find(wells->name[*rp]); if (i != wmap.end()) { @@ -603,11 +585,18 @@ namespace Opm well_controls_clear(ctrl); well_controls_assert_number_of_phases(ctrl, int(np)); - const int ok = + const int ok_resv = well_controls_add_new(RESERVOIR_RATE, target, & distr[0], ctrl); - if (ok != 0) { + // For WCONHIST/RESV the BHP limit is set to 1 atm. + // TODO: Make it possible to modify the BHP limit using + // the WELTARG keyword + const int ok_bhp = + well_controls_add_new(BHP, unit::convert::from(1.0, unit::atm), + NULL, ctrl); + + if (ok_resv != 0 && ok_bhp != 0) { xw.currentControls()[*rp] = 0; well_controls_set_current(ctrl, 0); }