diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 2d505ffb8..6b47a99d6 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -268,7 +268,7 @@ namespace Opm { /// \param[in] ReservoirState /// \param[in] FIPNUM for active cells not global cells. /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas. - std::vector + std::vector > computeFluidInPlace(const ReservoirState& x, const std::vector& fipnum); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 8d890438a..4aa32b7e3 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -2243,7 +2243,7 @@ typedef Eigen::Array - std::vector + std::vector > BlackoilModelBase:: computeFluidInPlace(const ReservoirState& x, const std::vector& fipnum) @@ -2282,7 +2282,10 @@ typedef Eigen::Array values(dims, V::Zero(7)); + std::vector > values(dims); + for (int i=0; i < dims; ++i) { + values[i].resize(7, 0.0); + } const V hydrocarbon = saturation[Oil].value() + saturation[Gas].value(); V hcpv; @@ -2356,7 +2359,11 @@ typedef Eigen::Array + std::vector > computeFluidInPlace(const ReservoirState& x, const std::vector& fipnum) const { diff --git a/opm/autodiff/NonlinearSolver.hpp b/opm/autodiff/NonlinearSolver.hpp index bb6670d80..34d0db906 100644 --- a/opm/autodiff/NonlinearSolver.hpp +++ b/opm/autodiff/NonlinearSolver.hpp @@ -21,11 +21,9 @@ #ifndef OPM_NONLINEARSOLVER_HEADER_INCLUDED #define OPM_NONLINEARSOLVER_HEADER_INCLUDED -#include #include #include #include -#include #include #include #include @@ -39,11 +37,6 @@ namespace Opm { class NonlinearSolver { public: - // --------- Types and enums --------- - typedef AutoDiffBlock ADB; - typedef ADB::V V; - typedef ADB::M M; - // Available relaxation scheme types. enum RelaxType { DAMPEN, SOR }; @@ -132,14 +125,13 @@ namespace Opm { /// \param[in] ReservoirState /// \param[in] FIPNUM for active cells not global cells. /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas. - std::vector - computeFluidInPlace(const ReservoirState& x, - const std::vector& fipnum) const + std::vector > + computeFluidInPlace(const ReservoirState& x, const std::vector& fipnum) const { return model_->computeFluidInPlace(x, fipnum); } - std::vector> + std::vector > computeFluidInPlace(const std::vector& fipnum) const { return model_->computeFluidInPlace(fipnum); @@ -156,9 +148,6 @@ namespace Opm { void detectOscillations(const std::vector>& residual_history, const int it, bool& oscillate, bool& stagnate) const; - /// Apply a stabilization to dx, depending on dxOld and relaxation parameters. - void stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const; - /// Apply a stabilization to dx, depending on dxOld and relaxation parameters. /// Implemention for Dune block vectors. template diff --git a/opm/autodiff/NonlinearSolver_impl.hpp b/opm/autodiff/NonlinearSolver_impl.hpp index d4d300f2a..62ccf24c7 100644 --- a/opm/autodiff/NonlinearSolver_impl.hpp +++ b/opm/autodiff/NonlinearSolver_impl.hpp @@ -248,36 +248,6 @@ namespace Opm } - template - void - NonlinearSolver::stabilizeNonlinearUpdate(V& dx, V& dxOld, const double omega) const - { - // The dxOld is updated with dx. - // If omega is equal to 1., no relaxtion will be appiled. - - const V tempDxOld = dxOld; - dxOld = dx; - - switch (relaxType()) { - case DAMPEN: - if (omega == 1.) { - return; - } - dx = dx*omega; - return; - case SOR: - if (omega == 1.) { - return; - } - dx = dx*omega + (1.-omega)*tempDxOld; - return; - default: - OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type."); - } - - return; - } - template template void @@ -294,15 +264,19 @@ namespace Opm if (omega == 1.) { return; } - dx *= omega; + for (size_t i = 0; i < dx.size(); ++i) { + dx[i] *= omega; + } return; case SOR: if (omega == 1.) { return; } - dx *= omega; - tempDxOld *= (1.-omega); - dx += tempDxOld; + for (size_t i = 0; i < dx.size(); ++i) { + dx[i] *= omega; + tempDxOld[i] *= (1.-omega); + dx[i] += tempDxOld[i]; + } return; default: OPM_THROW(std::runtime_error, "Can only handle DAMPEN and SOR relaxation type."); @@ -310,8 +284,6 @@ namespace Opm return; } - - } // namespace Opm diff --git a/opm/autodiff/SimulatorBase.hpp b/opm/autodiff/SimulatorBase.hpp index df852de4d..f7322c88c 100644 --- a/opm/autodiff/SimulatorBase.hpp +++ b/opm/autodiff/SimulatorBase.hpp @@ -21,6 +21,9 @@ #ifndef OPM_SIMULATORBASE_HEADER_INCLUDED #define OPM_SIMULATORBASE_HEADER_INCLUDED +#include +#include + #include #include #include @@ -159,17 +162,18 @@ namespace Opm WellState& xw); void - FIPUnitConvert(const UnitSystem& units, V& fip); + FIPUnitConvert(const UnitSystem& units, + std::vector >& fip); void FIPUnitConvert(const UnitSystem& units, - std::vector& fip); - - V - FIPTotals(const std::vector& fip, const ReservoirState& state); + std::vector& fip); + + std::vector + FIPTotals(const std::vector >& fip, const ReservoirState& state); void - outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg); + outputFluidInPlace(const std::vector& oip, const std::vector& cip, const UnitSystem& units, const int reg); void computeWellPotentials(const Wells* wells, const WellState& xw, diff --git a/opm/autodiff/SimulatorBase_impl.hpp b/opm/autodiff/SimulatorBase_impl.hpp index f44f4726c..cc760a7bc 100644 --- a/opm/autodiff/SimulatorBase_impl.hpp +++ b/opm/autodiff/SimulatorBase_impl.hpp @@ -149,7 +149,7 @@ namespace Opm fipnum[c] = fipnum_global[AutoDiffGrid::globalCell(grid_)[c]]; } } - std::vector OOIP; + std::vector > OOIP; // Main simulation loop. while (!timer.done()) { // Report timestep. @@ -275,10 +275,10 @@ namespace Opm report.solver_time += solver_timer.secsSinceStart(); // Compute current FIP. - std::vector COIP; + std::vector > COIP; COIP = solver->computeFluidInPlace(state, fipnum); - V OOIP_totals = FIPTotals(OOIP, state); - V COIP_totals = FIPTotals(COIP, state); + std::vector OOIP_totals = FIPTotals(OOIP, state); + std::vector COIP_totals = FIPTotals(COIP, state); //Convert to correct units FIPUnitConvert(eclipse_state_->getUnits(), COIP); @@ -665,21 +665,19 @@ namespace Opm } } - template void SimulatorBase::FIPUnitConvert(const UnitSystem& units, - std::vector& fip) + std::vector >& fip) { for (size_t i = 0; i < fip.size(); ++i) { FIPUnitConvert(units, fip[i]); } } - template void - SimulatorBase::FIPUnitConvert(const UnitSystem& units, V& fip) + SimulatorBase::FIPUnitConvert(const UnitSystem& units, std::vector& fip) { if (units.getType() == UnitSystem::UnitType::UNIT_TYPE_FIELD) { fip[0] = unit::convert::to(fip[0], unit::stb); @@ -700,10 +698,10 @@ namespace Opm template - V - SimulatorBase::FIPTotals(const std::vector& fip, const ReservoirState& state) + std::vector + SimulatorBase::FIPTotals(const std::vector >& fip, const ReservoirState& state) { - V totals(V::Zero(7)); + std::vector totals(7, 0.0); for (int i = 0; i < 5; ++i) { for (size_t reg = 0; reg < fip.size(); ++reg) { totals[i] += fip[reg][i]; @@ -713,15 +711,34 @@ namespace Opm const int np = state.numPhases(); const PhaseUsage& pu = props_.phaseUsage(); const DataBlock s = Eigen::Map(& state.saturation()[0], nc, np); - const V so = pu.phase_used[BlackoilPhases::Liquid] ? V(s.col(BlackoilPhases::Liquid)) : V::Zero(nc); - const V sg = pu.phase_used[BlackoilPhases::Vapour] ? V(s.col(BlackoilPhases::Vapour)) : V::Zero(nc); - const V hydrocarbon = so + sg; - const V p = Eigen::Map(& state.pressure()[0], nc); + std::vector so(nc); + std::vector sg(nc); + std::vector hydrocarbon(nc); + const auto& soCol = s.col(BlackoilPhases::Liquid); + const auto& sgCol = s.col(BlackoilPhases::Vapour); + for (unsigned c = 0; c < so.size(); ++ c) { + double mySo = 0.0; + if (pu.phase_used[BlackoilPhases::Liquid]) + mySo = soCol[c]; + double mySg = 0.0; + if (pu.phase_used[BlackoilPhases::Vapour]) + mySg = sgCol[c]; + so[c] = mySo; + sg[c] = mySg; + hydrocarbon[c] = mySo + mySg; + } + const std::vector p = state.pressure(); if ( ! is_parallel_run_ ) { + double tmp = 0.0; + double tmp2 = 0.0; + for (unsigned i = 0; i < p.size(); ++i) { + tmp += p[i] * geo_.poreVolume()[i] * hydrocarbon[i]; + tmp2 += geo_.poreVolume()[i] * hydrocarbon[i]; + } totals[5] = geo_.poreVolume().sum(); - totals[6] = (p * geo_.poreVolume() * hydrocarbon).sum() / ((geo_.poreVolume() * hydrocarbon).sum()); - } + totals[6] = unit::convert::to(tmp/tmp2, unit::barsa); + } else { #if HAVE_MPI @@ -730,8 +747,12 @@ namespace Opm auto operators = std::make_tuple(Opm::Reduction::makeGlobalSumFunctor(), Opm::Reduction::makeGlobalSumFunctor(), Opm::Reduction::makeGlobalSumFunctor()); - auto pav_nom = p * geo_.poreVolume() * hydrocarbon; - auto pav_denom = geo_.poreVolume() * hydrocarbon; + std::vector pav_nom(p.size()); + std::vector pav_denom(pav_nom.size()); + for (unsigned i = 0; i < p.size(); ++i) { + pav_nom[i] = p[i] * geo_.poreVolume()[i] * hydrocarbon[i]; + pav_denom[i] = geo_.poreVolume()[i] * hydrocarbon[i]; + } // using ref cref to prevent copying auto inputs = std::make_tuple(std::cref(geo_.poreVolume()), @@ -756,7 +777,7 @@ namespace Opm template void - SimulatorBase::outputFluidInPlace(const V& oip, const V& cip, const UnitSystem& units, const int reg) + SimulatorBase::outputFluidInPlace(const std::vector& oip, const std::vector& cip, const UnitSystem& units, const int reg) { std::ostringstream ss; if (!reg) { diff --git a/opm/autodiff/VFPInjProperties.cpp b/opm/autodiff/VFPInjProperties.cpp index d029b495f..95ede87b4 100644 --- a/opm/autodiff/VFPInjProperties.cpp +++ b/opm/autodiff/VFPInjProperties.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include @@ -67,14 +68,10 @@ VFPInjProperties::VFPInjProperties(const std::map& tables) { } } - - - - VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, - const Wells& wells, - const ADB& qs, - const ADB& thp_val) const { + const Wells& wells, + const ADB& qs, + const ADB& thp_val) const { const int nw = wells.number_of_wells; //Short-hands for water / oil / gas phases @@ -87,14 +84,11 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, return bhp(table_id, w, o, g, thp_val); } - - - VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, - const ADB& aqua, - const ADB& liquid, - const ADB& vapour, - const ADB& thp_arg) const { + const ADB& aqua, + const ADB& liquid, + const ADB& vapour, + const ADB& thp_arg) const { const int nw = thp_arg.size(); std::vector block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg); diff --git a/opm/autodiff/VFPInjProperties.hpp b/opm/autodiff/VFPInjProperties.hpp index 1f69f006d..4987fdcd7 100644 --- a/opm/autodiff/VFPInjProperties.hpp +++ b/opm/autodiff/VFPInjProperties.hpp @@ -22,7 +22,6 @@ #include #include -#include #include #include #include @@ -34,6 +33,8 @@ namespace Opm { +template +class AutoDiffBlock; class VFPInjProperties { public: @@ -154,10 +155,10 @@ public: * the above parameters from the values in the input table. */ double bhp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& thp) const; + const double& aqua, + const double& liquid, + const double& vapour, + const double& thp) const; /** @@ -172,10 +173,10 @@ public: * the above parameters from the values in the input table. */ double thp(int table_id, - const double& aqua, - const double& liquid, - const double& vapour, - const double& bhp) const; + const double& aqua, + const double& liquid, + const double& vapour, + const double& bhp) const; /** * Returns the table associated with the ID, or throws an exception if diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index 2af02f812..e2a9063da 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include #include #include @@ -54,14 +55,11 @@ VFPProdProperties::VFPProdProperties(const std::map& tables) } } - - - VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, - const Wells& wells, - const ADB& qs, - const ADB& thp_arg, - const ADB& alq) const { + const Wells& wells, + const ADB& qs, + const ADB& thp_arg, + const ADB& alq) const { const int nw = wells.number_of_wells; //Short-hands for water / oil / gas phases @@ -76,11 +74,11 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, - const ADB& aqua, - const ADB& liquid, - const ADB& vapour, - const ADB& thp_arg, - const ADB& alq) const { + const ADB& aqua, + const ADB& liquid, + const ADB& vapour, + const ADB& thp_arg, + const ADB& alq) const { const int nw = thp_arg.size(); std::vector block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp_arg, alq); diff --git a/opm/autodiff/VFPProdProperties.hpp b/opm/autodiff/VFPProdProperties.hpp index 31b2238c8..7d278e420 100644 --- a/opm/autodiff/VFPProdProperties.hpp +++ b/opm/autodiff/VFPProdProperties.hpp @@ -22,7 +22,6 @@ #include #include -#include #include #include #include @@ -33,6 +32,8 @@ namespace Opm { +template +class AutoDiffBlock; /** * Class which linearly interpolates BHP as a function of rate, tubing head pressure, diff --git a/opm/autodiff/WellHelpers.hpp b/opm/autodiff/WellHelpers.hpp index 6e0728c88..c26168ef3 100644 --- a/opm/autodiff/WellHelpers.hpp +++ b/opm/autodiff/WellHelpers.hpp @@ -23,7 +23,6 @@ #define OPM_WELLHELPERS_HEADER_INCLUDED #include -#include // #include #include @@ -121,7 +120,6 @@ namespace Opm { // --------- Types --------- - using Vector = AutoDiffBlock::V; /** * Simple hydrostatic correction for VFP table @@ -149,7 +147,7 @@ namespace Opm { return dp; } - + template inline Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth, const Vector& well_perforation_densities, const double gravity) { diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 4badc6b00..d376885a5 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -536,7 +536,7 @@ namespace { - std::vector + std::vector > FullyImplicitCompressiblePolymerSolver::computeFluidInPlace(const PolymerBlackoilState& x, const std::vector& fipnum) { @@ -568,7 +568,11 @@ namespace { const int dims = *std::max_element(fipnum.begin(), fipnum.end()); - std::vector values(dims, V::Zero(7)); + std::vector > values(dims); + for (int i=0; i < dims; ++i) { + values[i].resize(7, 0.0); + } + V hcpv = V::Zero(nc); V pres = V::Zero(nc); for (int i = 0; i < 5; ++i) { diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp index 096e8202a..beb0e9b0c 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp @@ -158,7 +158,7 @@ namespace Opm { /// \param[in] WellState /// \param[in] FIPNUM for active cells not global cells. /// \return fluid in place, number of fip regions, each region contains 5 values which are liquid, vapour, water, free gas and dissolved gas. - std::vector + std::vector > computeFluidInPlace(const PolymerBlackoilState& x, const std::vector& fipnum);