From 948d985f56a87f57e811e81e4280d9a4652340dc Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 4 Mar 2016 11:15:46 +0100 Subject: [PATCH 01/24] Add support for PMISC Pressure effects are added to relative permeability, capillary pressure, viscosity and density miscibility --- opm/autodiff/BlackoilSolventModel.hpp | 10 +- opm/autodiff/BlackoilSolventModel_impl.hpp | 122 ++++++++++++++++++--- opm/autodiff/SolventPropsAdFromDeck.cpp | 29 +++++ opm/autodiff/SolventPropsAdFromDeck.hpp | 8 ++ 4 files changed, 152 insertions(+), 17 deletions(-) diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index e8cc9a063..bc2a72ebf 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -144,7 +144,6 @@ namespace Opm { using Base::wellsActive; using Base::wells; using Base::variableState; - using Base::computePressures; using Base::computeGasPressure; using Base::applyThresholdPressures; using Base::fluidRsSat; @@ -157,7 +156,6 @@ namespace Opm { using Base::dsMax; using Base::drMaxRel; using Base::maxResidualAllowed; - using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; @@ -237,6 +235,14 @@ namespace Opm { // compute density and viscosity using the ToddLongstaff mixing model void computeToddLongstaffMixing(std::vector& viscosity, std::vector& density, const std::vector& saturations, const Opm::PhaseUsage pu); + // compute phase pressures. + std::vector + computePressures(const ADB& po, + const ADB& sw, + const ADB& so, + const ADB& sg, + const ADB& ss) const; + }; diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index f7fda17b0..3f703a045 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -177,17 +177,72 @@ namespace Opm { BlackoilSolventModel::variableStateExtractVars(const ReservoirState& x, const std::vector& indices, std::vector& vars) const - { - SolutionState state = Base::variableStateExtractVars(x, indices, vars); - if (has_solvent_) { - state.solvent_saturation = std::move(vars[indices[Solvent]]); + { + // This is more or less a copy of the base class. Refactoring is needed in the base class + // to avoid unnecessary copying. + + //using namespace Opm::AutoDiffGrid; + const int nc = Opm::AutoDiffGrid::numCells(grid_); + const Opm::PhaseUsage pu = fluid_.phaseUsage(); + + SolutionState state(fluid_.numPhases()); + + // Pressure. + state.pressure = std::move(vars[indices[Pressure]]); + + // Temperature cannot be a variable at this time (only constant). + const V temp = Eigen::Map(& x.temperature()[0], x.temperature().size()); + state.temperature = ADB::constant(temp); + + // Saturations + { + ADB so = ADB::constant(V::Ones(nc, 1)); + + if (active_[ Water ]) { + state.saturation[pu.phase_pos[ Water ]] = std::move(vars[indices[Sw]]); + const ADB& sw = state.saturation[pu.phase_pos[ Water ]]; + so -= sw; + } + if (has_solvent_) { + state.solvent_saturation = std::move(vars[indices[Solvent]]); + so -= state.solvent_saturation; + } + + 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[indices[Xvar]]; + ADB& sg = state.saturation[ pu.phase_pos[ Gas ] ]; + sg = Base::isSg_*xvar + Base::isRv_*so; + so -= sg; + + if (active_[ Oil ]) { + // RS and RV is only defined if both oil and gas phase are active. + state.canonical_phase_pressures = computePressures(state.pressure, state.saturation[pu.phase_pos[ Water ]], so, sg, state.solvent_saturation); + const ADB rsSat = fluidRsSat(state.canonical_phase_pressures[ Oil ], so , cells_); + if (has_disgas_) { + state.rs = (1-Base::isRs_)*rsSat + Base::isRs_*xvar; + } else { + state.rs = rsSat; + } + const ADB rvSat = fluidRvSat(state.canonical_phase_pressures[ Gas ], so , cells_); + if (has_vapoil_) { + state.rv = (1-Base::isRv_)*rvSat + Base::isRv_*xvar; + } else { + state.rv = rvSat; + } + } + } + if (active_[ Oil ]) { // Note that so is never a primary variable. - const Opm::PhaseUsage pu = fluid_.phaseUsage(); - state.saturation[pu.phase_pos[ Oil ]] -= state.solvent_saturation; + state.saturation[pu.phase_pos[ Oil ]] = std::move(so); } } + // wells + Base::variableStateExtractWellsVars(indices, vars, state); return state; + } @@ -665,7 +720,9 @@ namespace Opm { Selector zero_selector(ss.value() + sg.value(), Selector::Zero); ADB F_solvent = zero_selector.select(ss, ss / (ss + sg)); - const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_); + const ADB& po = state.canonical_phase_pressures[ Oil ]; + const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_) + * solvent_props_.pressureMiscibilityFunction(po, cells_); assert(active_[ Oil ]); assert(active_[ Gas ]); @@ -769,16 +826,27 @@ namespace Opm { // Compute effective viscosities and densities computeToddLongstaffMixing(viscosity, density, effective_saturations, pu); - // Store the computed volume factors and viscosities - b_eff_[pu.phase_pos[ Water ]] = bw; - b_eff_[pu.phase_pos[ Oil ]] = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs); - b_eff_[pu.phase_pos[ Gas ]] = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv); - b_eff_[solvent_pos_] = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_); + // compute the volume factors from the densities + const ADB b_eff_o = density[pu.phase_pos[ Oil ]] / (fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) * state.rs); + const ADB b_eff_g = density[pu.phase_pos[ Gas ]] / (fluid_.surfaceDensity(pu.phase_pos[ Gas ], cells_) + fluid_.surfaceDensity(pu.phase_pos[ Oil ], cells_) * state.rv); + const ADB b_eff_s = density[solvent_pos_] / solvent_props_.solventSurfaceDensity(cells_); + // account for pressure effects and store the computed volume factors and viscosities + const V ones = V::Constant(nc, 1.0); + const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_); + + b_eff_[pu.phase_pos[ Oil ]] = pmisc * b_eff_o + (ones - pmisc) * bo; + b_eff_[pu.phase_pos[ Gas ]] = pmisc * b_eff_g + (ones - pmisc) * bg; + b_eff_[solvent_pos_] = pmisc * b_eff_s + (ones - pmisc) * bs; + + // keep the mu*b interpolation + mu_eff_[pu.phase_pos[ Oil ]] = b_eff_[pu.phase_pos[ Oil ]] / (pmisc * b_eff_o / viscosity[pu.phase_pos[ Oil ]] + (ones - pmisc) * bo / mu_o); + mu_eff_[pu.phase_pos[ Gas ]] = b_eff_[pu.phase_pos[ Gas ]] / (pmisc * b_eff_g / viscosity[pu.phase_pos[ Gas ]] + (ones - pmisc) * bg / mu_g); + mu_eff_[solvent_pos_] = b_eff_[solvent_pos_] / (pmisc * b_eff_s / viscosity[solvent_pos_] + (ones - pmisc) * bs / mu_s); + + // for water the pure values are used mu_eff_[pu.phase_pos[ Water ]] = mu_w; - mu_eff_[pu.phase_pos[ Oil ]] = viscosity[pu.phase_pos[ Oil ]]; - mu_eff_[pu.phase_pos[ Gas ]] = viscosity[pu.phase_pos[ Gas ]]; - mu_eff_[solvent_pos_] = viscosity[solvent_pos_]; + b_eff_[pu.phase_pos[ Water ]] = bw; } template @@ -881,6 +949,30 @@ namespace Opm { } + template + std::vector + BlackoilSolventModel:: + computePressures(const ADB& po, + const ADB& sw, + const ADB& so, + const ADB& sg, + const ADB& ss) const + { + std::vector pressures = Base::computePressures(po, sw, so, sg); + // The imiscible capillary pressure is evaluated using the total gas saturation (sg + ss) + std::vector pressures_imisc = Base::computePressures(po, sw, so, sg + ss); + + // Pressure effects on capillary pressure miscibility + const ADB pmisc = solvent_props_.pressureMiscibilityFunction(po, cells_); + // Only the pcog is effected by the miscibility. Since pg = po + pcog, changing pg is eqvivalent + // to changing the gas pressure directly. + const int nc = cells_.size(); + const V ones = V::Constant(nc, 1.0); + pressures[Gas] = ( pmisc * pressures[Gas] + ((ones - pmisc) * pressures_imisc[Gas])); + return pressures; + } + + template void BlackoilSolventModel::assemble(const ReservoirState& reservoir_state, diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index 2fd967021..b94d1a64c 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/opm/autodiff/SolventPropsAdFromDeck.cpp @@ -172,6 +172,25 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n"); } + const TableContainer& pmiscTables = tables->getPmiscTables(); + if (!pmiscTables.empty()) { + + int numRegions = pmiscTables.size(); + + // resize the attributes of the object + pmisc_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::PmiscTable& pmiscTable = pmiscTables.getTable(regionIdx); + + // Copy data + const auto& po = pmiscTable.getOilPhasePressureColumn(); + const auto& pmisc = pmiscTable.getMiscibilityColumn(); + + pmisc_[regionIdx] = NonuniformTableLinear(po, pmisc); + + } + } + // miscible relative permeability multipleiers const TableContainer& msfnTables = tables->getMsfnTables(); if (!msfnTables.empty()) { @@ -339,6 +358,16 @@ ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction, return SolventPropsAdFromDeck::makeADBfromTables(solventFraction, cells, cellMiscRegionIdx_, misc_); } +ADB SolventPropsAdFromDeck::pressureMiscibilityFunction(const ADB& po, + const Cells& cells) const +{ + if (pmisc_.size() > 0) { + return SolventPropsAdFromDeck::makeADBfromTables(po, cells, cellMiscRegionIdx_, pmisc_); + } + // return ones if not specified i.e. no effect. + return ADB::constant(V::Constant(po.size(), 1.0)); +} + ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw, const Cells& cells) const { diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index 4aa69b4c4..cdc94022e 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.hpp +++ b/opm/autodiff/SolventPropsAdFromDeck.hpp @@ -108,6 +108,13 @@ public: ADB miscibilityFunction(const ADB& solventFraction, const Cells& cells) const; + /// Pressure dependent miscibility function + /// \param[in] solventFraction Array of n oil phase pressure . + /// \param[in] cells Array of n cell indices to be associated with the pressure values. + /// \return Array of n miscibility values + ADB pressureMiscibilityFunction(const ADB& po, + const Cells& cells) const; + /// Miscible critical gas saturation function /// \param[in] Sw Array of n water saturation values. /// \param[in] cells Array of n cell indices to be associated with the saturation values. @@ -177,6 +184,7 @@ private: std::vector > mkro_; std::vector > mkrsg_; std::vector > misc_; + std::vector > pmisc_; std::vector > sorwmis_; std::vector > sgcwmis_; std::vector mix_param_viscosity_; From 20819b8a037b6175ecab26649c4f70d5a8f0631f Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 7 Mar 2016 09:19:42 +0100 Subject: [PATCH 02/24] Skip the relperm diagnostics for the solvent model. --- opm/autodiff/FlowMainSolvent.hpp | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/opm/autodiff/FlowMainSolvent.hpp b/opm/autodiff/FlowMainSolvent.hpp index d3307f31e..5e1346996 100644 --- a/opm/autodiff/FlowMainSolvent.hpp +++ b/opm/autodiff/FlowMainSolvent.hpp @@ -107,6 +107,17 @@ namespace Opm Base::deck_->hasKeyword("SOLVENT"))); } + // run diagnostics + // Writes to: + // logFile_ + void runDiagnostics() + { + std::shared_ptr streamLog = std::make_shared(Base::logFile_, Opm::Log::DefaultMessageTypes); + const std::string msg = "Relperm diagnostic not implemented for the solvent model "; + std::cout << "\n" << msg << "\n" << std::endl; + streamLog->addMessage(Log::MessageType::Info, msg); + } + }; From b39e579197b248fcf801ecfd9edca39fa7bae9a3 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 17 Feb 2016 12:24:13 +0100 Subject: [PATCH 03/24] fix the fallout of the removal of the opm-core black-oil PVT classes --- opm/autodiff/BlackoilPropsAdFromDeck.cpp | 23 ++++++++--------------- opm/autodiff/SolventPropsAdFromDeck.cpp | 4 ++-- opm/autodiff/SolventPropsAdFromDeck.hpp | 4 +--- 3 files changed, 11 insertions(+), 20 deletions(-) diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index 904c4faad..97468642c 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -22,27 +22,20 @@ #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 Opm { // Making these typedef to make the code more readable. diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index 2fd967021..880c7fad6 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/opm/autodiff/SolventPropsAdFromDeck.cpp @@ -21,12 +21,13 @@ #include #include +#include + #include #include #include #include - namespace Opm { // Making these typedef to make the code more readable. @@ -383,7 +384,6 @@ ADB SolventPropsAdFromDeck::makeADBfromTables(const ADB& X_AD, } - V SolventPropsAdFromDeck::solventSurfaceDensity(const Cells& cells) const { const int n = cells.size(); V density(n); diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index 4aa69b4c4..e22f01847 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.hpp +++ b/opm/autodiff/SolventPropsAdFromDeck.hpp @@ -23,8 +23,7 @@ #include #include -#include -#include +#include #include #include @@ -139,7 +138,6 @@ public: private: - /// Makes ADB from table values /// \param[in] X Array of n table lookup values. /// \param[in] cells Array of n cell indices to be associated with the fraction values. From b2e02f6d2b465e734bec7300e77277cbc4b1bede Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Mon, 7 Mar 2016 14:23:45 +0100 Subject: [PATCH 04/24] Add power method for general f^g in the AutoDiffBlock A power method where both f and g are ADB variables is added using the general derivative rule (f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f' Tests are added to test_block.cpp --- opm/autodiff/AutoDiffBlock.hpp | 45 ++++++++++++++++++++++++++++++++++ tests/test_block.cpp | 35 +++++++++++++++++++++++++- 2 files changed, 79 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index 66901ea4e..b0c2a4d8a 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -1,5 +1,6 @@ /* Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2016 IRIS AS This file is part of the Open Porous Media project (OPM). @@ -676,6 +677,50 @@ namespace Opm return AutoDiffBlock::function( std::move(val), std::move(jac) ); } + /** + * @brief Computes the value of base raised to the power of exp + * + * @param base The base AD forward block + * @param exp The exponent AD forward block + * @return The value of base raised to the power of exp + */ template + AutoDiffBlock pow(const AutoDiffBlock& base, + const AutoDiffBlock& exp) + { + const int num_elem = base.value().size(); + assert(exp.value().size() == num_elem); + typename AutoDiffBlock::V val (num_elem); + for (int i = 0; i < num_elem; ++i) { + val[i] = std::pow(base.value()[i], exp.value()[i]); + } + + // (f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f' + typename AutoDiffBlock::V der1 = val; + for (int i = 0; i < num_elem; ++i) { + der1[i] *= std::log(base.value()[i]); + } + std::vector< typename AutoDiffBlock::M > jac1 (base.numBlocks()); + const typename AutoDiffBlock::M der1_diag(der1.matrix().asDiagonal()); + for (int block = 0; block < base.numBlocks(); block++) { + fastSparseProduct(der1_diag, exp.derivative()[block], jac1[block]); + } + typename AutoDiffBlock::V der2 = exp.value(); + for (int i = 0; i < num_elem; ++i) { + der2[i] *= std::pow(base.value()[i], exp.value()[i] - 1.0); + } + std::vector< typename AutoDiffBlock::M > jac2 (base.numBlocks()); + const typename AutoDiffBlock::M der2_diag(der2.matrix().asDiagonal()); + for (int block = 0; block < base.numBlocks(); block++) { + fastSparseProduct(der2_diag, base.derivative()[block], jac2[block]); + } + std::vector< typename AutoDiffBlock::M > jac (base.numBlocks()); + for (int block = 0; block < base.numBlocks(); block++) { + jac[block] = jac1[block] + jac2[block]; + } + + return AutoDiffBlock::function(std::move(val), std::move(jac)); + } + } // namespace Opm diff --git a/tests/test_block.cpp b/tests/test_block.cpp index 30928d400..48cbb3fce 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -1,5 +1,6 @@ /* Copyright 2013 SINTEF ICT, Applied Mathematics. + Copyright 2016 IRIS AS This file is part of the Open Porous Media project (OPM). @@ -293,7 +294,7 @@ BOOST_AUTO_TEST_CASE(Pow) vx << 0.2, 1.2, 13.4; ADB::V vy(3); - vy << 1.0, 2.2, 3.4; + vy << 2.0, 3.0, 0.5; std::vector vals{ vx, vy }; std::vector vars = ADB::variables(vals); @@ -303,6 +304,7 @@ BOOST_AUTO_TEST_CASE(Pow) const double tolerance = 1e-14; + // test exp = double const ADB xx = x * x; ADB xxpow2 = Opm::pow(x,2.0); checkClose(xxpow2, xx, tolerance); @@ -322,6 +324,37 @@ BOOST_AUTO_TEST_CASE(Pow) for (int i = 0 ; i < 3; ++i){ BOOST_CHECK_CLOSE(xpowhalf.value()[i], x_sqrt[i], 1e-4); } + + // test exp = ADB::V + ADB xpowyval = Opm::pow(x,y.value()); + + // each of the component of y is tested in the test above + // we compare with the results from the above tests. + ADB::V pick1(3); + pick1 << 1,0,0; + ADB::V pick2(3); + pick2 << 0,1,0; + ADB::V pick3(3); + pick3 << 0,0,1; + + ADB compare = pick1 * xx + pick2 * xxx + pick3 * xpowhalf; + checkClose(xpowyval, compare, tolerance); + + // test exp = ADB + ADB xpowy = Opm::pow(x,y); + + // the value and the first jacobian should be equal to the xpowyval + // the second jacobian is hand calculated + // log(0.2)*0.2^2.0, log(1.2) * 1.2^3.0, log(13.4) * 13.4^0.5 + ADB::V jac2(3); + jac2 << -0.0643775165 , 0.315051650 , 9.50019208855; + for (int i = 0 ; i < 3; ++i){ + BOOST_CHECK_CLOSE(xpowy.value()[i], xpowyval.value()[i], tolerance); + BOOST_CHECK_CLOSE(xpowy.derivative()[0].coeff(i,i), xpowyval.derivative()[0].coeff(i,i), tolerance); + BOOST_CHECK_CLOSE(xpowy.derivative()[1].coeff(i,i), jac2[i], 1e-4); + } + + } From 35b34f1b3af5e7ed28e5123b9071aaea635e30af Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Tue, 8 Mar 2016 10:04:19 +0100 Subject: [PATCH 05/24] Add pow() for constant base raised to variable exponent in AutoDiffBlock - associated tests are added - this PR also contains some clean up --- opm/autodiff/AutoDiffBlock.hpp | 144 ++++++++++++++++++--------------- tests/test_block.cpp | 23 ++++-- 2 files changed, 98 insertions(+), 69 deletions(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index b0c2a4d8a..71c42b746 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -624,49 +624,20 @@ namespace Opm return rhs * lhs; // Commutative operation. } + /** - * @brief Computes the value of base raised to the power of exp elementwise + * @brief Computes the value of base raised to the power of exponent * * @param base The AD forward block - * @param exp array of exponents - * @return The value of base raised to the power of exp elementwise + * @param exponent double + * @return The value of base raised to the power of exponent */ template AutoDiffBlock pow(const AutoDiffBlock& base, - const typename AutoDiffBlock::V& exp) + const double exponent) { - const int num_elem = base.value().size(); - typename AutoDiffBlock::V val (num_elem); - typename AutoDiffBlock::V derivative = exp; - assert(exp.size() == num_elem); - for (int i = 0; i < num_elem; ++i) { - val[i] = std::pow(base.value()[i], exp[i]); - derivative[i] *= std::pow(base.value()[i], exp[i] - 1.0); - } - const typename AutoDiffBlock::M derivative_diag(derivative.matrix().asDiagonal()); - - std::vector< typename AutoDiffBlock::M > jac (base.numBlocks()); - for (int block = 0; block < base.numBlocks(); block++) { - fastSparseProduct(derivative_diag, base.derivative()[block], jac[block]); - } - - return AutoDiffBlock::function( std::move(val), std::move(jac) ); - } - - - /** - * @brief Computes the value of base raised to the power of exp - * - * @param base The AD forward block - * @param exp exponent - * @return The value of base raised to the power of exp - */ - template - AutoDiffBlock pow(const AutoDiffBlock& base, - const double exp) - { - const typename AutoDiffBlock::V val = base.value().pow(exp); - const typename AutoDiffBlock::V derivative = exp * base.value().pow(exp - 1.0); + const typename AutoDiffBlock::V val = base.value().pow(exponent); + const typename AutoDiffBlock::V derivative = exponent * base.value().pow(exponent - 1.0); const typename AutoDiffBlock::M derivative_diag(derivative.matrix().asDiagonal()); std::vector< typename AutoDiffBlock::M > jac (base.numBlocks()); @@ -678,44 +649,91 @@ namespace Opm } /** - * @brief Computes the value of base raised to the power of exp + * @brief Computes the value of base raised to the power of exponent + * + * @param base The AD forward block + * @param exponent Array of exponents + * @return The value of base raised to the power of exponent elementwise + */ + template + AutoDiffBlock pow(const AutoDiffBlock& base, + const typename AutoDiffBlock::V& exponent) + { + // Add trivial derivatives and use the AD pow function + return pow(base,AutoDiffBlock::constant(exponent)); + } + + /** + * @brief Computes the value of base raised to the power of exponent + * + * @param base Array of base values + * @param exponent The AD forward block + * @return The value of base raised to the power of exponent elementwise + */ + template + AutoDiffBlock pow(const typename AutoDiffBlock::V& base, + const AutoDiffBlock& exponent) + { + // Add trivial derivatives and use the AD pow function + return pow(AutoDiffBlock::constant(base),exponent); + } + + /** + * @brief Computes the value of base raised to the power of exponent * * @param base The base AD forward block - * @param exp The exponent AD forward block + * @param exponent The exponent AD forward block * @return The value of base raised to the power of exp - */ template + */ + template AutoDiffBlock pow(const AutoDiffBlock& base, - const AutoDiffBlock& exp) - { + const AutoDiffBlock& exponent) + { const int num_elem = base.value().size(); - assert(exp.value().size() == num_elem); + assert(exponent.size() == num_elem); typename AutoDiffBlock::V val (num_elem); for (int i = 0; i < num_elem; ++i) { - val[i] = std::pow(base.value()[i], exp.value()[i]); + val[i] = std::pow(base.value()[i], exponent.value()[i]); } - // (f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f' - typename AutoDiffBlock::V der1 = val; - for (int i = 0; i < num_elem; ++i) { - der1[i] *= std::log(base.value()[i]); + // (f^g)' = f^g * ln(f) * g' + g * f^(g-1) * f' = der1 + der2 + // if f' is empty only der1 is calculated + // if g' is empty only der2 is calculated + // if f' and g' are non empty they should have the same size + int num_blocks = std::max (base.numBlocks(), exponent.numBlocks()); + if (!base.derivative().empty() && !exponent.derivative().empty()) { + assert(exponent.numBlocks() == base.numBlocks()); } - std::vector< typename AutoDiffBlock::M > jac1 (base.numBlocks()); - const typename AutoDiffBlock::M der1_diag(der1.matrix().asDiagonal()); - for (int block = 0; block < base.numBlocks(); block++) { - fastSparseProduct(der1_diag, exp.derivative()[block], jac1[block]); + std::vector< typename AutoDiffBlock::M > jac (num_blocks); + + if ( !exponent.derivative().empty() ) { + typename AutoDiffBlock::V der1 = val; + for (int i = 0; i < num_elem; ++i) { + der1[i] *= std::log(base.value()[i]); + } + std::vector< typename AutoDiffBlock::M > jac1 (exponent.numBlocks()); + const typename AutoDiffBlock::M der1_diag(der1.matrix().asDiagonal()); + for (int block = 0; block < exponent.numBlocks(); block++) { + fastSparseProduct(der1_diag, exponent.derivative()[block], jac1[block]); + jac[block] = jac1[block]; + } } - typename AutoDiffBlock::V der2 = exp.value(); - for (int i = 0; i < num_elem; ++i) { - der2[i] *= std::pow(base.value()[i], exp.value()[i] - 1.0); - } - std::vector< typename AutoDiffBlock::M > jac2 (base.numBlocks()); - const typename AutoDiffBlock::M der2_diag(der2.matrix().asDiagonal()); - for (int block = 0; block < base.numBlocks(); block++) { - fastSparseProduct(der2_diag, base.derivative()[block], jac2[block]); - } - std::vector< typename AutoDiffBlock::M > jac (base.numBlocks()); - for (int block = 0; block < base.numBlocks(); block++) { - jac[block] = jac1[block] + jac2[block]; + + if ( !base.derivative().empty() ) { + typename AutoDiffBlock::V der2 = exponent.value(); + for (int i = 0; i < num_elem; ++i) { + der2[i] *= std::pow(base.value()[i], exponent.value()[i] - 1.0); + } + std::vector< typename AutoDiffBlock::M > jac2 (base.numBlocks()); + const typename AutoDiffBlock::M der2_diag(der2.matrix().asDiagonal()); + for (int block = 0; block < base.numBlocks(); block++) { + fastSparseProduct(der2_diag, base.derivative()[block], jac2[block]); + if (!exponent.derivative().empty()) { + jac[block] += jac2[block]; + } else { + jac[block] = jac2[block]; + } + } } return AutoDiffBlock::function(std::move(val), std::move(jac)); diff --git a/tests/test_block.cpp b/tests/test_block.cpp index 48cbb3fce..ad3eda5fc 100644 --- a/tests/test_block.cpp +++ b/tests/test_block.cpp @@ -340,20 +340,31 @@ BOOST_AUTO_TEST_CASE(Pow) ADB compare = pick1 * xx + pick2 * xxx + pick3 * xpowhalf; checkClose(xpowyval, compare, tolerance); - // test exp = ADB - ADB xpowy = Opm::pow(x,y); + // test exponent = ADB::V and base = ADB + ADB xvalpowy = Opm::pow(x.value(),y); - // the value and the first jacobian should be equal to the xpowyval + // the value should be equal to xpowyval + // the first jacobian should be trivial // the second jacobian is hand calculated // log(0.2)*0.2^2.0, log(1.2) * 1.2^3.0, log(13.4) * 13.4^0.5 ADB::V jac2(3); jac2 << -0.0643775165 , 0.315051650 , 9.50019208855; for (int i = 0 ; i < 3; ++i){ - BOOST_CHECK_CLOSE(xpowy.value()[i], xpowyval.value()[i], tolerance); - BOOST_CHECK_CLOSE(xpowy.derivative()[0].coeff(i,i), xpowyval.derivative()[0].coeff(i,i), tolerance); - BOOST_CHECK_CLOSE(xpowy.derivative()[1].coeff(i,i), jac2[i], 1e-4); + BOOST_CHECK_CLOSE(xvalpowy.value()[i], xpowyval.value()[i], tolerance); + BOOST_CHECK_CLOSE(xvalpowy.derivative()[0].coeff(i,i), 0.0, tolerance); + BOOST_CHECK_CLOSE(xvalpowy.derivative()[1].coeff(i,i), jac2[i], 1e-4); } + // test exp = ADB + ADB xpowy = Opm::pow(x,y); + + // the first jacobian should be equal to the xpowyval + // the second jacobian should be equal to the xvalpowy + for (int i = 0 ; i < 3; ++i){ + BOOST_CHECK_CLOSE(xpowy.value()[i], xpowyval.value()[i], tolerance); + BOOST_CHECK_CLOSE(xpowy.derivative()[0].coeff(i,i), xpowyval.derivative()[0].coeff(i,i), tolerance); + BOOST_CHECK_CLOSE(xpowy.derivative()[1].coeff(i,i), xvalpowy.derivative()[1].coeff(i,i), tolerance); + } } From 8115d918dd728dbf99b15756b11efcc4e2fb604a Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 8 Mar 2016 15:44:53 +0100 Subject: [PATCH 06/24] adding flag to wops_ms_ to indicate if MSW is involved. to reduce some extra cost when MSW is not involved. --- opm/autodiff/BlackoilMultiSegmentModel.hpp | 1 + .../BlackoilMultiSegmentModel_impl.hpp | 91 ++++++++++++++----- 2 files changed, 70 insertions(+), 22 deletions(-) diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index 535537e91..e53647156 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -205,6 +205,7 @@ namespace Opm { AutoDiffMatrix eliminate_topseg; // change the top segment related to be zero std::vector well_cells; // the set of perforated cells V conn_trans_factors; // connection transmissibility factors + bool has_multisegment_wells; // flag indicating whether there is any muli-segment well }; MultiSegmentWellOps wops_ms_; diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index b09a64fe1..f58ac64c2 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -109,6 +109,9 @@ namespace Opm { BlackoilMultiSegmentModel:: MultiSegmentWellOps::MultiSegmentWellOps(const std::vector& wells_ms) { + // no multi-segment wells are involved by default. + has_multisegment_wells = false; + if (wells_ms.empty()) { return; } @@ -120,6 +123,9 @@ namespace Opm { for (int w = 0; w < nw; ++w) { total_nperf += wells_ms[w]->numberOfPerforations(); total_nseg += wells_ms[w]->numberOfSegments(); + if (wells_ms[w]->isMultiSegmented()) { + has_multisegment_wells = true; + } } // Create well_cells and conn_trans_factors. @@ -230,6 +236,12 @@ namespace Opm { top_well_segments_ = well_state.topSegmentLoc(); const int nw = wellsMultiSegment().size(); + + if ( !wops_ms_.has_multisegment_wells ) { + segvdt_ = V::Zero(nw); + return; + } + const int nseg_total = well_state.numSegments(); std::vector segment_volume; segment_volume.reserve(nseg_total); @@ -462,6 +474,12 @@ namespace Opm { well_perforation_densities_ = Eigen::Map(cd.data(), nperf_total); // This one is not useful for segmented wells at all well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf_total); + if ( !wops_ms_.has_multisegment_wells ) { + well_perforation_cell_densities_ = V::Zero(nperf_total); + well_perforation_cell_pressure_diffs_ = V::Zero(nperf_total); + return; + } + // compute the average of the fluid densites in the well blocks. // the average is weighted according to the fluid relative permeabilities. const std::vector kr_adb = Base::computeRelPerm(state); @@ -680,7 +698,7 @@ namespace Opm { // Compute drawdown. ADB h_nc = msperf_selector.select(well_segment_perforation_pressure_diffs_, - ADB::constant(well_perforation_pressure_diffs_)); + ADB::constant(well_perforation_pressure_diffs_)); const V h_cj = msperf_selector.select(well_perforation_cell_pressure_diffs_, V::Zero(nperf)); // Special handling for when we are called from solveWellEq(). @@ -890,23 +908,27 @@ namespace Opm { std::vector segment_volume_change_dt(np, ADB::null()); for (int phase = 0; phase < np; ++phase) { - // Gain of the surface volume of each component in the segment by dt - segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] - - segment_comp_surf_volume_initial_[phase]; + if ( wops_ms_.has_multisegment_wells ) { + // Gain of the surface volume of each component in the segment by dt + segment_volume_change_dt[phase] = segment_comp_surf_volume_current_[phase] - + segment_comp_surf_volume_initial_[phase]; - // Special handling for when we are called from solveWellEq(). - // TODO: restructure to eliminate need for special treatmemt. - if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) { - assert(segment_volume_change_dt[phase].numBlocks() > 2); - assert(segqs.numBlocks() == 2); - segment_volume_change_dt[phase] = detail::onlyWellDerivs(segment_volume_change_dt[phase]); - assert(segment_volume_change_dt[phase].numBlocks() == 2); + // Special handling for when we are called from solveWellEq(). + // TODO: restructure to eliminate need for special treatmemt. + if (segment_volume_change_dt[phase].numBlocks() != segqs.numBlocks()) { + assert(segment_volume_change_dt[phase].numBlocks() > 2); + assert(segqs.numBlocks() == 2); + segment_volume_change_dt[phase] = detail::onlyWellDerivs(segment_volume_change_dt[phase]); + assert(segment_volume_change_dt[phase].numBlocks() == 2); + } + + const ADB cq_s_seg = wops_ms_.p2s * cq_s[phase]; + const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total)); + segqs -= superset(cq_s_seg + wops_ms_.s2s_inlets * segqs_phase + segment_volume_change_dt[phase], + Span(nseg_total, 1, phase * nseg_total), np * nseg_total); + } else { + segqs -= superset(wops_ms_.p2s * cq_s[phase], Span(nseg_total, 1, phase * nseg_total), np * nseg_total); } - - const ADB cq_s_seg = wops_ms_.p2s * cq_s[phase]; - const ADB segqs_phase = subset(segqs, Span(nseg_total, 1, phase * nseg_total)); - segqs -= superset(cq_s_seg + wops_ms_.s2s_inlets * segqs_phase + segment_volume_change_dt[phase], - Span(nseg_total, 1, phase * nseg_total), np * nseg_total); } residual_.well_flux_eq = segqs; @@ -1169,13 +1191,17 @@ namespace Opm { ADB others_residual = ADB::constant(V::Zero(nseg_total)); - // Special handling for when we are called from solveWellEq(). - // TODO: restructure to eliminate need for special treatmemt. - ADB wspd = (state.segp.numBlocks() == 2) - ? detail::onlyWellDerivs(well_segment_pressures_delta_) - : well_segment_pressures_delta_; + if ( wops_ms_.has_multisegment_wells ) { + // Special handling for when we are called from solveWellEq(). + // TODO: restructure to eliminate need for special treatmemt. + ADB wspd = (state.segp.numBlocks() == 2) + ? detail::onlyWellDerivs(well_segment_pressures_delta_) + : well_segment_pressures_delta_; - others_residual = wops_ms_.eliminate_topseg * (state.segp - wops_ms_.s2s_outlet * state.segp + wspd); + others_residual = wops_ms_.eliminate_topseg * (state.segp - wops_ms_.s2s_outlet * state.segp + wspd); + } else { + others_residual = wops_ms_.eliminate_topseg * (state.segp - wops_ms_.s2s_outlet * state.segp); + } // all the control equations // TODO: can be optimized better @@ -1284,10 +1310,25 @@ namespace Opm { void BlackoilMultiSegmentModel::computeSegmentFluidProperties(const SolutionState& state) { + const int nw = wellsMultiSegment().size(); const int nseg_total = state.segp.size(); const int np = numPhases(); + if ( !wops_ms_.has_multisegment_wells ){ + // not sure if this is needed actually + // TODO: to check later if this is really necessary. + segment_mass_flow_rates_ = ADB::constant(V::Zero(nseg_total)); + well_segment_densities_ = ADB::constant(V::Zero(nseg_total)); + segment_mass_flow_rates_ = ADB::constant(V::Zero(nseg_total)); + segment_viscosities_ = ADB::constant(V::Zero(nseg_total)); + for (int phase = 0; phase < np; ++phase) { + segment_comp_surf_volume_current_[phase] = ADB::constant(V::Zero(nseg_total)); + segment_comp_surf_volume_initial_[phase] = V::Zero(nseg_total); + } + return; + } + // although we will calculate segment density for non-segmented wells at the same time, // while under most of the cases, they will not be used, // since for most of the cases, the density calculation for non-segment wells are @@ -1475,6 +1516,12 @@ namespace Opm { const int nw = wellsMultiSegment().size(); const int nseg_total = state.segp.size(); + if ( !wops_ms_.has_multisegment_wells ) { + well_segment_pressures_delta_ = ADB::constant(V::Zero(nseg_total)); + well_segment_perforation_pressure_diffs_ = wops_ms_.s2p * well_segment_pressures_delta_; + return; + } + // calculate the depth difference of the segments // TODO: we need to store the following values somewhere to avoid recomputation. V segment_depth_delta = V::Zero(nseg_total); From b4b821521e6a6c7607af6e9073bb65a59016d2b6 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Thu, 10 Mar 2016 14:41:10 +0800 Subject: [PATCH 07/24] enable relperm diagnostics for solvent model. --- opm/autodiff/FlowMainSolvent.hpp | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/opm/autodiff/FlowMainSolvent.hpp b/opm/autodiff/FlowMainSolvent.hpp index 5e1346996..eaac004a8 100644 --- a/opm/autodiff/FlowMainSolvent.hpp +++ b/opm/autodiff/FlowMainSolvent.hpp @@ -107,18 +107,6 @@ namespace Opm Base::deck_->hasKeyword("SOLVENT"))); } - // run diagnostics - // Writes to: - // logFile_ - void runDiagnostics() - { - std::shared_ptr streamLog = std::make_shared(Base::logFile_, Opm::Log::DefaultMessageTypes); - const std::string msg = "Relperm diagnostic not implemented for the solvent model "; - std::cout << "\n" << msg << "\n" << std::endl; - streamLog->addMessage(Log::MessageType::Info, msg); - } - - }; From f0a2959a6dff057bec8db3d2707f40c009a2a90d Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Tue, 15 Mar 2016 09:43:18 +0800 Subject: [PATCH 08/24] use EclipsePRTLog. --- opm/autodiff/FlowMain.hpp | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 527ec85ae..088a31aa4 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -71,8 +71,7 @@ #include #include -#include -#include +#include #include #include #include @@ -352,11 +351,8 @@ namespace Opm // Create Parser ParserPtr parser(new Parser()); { - std::shared_ptr streamLog = std::make_shared(logFile_ , Log::DefaultMessageTypes); - std::shared_ptr counterLog = std::make_shared(Log::DefaultMessageTypes); - - OpmLog::addBackend( "STREAM" , streamLog ); - OpmLog::addBackend( "COUNTER" , counterLog ); + std::shared_ptr prtLog = std::make_shared(logFile_ , Log::DefaultMessageTypes); + OpmLog::addBackend( "ECLIPSEPRTLOG" , prtLog ); } // Create Deck and EclipseState. From debf039175880701e1baf87beb34a58f961e1875 Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Thu, 17 Mar 2016 08:55:52 +0800 Subject: [PATCH 09/24] rename ParseMode as ParseContext. --- opm/autodiff/FlowMain.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 088a31aa4..b3095fd2b 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -74,7 +74,7 @@ #include #include #include -#include +#include #include #include #include @@ -357,10 +357,10 @@ namespace Opm // Create Deck and EclipseState. try { - ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); - deck_ = parser->parseFile(deck_filename, parseMode); + ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }}); + deck_ = parser->parseFile(deck_filename, parseContext); checkDeck(deck_, parser); - eclipse_state_.reset(new EclipseState(deck_, parseMode)); + eclipse_state_.reset(new EclipseState(deck_, parseContext)); } catch (const std::invalid_argument& e) { std::cerr << "Failed to create valid EclipseState object. See logfile: " << logFile_ << std::endl; From 17e5eb2dee051eece2e2bcd3c1b935509a46493f Mon Sep 17 00:00:00 2001 From: Liu Ming Date: Thu, 17 Mar 2016 10:02:50 +0800 Subject: [PATCH 10/24] rename ParseMode as ParseContext. --- examples/opm_init_check.cpp | 8 ++--- examples/sim_2p_incomp_ad.cpp | 8 ++--- examples/sim_poly2p_comp_reorder.cpp | 8 ++--- examples/sim_poly2p_incomp_reorder.cpp | 8 ++--- examples/sim_poly_fi2p_comp_ad.cpp | 8 ++--- tests/test_boprops_ad.cpp | 8 ++--- tests/test_rateconverter.cpp | 8 ++--- tests/test_transmissibilitymultipliers.cpp | 38 +++++++++++----------- tests/test_vfpproperties.cpp | 6 ++-- 9 files changed, 50 insertions(+), 50 deletions(-) diff --git a/examples/opm_init_check.cpp b/examples/opm_init_check.cpp index 233bbe2b9..a3ab8787a 100644 --- a/examples/opm_init_check.cpp +++ b/examples/opm_init_check.cpp @@ -31,7 +31,7 @@ #include #include -#include +#include #include #include #include @@ -345,10 +345,10 @@ int main(int argc, char** argv) { ParserPtr parser(new Parser()); - ParseMode parseMode; + ParseContext parseContext; std::cout << "Parsing input file ............: " << input_file << std::endl; - DeckConstPtr deck = parser->parseFile(input_file, parseMode); - std::shared_ptr state = std::make_shared( deck , parseMode ); + DeckConstPtr deck = parser->parseFile(input_file, parseContext); + std::shared_ptr state = std::make_shared( deck , parseContext ); std::cout << "Loading eclipse INIT file .....: " << init_file << std::endl; ecl_file_type * ecl_init = ecl_file_open( init_file.c_str() , 0 ); diff --git a/examples/sim_2p_incomp_ad.cpp b/examples/sim_2p_incomp_ad.cpp index dea22d666..e4e203228 100644 --- a/examples/sim_2p_incomp_ad.cpp +++ b/examples/sim_2p_incomp_ad.cpp @@ -47,7 +47,7 @@ #include #include -#include +#include #include #include @@ -112,9 +112,9 @@ try double gravity[3] = { 0.0 }; if (use_deck) { std::string deck_filename = param.get("deck_filename"); - Opm::ParseMode parseMode; - deck = parser->parseFile(deck_filename, parseMode); - eclipseState.reset(new EclipseState(deck , parseMode)); + Opm::ParseContext parseContext; + deck = parser->parseFile(deck_filename, parseContext); + eclipseState.reset(new EclipseState(deck , parseContext)); // Grid init grid.reset(new GridManager(deck)); diff --git a/examples/sim_poly2p_comp_reorder.cpp b/examples/sim_poly2p_comp_reorder.cpp index 72b84f827..175c7019b 100644 --- a/examples/sim_poly2p_comp_reorder.cpp +++ b/examples/sim_poly2p_comp_reorder.cpp @@ -48,7 +48,7 @@ #include #include -#include +#include #include #include @@ -100,9 +100,9 @@ try if (use_deck) { std::string deck_filename = param.get("deck_filename"); ParserPtr parser(new Opm::Parser()); - Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); - deck = parser->parseFile(deck_filename , parseMode); - eclipseState.reset(new Opm::EclipseState(deck , parseMode)); + Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }}); + deck = parser->parseFile(deck_filename , parseContext); + eclipseState.reset(new Opm::EclipseState(deck , parseContext)); // Grid init grid.reset(new GridManager(deck)); diff --git a/examples/sim_poly2p_incomp_reorder.cpp b/examples/sim_poly2p_incomp_reorder.cpp index 5ab15aa9b..cbb234eea 100644 --- a/examples/sim_poly2p_incomp_reorder.cpp +++ b/examples/sim_poly2p_incomp_reorder.cpp @@ -48,7 +48,7 @@ #include #include -#include +#include #include #include @@ -99,10 +99,10 @@ try double gravity[3] = { 0.0 }; if (use_deck) { std::string deck_filename = param.get("deck_filename"); - Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); + Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }}); ParserPtr parser(new Opm::Parser()); - deck = parser->parseFile(deck_filename , parseMode); - eclipseState.reset(new Opm::EclipseState(deck , parseMode)); + deck = parser->parseFile(deck_filename , parseContext); + eclipseState.reset(new Opm::EclipseState(deck , parseContext)); // Grid init grid.reset(new GridManager(deck)); diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp index 55a577787..4713fbcdd 100644 --- a/examples/sim_poly_fi2p_comp_ad.cpp +++ b/examples/sim_poly_fi2p_comp_ad.cpp @@ -66,7 +66,7 @@ #include #include #include -#include +#include #include #include @@ -141,7 +141,7 @@ try } std::string logFile = output_dir + "/LOGFILE.txt"; - Opm::ParseMode parseMode({{ ParseMode::PARSE_RANDOM_SLASH , InputError::IGNORE }}); + Opm::ParseContext parseContext({{ ParseContext::PARSE_RANDOM_SLASH , InputError::IGNORE }}); Opm::ParserPtr parser(new Opm::Parser()); { std::shared_ptr streamLog = std::make_shared(logFile , Opm::Log::DefaultMessageTypes); @@ -154,9 +154,9 @@ try Opm::DeckConstPtr deck; std::shared_ptr eclipseState; try { - deck = parser->parseFile(deck_filename , parseMode); + deck = parser->parseFile(deck_filename , parseContext); Opm::checkDeck(deck, parser); - eclipseState.reset(new Opm::EclipseState(deck , parseMode)); + eclipseState.reset(new Opm::EclipseState(deck , parseContext)); } catch (const std::invalid_argument& e) { std::cerr << "Failed to create valid ECLIPSESTATE object. See logfile: " << logFile << std::endl; diff --git a/tests/test_boprops_ad.cpp b/tests/test_boprops_ad.cpp index d4a036a7a..20f0a21cf 100644 --- a/tests/test_boprops_ad.cpp +++ b/tests/test_boprops_ad.cpp @@ -37,7 +37,7 @@ #include #include -#include +#include #include #include @@ -48,10 +48,10 @@ struct SetupSimple { SetupSimple() { - Opm::ParseMode parseMode; + Opm::ParseContext parseContext; Opm::ParserPtr parser(new Opm::Parser()); - deck = parser->parseFile("fluid.data", parseMode); - eclState.reset(new Opm::EclipseState(deck , parseMode)); + deck = parser->parseFile("fluid.data", parseContext); + eclState.reset(new Opm::EclipseState(deck , parseContext)); param.disableOutput(); param.insertParameter("init_rock" , "false" ); diff --git a/tests/test_rateconverter.cpp b/tests/test_rateconverter.cpp index 81dd53137..6c9e28898 100644 --- a/tests/test_rateconverter.cpp +++ b/tests/test_rateconverter.cpp @@ -39,7 +39,7 @@ #include #include -#include +#include #include #include @@ -48,9 +48,9 @@ struct SetupSimple { SetupSimple() { Opm::ParserPtr parser(new Opm::Parser()); - Opm::ParseMode parseMode; - deck = parser->parseFile("fluid.data" , parseMode); - eclState.reset(new Opm::EclipseState(deck , parseMode)); + Opm::ParseContext parseContext; + deck = parser->parseFile("fluid.data" , parseContext); + eclState.reset(new Opm::EclipseState(deck , parseContext)); param.disableOutput(); param.insertParameter("init_rock" , "false" ); diff --git a/tests/test_transmissibilitymultipliers.cpp b/tests/test_transmissibilitymultipliers.cpp index 90c77ca16..0e6ec29da 100644 --- a/tests/test_transmissibilitymultipliers.cpp +++ b/tests/test_transmissibilitymultipliers.cpp @@ -31,7 +31,7 @@ #include #include -#include +#include #include #include @@ -158,12 +158,12 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface) { Opm::parameter::ParameterGroup param; Opm::ParserPtr parser(new Opm::Parser() ); - Opm::ParseMode parseMode; + Opm::ParseContext parseContext; ///// // create a DerivedGeology object without any multipliers involved - Opm::DeckConstPtr origDeck = parser->parseString(origDeckString, parseMode); - Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseMode)); + Opm::DeckConstPtr origDeck = parser->parseString(origDeckString, parseContext); + Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseContext)); auto origGridManager = std::make_shared(origEclipseState->getEclipseGrid()); auto origProps = std::make_shared(origDeck, origEclipseState, *(origGridManager->c_grid())); @@ -173,8 +173,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface) ///// // create a DerivedGeology object _with_ transmissibility multipliers involved - Opm::DeckConstPtr multDeck = parser->parseString(multDeckString, parseMode); - Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseMode)); + Opm::DeckConstPtr multDeck = parser->parseString(multDeckString, parseContext); + Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseContext)); auto multGridManager = std::make_shared(multEclipseState->getEclipseGrid()); auto multProps = std::make_shared(multDeck, multEclipseState, *(multGridManager->c_grid())); @@ -185,8 +185,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface) ///// // create a DerivedGeology object _with_ transmissibility multipliers involved for // the negative faces - Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString, parseMode); - Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck , parseMode)); + Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString, parseContext); + Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck , parseContext)); auto multMinusGridManager = std::make_shared(multMinusEclipseState->getEclipseGrid()); auto multMinusProps = std::make_shared(multMinusDeck, multMinusEclipseState, *(multMinusGridManager->c_grid())); @@ -196,8 +196,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersLegacyGridInterface) ///// // create a DerivedGeology object with the NTG keyword involved - Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseMode); - Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseMode)); + Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseContext); + Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseContext)); auto ntgGridManager = std::make_shared(ntgEclipseState->getEclipseGrid()); auto ntgProps = std::make_shared(ntgDeck, ntgEclipseState, *(ntgGridManager->c_grid())); @@ -273,12 +273,12 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid) Opm::parameter::ParameterGroup param; Opm::ParserPtr parser(new Opm::Parser() ); - Opm::ParseMode parseMode; + Opm::ParseContext parseContext; ///// // create a DerivedGeology object without any multipliers involved - Opm::DeckConstPtr origDeck = parser->parseString(origDeckString , parseMode); - Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseMode)); + Opm::DeckConstPtr origDeck = parser->parseString(origDeckString , parseContext); + Opm::EclipseStateConstPtr origEclipseState(new Opm::EclipseState(origDeck , parseContext)); auto origGrid = std::make_shared(); origGrid->processEclipseFormat(origEclipseState->getEclipseGrid(), 0.0, false); @@ -292,8 +292,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid) ///// // create a DerivedGeology object _with_ transmissibility multipliers involved - Opm::DeckConstPtr multDeck = parser->parseString(multDeckString,parseMode); - Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseMode)); + Opm::DeckConstPtr multDeck = parser->parseString(multDeckString,parseContext); + Opm::EclipseStateConstPtr multEclipseState(new Opm::EclipseState(multDeck, parseContext)); auto multGrid = std::make_shared(); multGrid->processEclipseFormat(multEclipseState->getEclipseGrid(), 0.0, false); @@ -306,8 +306,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid) ///// // create a DerivedGeology object _with_ transmissibility multipliers involved for // the negative faces - Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString , parseMode); - Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck, parseMode)); + Opm::DeckConstPtr multMinusDeck = parser->parseString(multMinusDeckString , parseContext); + Opm::EclipseStateConstPtr multMinusEclipseState(new Opm::EclipseState(multMinusDeck, parseContext)); auto multMinusGrid = std::make_shared(); multMinusGrid->processEclipseFormat(multMinusEclipseState->getEclipseGrid(), 0.0, false); @@ -320,8 +320,8 @@ BOOST_AUTO_TEST_CASE(TransmissibilityMultipliersCpGrid) ///// // create a DerivedGeology object with the NTG keyword involved - Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseMode); - Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseMode)); + Opm::DeckConstPtr ntgDeck = parser->parseString(ntgDeckString, parseContext); + Opm::EclipseStateConstPtr ntgEclipseState(new Opm::EclipseState(ntgDeck, parseContext)); auto ntgGrid = std::make_shared(); ntgGrid->processEclipseFormat(ntgEclipseState->getEclipseGrid(), 0.0, false); diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 3a1b7a164..e8e1d9415 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -39,7 +39,7 @@ #include #include -#include +#include #include #include #include @@ -1047,7 +1047,7 @@ VFPPROD \n\ std::shared_ptr units(Opm::UnitSystem::newFIELD()); Opm::ParserPtr parser(new Opm::Parser()); - Opm::ParseMode parse_mode; + Opm::ParseContext parse_mode; deck = parser->parseString(table_str, parse_mode); BOOST_REQUIRE(deck->hasKeyword("VFPPROD")); @@ -1108,7 +1108,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD) std::shared_ptr units(Opm::UnitSystem::newMETRIC()); Opm::ParserPtr parser(new Opm::Parser()); - Opm::ParseMode parse_mode; + Opm::ParseContext parse_mode; boost::filesystem::path file("VFPPROD2"); deck = parser->parseFile(file.string(), parse_mode); From ea9e572a656c9ca1089723fe33a55b5e250b0de7 Mon Sep 17 00:00:00 2001 From: "Kjell W. Kongsvik" Date: Thu, 17 Mar 2016 14:45:30 +0100 Subject: [PATCH 11/24] Refactored to use OutputWriter from opm-output As OutputWriter has been deleted from opm-core Only changes in includes --- examples/sim_poly_fi2p_comp_ad.cpp | 2 +- opm/autodiff/SimulatorBase.hpp | 2 +- opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp | 2 +- opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp | 6 +++--- opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp | 2 +- opm/autodiff/SimulatorIncompTwophaseAd.cpp | 2 +- opm/polymer/SimulatorCompressiblePolymer.cpp | 2 +- opm/polymer/SimulatorPolymer.cpp | 4 ++-- .../fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp | 2 +- .../SimulatorFullyImplicitCompressiblePolymer.hpp | 4 ++-- 10 files changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp index 4713fbcdd..24f7ab0c5 100644 --- a/examples/sim_poly_fi2p_comp_ad.cpp +++ b/examples/sim_poly_fi2p_comp_ad.cpp @@ -37,7 +37,7 @@ #include #include -#include +#include #include #include #include diff --git a/opm/autodiff/SimulatorBase.hpp b/opm/autodiff/SimulatorBase.hpp index 1680ace03..50de46032 100644 --- a/opm/autodiff/SimulatorBase.hpp +++ b/opm/autodiff/SimulatorBase.hpp @@ -40,7 +40,7 @@ #include #include #include -#include +#include #include #include diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index a57fcdb8a..1bc37aa19 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -22,7 +22,7 @@ #include "SimulatorFullyImplicitBlackoilOutput.hpp" #include -#include +#include #include #include #include diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index a4e8e879c..9d09fd7f6 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -24,12 +24,12 @@ #include #include #include -#include +#include #include #include -#include -#include +#include +#include #include #include diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp index bad80a289..9ad4492db 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilSolvent.hpp @@ -44,7 +44,7 @@ #include //#include #include -#include +#include #include #include diff --git a/opm/autodiff/SimulatorIncompTwophaseAd.cpp b/opm/autodiff/SimulatorIncompTwophaseAd.cpp index 49c36a985..c4dbd47ee 100644 --- a/opm/autodiff/SimulatorIncompTwophaseAd.cpp +++ b/opm/autodiff/SimulatorIncompTwophaseAd.cpp @@ -37,7 +37,7 @@ #include #include #include -#include +#include #include #include diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp index a483cda58..eee04dde1 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.cpp +++ b/opm/polymer/SimulatorCompressiblePolymer.cpp @@ -34,7 +34,7 @@ #include #include #include -#include +#include #include #include diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp index 615bbd4eb..c6d54ea29 100644 --- a/opm/polymer/SimulatorPolymer.cpp +++ b/opm/polymer/SimulatorPolymer.cpp @@ -36,7 +36,7 @@ #include #include #include -#include +#include #include #include @@ -62,7 +62,7 @@ #include #ifdef HAVE_ERT -#include +#include #endif diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp index 091e22f8b..76d3bc7d0 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitBlackoilPolymer.hpp @@ -45,7 +45,7 @@ #include //#include #include -#include +#include #include #include diff --git a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp index fdae6d0d5..136cdc62c 100644 --- a/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp +++ b/opm/polymer/fullyimplicit/SimulatorFullyImplicitCompressiblePolymer.hpp @@ -38,8 +38,8 @@ #include #include #include -#include -#include +#include +#include #include #include From e63bf7aeadb5473c15fe78bd6e698e6c3254825b Mon Sep 17 00:00:00 2001 From: chflo Date: Thu, 17 Mar 2016 16:19:07 +0100 Subject: [PATCH 12/24] Changed includes due to OpmLog moved from parser to common --- examples/sim_poly_fi2p_comp_ad.cpp | 6 +++--- opm/autodiff/FlowMain.hpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp index 4713fbcdd..eb266ef6a 100644 --- a/examples/sim_poly_fi2p_comp_ad.cpp +++ b/examples/sim_poly_fi2p_comp_ad.cpp @@ -61,9 +61,9 @@ #include #include -#include -#include -#include +#include +#include +#include #include #include #include diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index b3095fd2b..b28703ffc 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -70,8 +70,8 @@ #include -#include -#include +#include +#include #include #include #include From 94aa360fe610273c5036f9c3cd1dd388238cf90a Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 18 Mar 2016 13:16:23 +0100 Subject: [PATCH 13/24] Bugfix. Fix setting initial rates in updateWellControls - The initial rates are only set to target values for single phase producers (orat, wrat, grat). - For injectors compi is used to determine the initial target rates. --- opm/autodiff/BlackoilModelBase_impl.hpp | 28 ++++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index c0b3fefad..a7c90b017 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1610,11 +1610,33 @@ namespace detail { break; case SURFACE_RATE: - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0) { - xw.wellRates()[np*w + phase] = target * distr[phase]; + // assign target value as initial guess for injectors and + // single phase producers (orat, grat, wrat) + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + for (int phase = 0; phase < np; ++phase) { + const double& compi = wells().comp_frac[np * w + phase]; + if (compi > 0.0) { + xw.wellRates()[np*w + phase] = target * compi; + } } + } else if (well_type == PRODUCER) { + // for single phase producers sum of distr should be 1.0 + double sumdistr = distr[0]; + for (int phase = 1; phase < np; ++phase) { + sumdistr += distr[phase]; + } + std::cout << sumdistr << std::endl; + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && sumdistr == 1.0 ) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); } + + break; } From d91831b971d94626930847565488bd17b8fb0c86 Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Fri, 18 Mar 2016 14:09:06 +0100 Subject: [PATCH 14/24] Avoid comparsion of floating point numbers --- opm/autodiff/BlackoilModelBase_impl.hpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index a7c90b017..48ab339a1 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1621,14 +1621,18 @@ namespace detail { } } } else if (well_type == PRODUCER) { - // for single phase producers sum of distr should be 1.0 - double sumdistr = distr[0]; - for (int phase = 1; phase < np; ++phase) { - sumdistr += distr[phase]; - } - std::cout << sumdistr << std::endl; + + // only set target as initial rates for single phase + // producers. (orat, grat and wrat, and not lrat) + // lrat will result in numPhasesWithTargetsUnderThisControl == 2 + int numPhasesWithTargetsUnderThisControl = 0; for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0 && sumdistr == 1.0 ) { + if (distr[phase] > 0.0) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { xw.wellRates()[np*w + phase] = target * distr[phase]; } } From 18c07d5d66ec7c5c4d3c733767877e794c3af0f3 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Sun, 28 Feb 2016 11:27:00 +0100 Subject: [PATCH 15/24] Replaced SimulatorState -> SimulationDataContainer --- CMakeLists_files.cmake | 3 + examples/sim_2p_incomp_ad.cpp | 67 ++++++----- examples/sim_poly2p_comp_reorder.cpp | 106 ++++++++++-------- examples/sim_poly2p_incomp_reorder.cpp | 101 ++++++++++------- examples/sim_poly_fi2p_comp_ad.cpp | 12 +- opm/autodiff/BackupRestore.hpp | 17 +-- opm/autodiff/BlackoilModelBase.hpp | 4 +- opm/autodiff/BlackoilModelBase_impl.hpp | 5 +- opm/autodiff/BlackoilSolventModel_impl.hpp | 14 ++- opm/autodiff/BlackoilSolventState.cpp | 34 ++++++ opm/autodiff/BlackoilSolventState.hpp | 23 +--- opm/autodiff/FlowMain.hpp | 46 +++++--- opm/autodiff/ParallelDebugOutput.hpp | 95 ++++++++-------- .../SimulatorFullyImplicitBlackoilOutput.cpp | 27 +++-- .../SimulatorFullyImplicitBlackoilOutput.hpp | 28 ++--- opm/autodiff/TransportSolverTwophaseAd.cpp | 1 + opm/polymer/CompressibleTpfaPolymer.cpp | 4 +- opm/polymer/IncompTpfaPolymer.cpp | 8 +- opm/polymer/IncompTpfaPolymer.hpp | 3 +- opm/polymer/PolymerBlackoilState.cpp | 42 +++++++ opm/polymer/PolymerBlackoilState.hpp | 23 +--- opm/polymer/PolymerState.cpp | 42 +++++++ opm/polymer/PolymerState.hpp | 31 +---- opm/polymer/SimulatorCompressiblePolymer.cpp | 18 +-- opm/polymer/SimulatorPolymer.cpp | 26 ++--- .../BlackoilPolymerModel_impl.hpp | 36 +++--- ...FullyImplicitCompressiblePolymerSolver.cpp | 46 +++++--- opm/polymer/polymerUtilities.cpp | 10 +- tests/test_rateconverter.cpp | 4 +- tests/test_singlecellsolves.cpp | 76 ++++++++----- 30 files changed, 569 insertions(+), 383 deletions(-) create mode 100644 opm/autodiff/BlackoilSolventState.cpp create mode 100644 opm/polymer/PolymerBlackoilState.cpp create mode 100644 opm/polymer/PolymerState.cpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f0b3cd2f7..d18971f0c 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -47,6 +47,9 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/VFPProdProperties.cpp opm/autodiff/VFPInjProperties.cpp opm/autodiff/WellMultiSegment.cpp + opm/autodiff/BlackoilSolventState.cpp + opm/polymer/PolymerState.cpp + opm/polymer/PolymerBlackoilState.cpp opm/polymer/CompressibleTpfaPolymer.cpp opm/polymer/IncompTpfaPolymer.cpp opm/polymer/PolymerInflow.cpp diff --git a/examples/sim_2p_incomp_ad.cpp b/examples/sim_2p_incomp_ad.cpp index e4e203228..8e671a6ee 100644 --- a/examples/sim_2p_incomp_ad.cpp +++ b/examples/sim_2p_incomp_ad.cpp @@ -105,8 +105,9 @@ try std::unique_ptr grid; std::unique_ptr props; std::unique_ptr rock_comp; + std::unique_ptr state; EclipseStateConstPtr eclipseState; - TwophaseState state; + // bool check_well_controls = false; // int max_well_control_iterations = 0; double gravity[3] = { 0.0 }; @@ -118,19 +119,25 @@ try // Grid init grid.reset(new GridManager(deck)); - // Rock and fluid init - props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid())); - // 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); - } else { - initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + { + const UnstructuredGrid& ug_grid = *(grid->c_grid()); + + // Rock and fluid init + props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid)); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + + state.reset( new TwophaseState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ))); + // 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(ug_grid, *props, param, gravity[2], *state); + } else { + initStateFromDeck(ug_grid, *props, deck, gravity[2], *state); + } } } else { // Grid init. @@ -141,14 +148,20 @@ try const double dy = param.getDefault("dy", 1.0); const double dz = param.getDefault("dz", 1.0); grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); - // Rock and fluid init. - props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - // Rock compressibility. - rock_comp.reset(new RockCompressibility(param)); - // Gravity. - gravity[2] = param.getDefault("gravity", 0.0); - // Init state variables (saturation and pressure). - initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); + { + const UnstructuredGrid& ug_grid = *(grid->c_grid()); + + // Rock and fluid init. + props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid ))); + + state.reset( new TwophaseState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ))); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(ug_grid, *props, param, gravity[2], *state); + } } // Warn if gravity but no density difference. @@ -170,7 +183,7 @@ try // terms of total pore volume. std::vector porevol; if (rock_comp->isActive()) { - computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol); } else { computePorevolume(*grid->c_grid(), props->porosity(), porevol); } @@ -236,8 +249,8 @@ try simtimer.init(param); warnIfUnusedParams(param); WellState well_state; - well_state.init(0, state); - rep = simulator.run(simtimer, state, well_state); + well_state.init(0, *state); + rep = simulator.run(simtimer, *state, well_state); } else { // With a deck, we may have more report steps etc. WellState well_state; @@ -255,7 +268,7 @@ try // properly transfer old well state to it every report step, // since number of wells may change etc. if (reportStepIdx == 0) { - well_state.init(wells.c_wells(), state); + well_state.init(wells.c_wells(), *state); } simtimer.setCurrentStepNum(reportStepIdx); @@ -273,7 +286,7 @@ try if (reportStepIdx == 0) { warnIfUnusedParams(param); } - SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); + SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state); if (output) { epoch_rep.reportParam(epoch_os); } diff --git a/examples/sim_poly2p_comp_reorder.cpp b/examples/sim_poly2p_comp_reorder.cpp index 175c7019b..b4fd352ce 100644 --- a/examples/sim_poly2p_comp_reorder.cpp +++ b/examples/sim_poly2p_comp_reorder.cpp @@ -92,7 +92,7 @@ try boost::scoped_ptr rock_comp; Opm::DeckConstPtr deck; EclipseStateConstPtr eclipseState; - PolymerBlackoilState state; + std::unique_ptr state; Opm::PolymerProperties poly_props; // bool check_well_controls = false; // int max_well_control_iterations = 0; @@ -106,23 +106,29 @@ try // Grid init grid.reset(new GridManager(deck)); - // Rock and fluid init - props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, *grid->c_grid())); - // 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); - } else { - initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + { + const UnstructuredGrid& ug_grid = *(grid->c_grid()); + + // Rock and fluid init + props.reset(new BlackoilPropertiesFromDeck(deck, eclipseState, ug_grid)); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + + state.reset( new PolymerBlackoilState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2)); + // 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(ug_grid, *props, param, gravity[2], *state); + } else { + initStateFromDeck(ug_grid, *props, deck, gravity[2], *state); + } + initBlackoilSurfvol(ug_grid, *props, *state); + // Init polymer properties. + poly_props.readFromDeck(deck, eclipseState); } - initBlackoilSurfvol(*grid->c_grid(), *props, state); - // Init polymer properties. - poly_props.readFromDeck(deck, eclipseState); } else { // Grid init. const int nx = param.getDefault("nx", 100); @@ -132,31 +138,43 @@ try const double dy = param.getDefault("dy", 1.0); const double dz = param.getDefault("dz", 1.0); grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); - // Rock and fluid init. - props.reset(new BlackoilPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - // Rock compressibility. - rock_comp.reset(new RockCompressibility(param)); - // Gravity. - gravity[2] = param.getDefault("gravity", 0.0); - // Init state variables (saturation and pressure). - initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); - initBlackoilSurfvol(*grid->c_grid(), *props, state); - // Init Polymer state - if (param.has("poly_init")) { - double poly_init = param.getDefault("poly_init", 0.0); - for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) { - double smin[2], smax[2]; - props->satRange(1, &cell, smin, smax); - if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) { - state.concentration()[cell] = poly_init; - state.maxconcentration()[cell] = poly_init; - } else { - state.saturation()[2*cell + 0] = 0.; - state.saturation()[2*cell + 1] = 1.; - state.concentration()[cell] = 0.; - state.maxconcentration()[cell] = 0.; + { + const UnstructuredGrid& ug_grid = *(grid->c_grid()); + + // Rock and fluid init. + props.reset(new BlackoilPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid ))); + state.reset( new PolymerBlackoilState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2)); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(ug_grid, *props, param, gravity[2], *state); + initBlackoilSurfvol(ug_grid, *props, *state); + // Init Polymer state + + if (param.has("poly_init")) { + double poly_init = param.getDefault("poly_init", 0.0); + for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) { + double smin[2], smax[2]; + + auto& saturation = state->saturation(); + auto& concentration = state->getCellData( state->CONCENTRATION ); + auto& max_concentration = state->getCellData( state->CMAX ); + + props->satRange(1, &cell, smin, smax); + if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) { + concentration[cell] = poly_init; + max_concentration[cell] = poly_init; + } else { + saturation[2*cell + 0] = 0.; + saturation[2*cell + 1] = 1.; + concentration[cell] = 0.; + max_concentration[cell] = 0.; + } } } + } // Init polymer properties. // Setting defaults to provide a simple example case. @@ -240,8 +258,8 @@ try simtimer.init(param); warnIfUnusedParams(param); WellState well_state; - well_state.init(0, state); - rep = simulator.run(simtimer, state, well_state); + well_state.init(0, *state); + rep = simulator.run(simtimer, *state, well_state); } else { // With a deck, we may have more epochs etc. WellState well_state; @@ -284,7 +302,7 @@ try // properly transfer old well state to it every report step, // since number of wells may change etc. if (reportStepIdx == 0) { - well_state.init(wells.c_wells(), state); + well_state.init(wells.c_wells(), *state); } // Create and run simulator. @@ -300,7 +318,7 @@ try if (reportStepIdx == 0) { warnIfUnusedParams(param); } - SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); + SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state); // Update total timing report and remember step number. rep += epoch_rep; diff --git a/examples/sim_poly2p_incomp_reorder.cpp b/examples/sim_poly2p_incomp_reorder.cpp index cbb234eea..3826a06d8 100644 --- a/examples/sim_poly2p_incomp_reorder.cpp +++ b/examples/sim_poly2p_incomp_reorder.cpp @@ -92,7 +92,7 @@ try boost::scoped_ptr props; boost::scoped_ptr rock_comp; EclipseStateConstPtr eclipseState; - PolymerState state; + std::unique_ptr state; Opm::PolymerProperties poly_props; // bool check_well_controls = false; // int max_well_control_iterations = 0; @@ -106,22 +106,27 @@ try // Grid init grid.reset(new GridManager(deck)); - // Rock and fluid init - props.reset(new IncompPropertiesFromDeck(deck, eclipseState, *grid->c_grid())); - // 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); - } else { - initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + { + const UnstructuredGrid& ug_grid = *(grid->c_grid()); + // Rock and fluid init + props.reset(new IncompPropertiesFromDeck(deck, eclipseState, ug_grid )); + // check_well_controls = param.getDefault("check_well_controls", false); + // max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); + state.reset( new PolymerState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2)); + + // 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(ug_grid, *props, param, gravity[2], *state); + } else { + initStateFromDeck(ug_grid, *props, deck, gravity[2], *state); + } + // Init polymer properties. + poly_props.readFromDeck(deck, eclipseState); } - // Init polymer properties. - poly_props.readFromDeck(deck, eclipseState); } else { // Grid init. const int nx = param.getDefault("nx", 100); @@ -131,28 +136,38 @@ try const double dy = param.getDefault("dy", 1.0); const double dz = param.getDefault("dz", 1.0); grid.reset(new GridManager(nx, ny, nz, dx, dy, dz)); - // Rock and fluid init. - props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - // Rock compressibility. - rock_comp.reset(new RockCompressibility(param)); - // Gravity. - gravity[2] = param.getDefault("gravity", 0.0); - // Init state variables (saturation and pressure). - initStateBasic(*grid->c_grid(), *props, param, gravity[2], state); - // Init Polymer state - if (param.has("poly_init")) { - double poly_init = param.getDefault("poly_init", 0.0); - for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) { - double smin[2], smax[2]; - props->satRange(1, &cell, smin, smax); - if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) { - state.concentration()[cell] = poly_init; - state.maxconcentration()[cell] = poly_init; - } else { - state.saturation()[2*cell + 0] = 0.; - state.saturation()[2*cell + 1] = 1.; - state.concentration()[cell] = 0.; - state.maxconcentration()[cell] = 0.; + { + const UnstructuredGrid& ug_grid = *(grid->c_grid()); + + // Rock and fluid init. + props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid )));; + state.reset( new PolymerState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ) , 2)); + // Rock compressibility. + rock_comp.reset(new RockCompressibility(param)); + // Gravity. + gravity[2] = param.getDefault("gravity", 0.0); + // Init state variables (saturation and pressure). + initStateBasic(ug_grid, *props, param, gravity[2], *state); + // Init Polymer state + if (param.has("poly_init")) { + double poly_init = param.getDefault("poly_init", 0.0); + for (int cell = 0; cell < UgGridHelpers::numCells( ug_grid ); ++cell) { + double smin[2], smax[2]; + + auto& saturation = state->saturation(); + auto& concentration = state->getCellData( state->CONCENTRATION ); + auto& max_concentration = state->getCellData( state->CMAX ); + + props->satRange(1, &cell, smin, smax); + if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) { + concentration[cell] = poly_init; + max_concentration[cell] = poly_init; + } else { + saturation[2*cell + 0] = 0.; + saturation[2*cell + 1] = 1.; + concentration[cell] = 0.; + max_concentration[cell] = 0.; + } } } } @@ -212,7 +227,7 @@ try // terms of total pore volume. std::vector porevol; if (rock_comp->isActive()) { - computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state.pressure(), porevol); + computePorevolume(*grid->c_grid(), props->porosity(), *rock_comp, state->pressure(), porevol); } else { computePorevolume(*grid->c_grid(), props->porosity(), porevol); } @@ -276,8 +291,8 @@ try simtimer.init(param); warnIfUnusedParams(param); WellState well_state; - well_state.init(0, state); - rep = simulator.run(simtimer, state, well_state); + well_state.init(0, *state); + rep = simulator.run(simtimer, *state, well_state); } else { // With a deck, we may have more epochs etc. @@ -321,7 +336,7 @@ try // properly transfer old well state to it every report step, // since number of wells may change etc. if (reportStepIdx == 0) { - well_state.init(wells.c_wells(), state); + well_state.init(wells.c_wells(), *state); } // Create and run simulator. @@ -339,7 +354,7 @@ try if (reportStepIdx == 0) { warnIfUnusedParams(param); } - SimulatorReport epoch_rep = simulator.run(simtimer, state, well_state); + SimulatorReport epoch_rep = simulator.run(simtimer, *state, well_state); // Update total timing report and remember step number. rep += epoch_rep; diff --git a/examples/sim_poly_fi2p_comp_ad.cpp b/examples/sim_poly_fi2p_comp_ad.cpp index eb266ef6a..4b17c59c7 100644 --- a/examples/sim_poly_fi2p_comp_ad.cpp +++ b/examples/sim_poly_fi2p_comp_ad.cpp @@ -115,7 +115,7 @@ try std::shared_ptr props; std::shared_ptr new_props; std::shared_ptr rock_comp; - PolymerBlackoilState state; + std::unique_ptr state; // bool check_well_controls = false; // int max_well_control_iterations = 0; double gravity[3] = { 0.0 }; @@ -187,6 +187,8 @@ try Opm::UgGridHelpers::globalCell(cGrid), Opm::UgGridHelpers::cartDims(cGrid), param)); + + state.reset( new PolymerBlackoilState( Opm::UgGridHelpers::numCells(cGrid), Opm::UgGridHelpers::numFaces(cGrid), 2)); new_props.reset(new BlackoilPropsAdFromDeck(deck, eclipseState, materialLawManager, cGrid)); PolymerProperties polymer_props(deck, eclipseState); PolymerPropsAd polymer_props_ad(polymer_props); @@ -199,10 +201,10 @@ try // 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); + initStateBasic(*grid->c_grid(), *props, param, gravity[2], *state); + initBlackoilSurfvol(*grid->c_grid(), *props, *state); } else { - initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state); + initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], *state); } bool use_gravity = (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0); @@ -254,7 +256,7 @@ try deck, *fis_solver, grav); - fullReport= simulator.run(simtimer, state); + fullReport= simulator.run(simtimer, *state); std::cout << "\n\n================ End of simulation ===============\n\n"; fullReport.report(std::cout); diff --git a/opm/autodiff/BackupRestore.hpp b/opm/autodiff/BackupRestore.hpp index 70fe4e9b6..0885257a3 100644 --- a/opm/autodiff/BackupRestore.hpp +++ b/opm/autodiff/BackupRestore.hpp @@ -21,8 +21,9 @@ #define OPM_BACKUPRESTORE_HEADER_INCLUDED #include +#include + -#include #include #include @@ -159,7 +160,7 @@ namespace Opm { BlackoilStateId = 2, WellStateFullyImplicitBackoilId = 3 }; - inline int objectId( const SimulatorState& /* state */) { + inline int objectId( const SimulationDataContainer& /* state */) { return SimulatorStateId; } inline int objectId( const WellState& /* state */) { @@ -184,9 +185,9 @@ namespace Opm { } } - // SimulatorState + // SimulationDataContainer inline - std::ostream& operator << (std::ostream& out, const SimulatorState& state ) + std::ostream& operator << (std::ostream& out, const SimulationDataContainer& state ) { // write id of object to stream writeValue( out, objectId( state ) ); @@ -205,7 +206,7 @@ namespace Opm { } inline - std::istream& operator >> (std::istream& in, SimulatorState& state ) + std::istream& operator >> (std::istream& in, SimulationDataContainer& state ) { // check id of stored object checkObjectId( in, state ); @@ -213,7 +214,7 @@ namespace Opm { int numPhases = 0; readValue( in, numPhases ); - if( numPhases != state.numPhases() ) + if( numPhases != (int) state.numPhases() ) OPM_THROW(std::logic_error,"num phases wrong"); // read variables @@ -234,7 +235,7 @@ namespace Opm { writeValue( out, objectId( state ) ); // backup simulator state - const SimulatorState& simstate = static_cast< const SimulatorState& > (state); + const SimulationDataContainer& simstate = static_cast< const SimulationDataContainer& > (state); out << simstate; // backup additional blackoil state variables @@ -252,7 +253,7 @@ namespace Opm { checkObjectId( in, state ); // restore simulator state - SimulatorState& simstate = static_cast< SimulatorState& > (state); + SimulationDataContainer& simstate = static_cast< SimulationDataContainer& > (state); in >> simstate; // restore additional blackoil state variables diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 9eceb3565..2551bd7fb 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -45,7 +45,7 @@ namespace Opm { class RockCompressibility; class NewtonIterationBlackoilInterface; class VFPProperties; - + class SimulationDataContainer; /// Struct for containing iteration variables. struct DefaultBlackoilSolutionState @@ -207,7 +207,7 @@ namespace Opm { /// \brief compute the relative change between to simulation states // \return || u^n+1 - u^n || / || u^n+1 || - double relativeChange( const SimulatorState& previous, const SimulatorState& current ) const; + double relativeChange( const SimulationDataContainer& previous, const SimulationDataContainer& current ) const; /// The size (number of unknowns) of the nonlinear system of equations. int sizeNonLinear() const; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index c0b3fefad..ad8a32562 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -48,6 +48,7 @@ #include #include +#include #include #include #include @@ -2520,8 +2521,8 @@ namespace detail { template double BlackoilModelBase:: - relativeChange(const SimulatorState& previous, - const SimulatorState& current ) const + relativeChange(const SimulationDataContainer& previous, + const SimulationDataContainer& current ) const { std::vector< double > p0 ( previous.pressure() ); std::vector< double > sat0( previous.saturation() ); diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 3f703a045..fb4688caf 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -137,11 +137,14 @@ namespace Opm { // Initial polymer concentration. if (has_solvent_) { - assert (not x.solvent_saturation().empty()); - const int nc = x.solvent_saturation().size(); - const V ss = Eigen::Map(&x.solvent_saturation()[0], nc); + const auto& solvent_saturation = x.getCellData( BlackoilSolventState::SSOL ); + const int nc = solvent_saturation.size(); + const V ss = Eigen::Map(solvent_saturation.data() , nc); + + // This is some insanely detailed flickety flackety code; // Solvent belongs after other reservoir vars but before well vars. auto solvent_pos = vars0.begin() + fluid_.numPhases(); + assert (not solvent_saturation.empty()); assert(solvent_pos == vars0.end() - 2); vars0.insert(solvent_pos, ss); } @@ -548,9 +551,10 @@ namespace Opm { Base::updateState(modified_dx, reservoir_state, well_state); // Update solvent. - const V ss_old = Eigen::Map(&reservoir_state.solvent_saturation()[0], nc, 1); + auto& solvent_saturation = reservoir_state.getCellData( reservoir_state.SSOL ); + const V ss_old = Eigen::Map(&solvent_saturation[0], nc, 1); const V ss = (ss_old - dss).max(zero); - std::copy(&ss[0], &ss[0] + nc, reservoir_state.solvent_saturation().begin()); + std::copy(&ss[0], &ss[0] + nc, solvent_saturation.begin()); // adjust oil saturation const Opm::PhaseUsage& pu = fluid_.phaseUsage(); diff --git a/opm/autodiff/BlackoilSolventState.cpp b/opm/autodiff/BlackoilSolventState.cpp new file mode 100644 index 000000000..681bd1ca2 --- /dev/null +++ b/opm/autodiff/BlackoilSolventState.cpp @@ -0,0 +1,34 @@ +/* + Copyright 2015 IRIS AS, Applied Mathematics. + + 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 + +namespace Opm +{ + const std::string BlackoilSolventState::SSOL = "SSOL"; + + BlackoilSolventState::BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases) + : BlackoilState( number_of_cells , number_of_faces , number_of_phases) + { + registerCellData( SSOL , 1 ); + }; + +} // namespace Opm + diff --git a/opm/autodiff/BlackoilSolventState.hpp b/opm/autodiff/BlackoilSolventState.hpp index cf5246443..0103e0904 100644 --- a/opm/autodiff/BlackoilSolventState.hpp +++ b/opm/autodiff/BlackoilSolventState.hpp @@ -20,9 +20,9 @@ #ifndef OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED #define OPM_BLACKOILSOLVENTSTATE_HEADER_INCLUDED -#include -#include -#include +#include + + #include namespace Opm @@ -33,22 +33,9 @@ namespace Opm class BlackoilSolventState : public BlackoilState { public: - void init(const UnstructuredGrid& g, int num_phases) - { - this->init(g.number_of_cells, g.number_of_faces, num_phases); - } + static const std::string SSOL; - void init(int number_of_cells, int number_of_faces, int num_phases) - { - BlackoilState::init(number_of_cells, number_of_faces, num_phases); - solventId_ = SimulatorState::registerCellData( "SSOL", 1 ); - } - - std::vector& solvent_saturation() { return cellData()[ solventId_ ]; } - const std::vector& solvent_saturation() const { return cellData()[ solventId_ ]; } - - private: - int solventId_; + BlackoilSolventState( int number_of_cells, int number_of_faces, int number_of_phases); }; } // namespace Opm diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index b28703ffc..1f5c7f4f0 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -183,7 +183,8 @@ namespace Opm bool use_local_perm_ = true; std::unique_ptr geoprops_; // setupState() - ReservoirState state_; + std::unique_ptr state_; + std::vector threshold_pressures_; // distributeData() boost::any parallel_information_; @@ -441,8 +442,13 @@ namespace Opm Opm::UgGridHelpers::cartDims(grid), param_); + // Init state variables (saturation and pressure). if (param_.has("init_saturation")) { + state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases() )); + initStateBasic(Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid), Opm::UgGridHelpers::cartDims(grid), @@ -451,25 +457,31 @@ namespace Opm Opm::UgGridHelpers::beginFaceCentroids(grid), Opm::UgGridHelpers::beginCellCentroids(grid), Opm::UgGridHelpers::dimensions(grid), - props, param_, gravity_[2], state_); + props, param_, gravity_[2], *state_); - initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, state_); + initBlackoilSurfvol(Opm::UgGridHelpers::numCells(grid), props, *state_); enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour }; if (pu.phase_used[Oil] && pu.phase_used[Gas]) { const int numPhases = props.numPhases(); const int numCells = Opm::UgGridHelpers::numCells(grid); + + // Uglyness 1: The state is a templated type, here we however make explicit use BlackoilState. + auto& gor = state_->getCellData( BlackoilState::GASOILRATIO ); + const auto& surface_vol = state_->getCellData( BlackoilState::SURFACEVOL ); for (int c = 0; c < numCells; ++c) { - state_.gasoilratio()[c] = state_.surfacevol()[c*numPhases + pu.phase_pos[Gas]] - / state_.surfacevol()[c*numPhases + pu.phase_pos[Oil]]; + // Uglyness 2: Here we explicitly use the layout of the saturation in the surface_vol field. + gor[c] = surface_vol[ c * numPhases + pu.phase_pos[Gas]] / surface_vol[ c * numPhases + pu.phase_pos[Oil]]; } } } else if (deck_->hasKeyword("EQUIL") && props.numPhases() == 3) { - state_.init(Opm::UgGridHelpers::numCells(grid), - Opm::UgGridHelpers::numFaces(grid), - props.numPhases()); - initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], state_); - state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); + // Which state class are we really using - what a f... mess? + state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases())); + + initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], *state_); + //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); } else { initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid), @@ -478,12 +490,12 @@ namespace Opm Opm::UgGridHelpers::beginFaceCentroids(grid), Opm::UgGridHelpers::beginCellCentroids(grid), Opm::UgGridHelpers::dimensions(grid), - props, deck_, gravity_[2], state_); + props, deck_, gravity_[2], *state_); } // Threshold pressures. std::map, double> maxDp; - computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), state_, props, gravity_[2]); + computeMaxDp(maxDp, deck_, eclipse_state_, grid_init_->grid(), *state_, props, gravity_[2]); threshold_pressures_ = thresholdPressures(deck_, eclipse_state_, grid, maxDp); std::vector threshold_pressures_nnc = thresholdPressuresNNC(eclipse_state_, geoprops_->nnc(), maxDp); threshold_pressures_.insert(threshold_pressures_.end(), threshold_pressures_nnc.begin(), threshold_pressures_nnc.end()); @@ -493,9 +505,9 @@ namespace Opm const int numCells = Opm::UgGridHelpers::numCells(grid); std::vector cells(numCells); for (int c = 0; c < numCells; ++c) { cells[c] = c; } - std::vector pc = state_.saturation(); - props.capPress(numCells, state_.saturation().data(), cells.data(), pc.data(), nullptr); - fluidprops_->setSwatInitScaling(state_.saturation(), pc); + std::vector pc = state_->saturation(); + props.capPress(numCells, state_->saturation().data(), cells.data(), pc.data(), nullptr); + fluidprops_->setSwatInitScaling(state_->saturation(), pc); } } @@ -518,7 +530,7 @@ namespace Opm // and initilialize new properties and states for it. if (must_distribute_) { distributeGridAndData(grid_init_->grid(), deck_, eclipse_state_, - state_, *fluidprops_, *geoprops_, + *state_, *fluidprops_, *geoprops_, material_law_manager_, threshold_pressures_, parallel_information_, use_local_perm_); } @@ -617,7 +629,7 @@ namespace Opm << std::flush; } - SimulatorReport fullReport = simulator_->run(simtimer, state_); + SimulatorReport fullReport = simulator_->run(simtimer, *state_); if (output_cout_) { std::cout << "\n\n================ End of simulation ===============\n\n"; diff --git a/opm/autodiff/ParallelDebugOutput.hpp b/opm/autodiff/ParallelDebugOutput.hpp index 2398a6f28..dabe9b4be 100644 --- a/opm/autodiff/ParallelDebugOutput.hpp +++ b/opm/autodiff/ParallelDebugOutput.hpp @@ -19,8 +19,10 @@ #ifndef OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED #define OPM_PARALLELDEBUGOUTPUT_HEADER_INCLUDED +#include + + #include -#include #include #include @@ -35,6 +37,7 @@ namespace Opm { class ParallelDebugOutputInterface + { protected: ParallelDebugOutputInterface () {} @@ -42,11 +45,11 @@ namespace Opm virtual ~ParallelDebugOutputInterface() {} // gather solution to rank 0 for EclipseWriter - virtual bool collectToIORank( const SimulatorState& localReservoirState, + virtual bool collectToIORank( const SimulationDataContainer& localReservoirState, const WellState& localWellState, const int reportStep ) = 0; - virtual const SimulatorState& globalReservoirState() const = 0 ; + virtual const SimulationDataContainer& globalReservoirState() const = 0 ; virtual const WellState& globalWellState() const = 0 ; virtual bool isIORank() const = 0; virtual bool isParallel() const = 0; @@ -60,7 +63,7 @@ namespace Opm protected: const GridImpl& grid_; - const SimulatorState* globalState_; + const SimulationDataContainer* globalState_; const WellState* wellState_; public: @@ -71,7 +74,7 @@ namespace Opm : grid_( grid ) {} // gather solution to rank 0 for EclipseWriter - virtual bool collectToIORank( const SimulatorState& localReservoirState, + virtual bool collectToIORank( const SimulationDataContainer& localReservoirState, const WellState& localWellState, const int /* reportStep */) { @@ -80,7 +83,7 @@ namespace Opm return true ; } - virtual const SimulatorState& globalReservoirState() const { return *globalState_; } + virtual const SimulationDataContainer& globalReservoirState() const { return *globalState_; } virtual const WellState& globalWellState() const { return *wellState_; } virtual bool isIORank () const { return true; } virtual bool isParallel () const { return false; } @@ -243,7 +246,7 @@ namespace Opm Dune::CpGrid& globalGrid = *grid_; // initialize global state with correct sizes - globalReservoirState_.init( globalGrid.numCells(), globalGrid.numFaces(), numPhases ); + globalReservoirState_.reset( new SimulationDataContainer( globalGrid.numCells(), globalGrid.numFaces(), numPhases )); // copy global cartesian index globalIndex_ = globalGrid.globalCell(); @@ -306,18 +309,18 @@ namespace Opm } } - class PackUnPackSimulatorState : public P2PCommunicatorType::DataHandleInterface + class PackUnPackSimulationDataContainer : public P2PCommunicatorType::DataHandleInterface { - const SimulatorState& localState_; - SimulatorState& globalState_; + const SimulationDataContainer& localState_; + SimulationDataContainer& globalState_; const WellState& localWellState_; WellState& globalWellState_; const IndexMapType& localIndexMap_; const IndexMapStorageType& indexMaps_; public: - PackUnPackSimulatorState( const SimulatorState& localState, - SimulatorState& globalState, + PackUnPackSimulationDataContainer( const SimulationDataContainer& localState, + SimulationDataContainer& globalState, const WellState& localWellState, WellState& globalWellState, const IndexMapType& localIndexMap, @@ -333,12 +336,11 @@ namespace Opm if( isIORank ) { // add missing data to global state - for( size_t i=globalState_.cellData().size(); - i& data = localState_.cellData()[ d ]; - const size_t stride = data.size() / numCells ; - assert( numCells * stride == data.size() ); + for (const auto& pair : localState_.cellData()) { + const std::string& key = pair.first; + const auto& data = pair.second; + const size_t stride = localState_.numCellDataComponents( key ); for( size_t i=0; i& data = globalState_.cellData()[ d ]; - const size_t stride = data.size() / numCells ; - assert( numCells * stride == data.size() ); + // write all cell data registered in local state + for (auto& pair : globalState_.cellData()) { + const std::string& key = pair.first; + auto& data = pair.second; + const size_t stride = globalState_.numCellDataComponents( key ); for( size_t i=0; i grid_; - Opm::EclipseStateConstPtr eclipseState_; - const double* permeability_; - P2PCommunicatorType toIORankComm_; - IndexMapType globalIndex_; - IndexMapType localIndexMap_; - IndexMapStorageType indexMaps_; - SimulatorState globalReservoirState_; + std::unique_ptr< Dune::CpGrid > grid_; + Opm::EclipseStateConstPtr eclipseState_; + const double* permeability_; + P2PCommunicatorType toIORankComm_; + IndexMapType globalIndex_; + IndexMapType localIndexMap_; + IndexMapStorageType indexMaps_; + std::unique_ptr globalReservoirState_; // this needs to be revised - WellStateFullyImplicitBlackoil globalWellState_; + WellStateFullyImplicitBlackoil globalWellState_; // true if we are on I/O rank - const bool isIORank_; + const bool isIORank_; }; #endif // #if HAVE_DUNE_CORNERPOINT diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index a57fcdb8a..d5f32b21c 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -21,6 +21,11 @@ #include "SimulatorFullyImplicitBlackoilOutput.hpp" +#include + +#include + +#include #include #include #include @@ -30,9 +35,6 @@ #include #include -#include - - #include #include #include @@ -50,7 +52,7 @@ namespace Opm void outputStateVtk(const UnstructuredGrid& grid, - const SimulatorState& state, + const SimulationDataContainer& state, const int step, const std::string& output_dir) { @@ -87,16 +89,15 @@ namespace Opm void outputStateMatlab(const UnstructuredGrid& grid, - const Opm::SimulatorState& state, + const Opm::SimulationDataContainer& state, const int step, const std::string& output_dir) { Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - for( unsigned int i=0; i cell_velocity; @@ -204,7 +205,7 @@ namespace Opm #ifdef HAVE_DUNE_CORNERPOINT void outputStateVtk(const Dune::CpGrid& grid, - const Opm::SimulatorState& state, + const Opm::SimulationDataContainer& state, const int step, const std::string& output_dir) { @@ -255,7 +256,7 @@ namespace Opm void BlackoilOutputWriter:: writeTimeStep(const SimulatorTimerInterface& timer, - const SimulatorState& localState, + const SimulationDataContainer& localState, const WellState& localWellState, bool substep) { @@ -271,7 +272,7 @@ namespace Opm isIORank = parallelOutput_->collectToIORank( localState, localWellState, timer.reportStepNum() ); } - const SimulatorState& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState; + const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState; const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState; // output is only done on I/O rank @@ -306,6 +307,7 @@ namespace Opm // write resport step number backupfile_.write( (const char *) &reportStep, sizeof(int) ); + /* try { const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state); backupfile_ << boState; @@ -317,6 +319,7 @@ namespace Opm { } + */ /* const WellStateFullyImplicitBlackoil* boWellState = diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index a4e8e879c..4e51690e6 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -20,7 +20,6 @@ #ifndef OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED #define OPM_SIMULATORFULLYIMPLICITBLACKOILOUTPUT_HEADER_INCLUDED #include -#include #include #include #include @@ -53,14 +52,17 @@ namespace Opm { + class SimulationDataContainer; + class BlackoilState; + void outputStateVtk(const UnstructuredGrid& grid, - const Opm::SimulatorState& state, + const Opm::SimulationDataContainer& state, const int step, const std::string& output_dir); void outputStateMatlab(const UnstructuredGrid& grid, - const Opm::SimulatorState& state, + const Opm::SimulationDataContainer& state, const int step, const std::string& output_dir); @@ -69,23 +71,23 @@ namespace Opm const std::string& output_dir); #ifdef HAVE_DUNE_CORNERPOINT void outputStateVtk(const Dune::CpGrid& grid, - const Opm::SimulatorState& state, + const Opm::SimulationDataContainer& state, const int step, const std::string& output_dir); #endif template void outputStateMatlab(const Grid& grid, - const Opm::SimulatorState& state, + const Opm::SimulationDataContainer& state, const int step, const std::string& output_dir) { Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - for( unsigned int i=0; i cell_velocity; @@ -154,7 +156,7 @@ namespace Opm /** \copydoc Opm::OutputWriter::writeTimeStep */ void writeTimeStep(const SimulatorTimerInterface& timer, - const SimulatorState& state, + const SimulationDataContainer& state, const WellState&, bool /*substep*/ = false) { @@ -184,7 +186,7 @@ namespace Opm /** \copydoc Opm::OutputWriter::writeTimeStep */ void writeTimeStep(const SimulatorTimerInterface& timer, - const SimulatorState& reservoirState, + const SimulationDataContainer& reservoirState, const WellState& wellState, bool /*substep*/ = false) { @@ -214,7 +216,7 @@ namespace Opm /** \copydoc Opm::OutputWriter::writeTimeStep */ void writeTimeStep(const SimulatorTimerInterface& timer, - const SimulatorState& reservoirState, + const SimulationDataContainer& reservoirState, const Opm::WellState& wellState, bool substep = false); @@ -235,7 +237,7 @@ namespace Opm void initFromRestartFile(const PhaseUsage& phaseusage, const double* permeability, const Grid& grid, - SimulatorState& simulatorstate, + SimulationDataContainer& simulatorstate, WellStateFullyImplicitBlackoil& wellstate); bool isRestart() const; @@ -315,7 +317,7 @@ namespace Opm initFromRestartFile( const PhaseUsage& phaseusage, const double* permeability, const Grid& grid, - SimulatorState& simulatorstate, + SimulationDataContainer& simulatorstate, WellStateFullyImplicitBlackoil& wellstate) { WellsManager wellsmanager(eclipseState_, diff --git a/opm/autodiff/TransportSolverTwophaseAd.cpp b/opm/autodiff/TransportSolverTwophaseAd.cpp index 65cc55b30..61ba63ab0 100644 --- a/opm/autodiff/TransportSolverTwophaseAd.cpp +++ b/opm/autodiff/TransportSolverTwophaseAd.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include diff --git a/opm/polymer/CompressibleTpfaPolymer.cpp b/opm/polymer/CompressibleTpfaPolymer.cpp index 67feef107..604ba0334 100644 --- a/opm/polymer/CompressibleTpfaPolymer.cpp +++ b/opm/polymer/CompressibleTpfaPolymer.cpp @@ -88,8 +88,8 @@ namespace Opm PolymerBlackoilState& state, WellState& well_state) { - c_ = &state.concentration(); - cmax_ = &state.maxconcentration(); + c_ = &state.getCellData( state.CONCENTRATION ); + cmax_ = &state.getCellData( state.CMAX ); CompressibleTpfa::solve(dt, state, well_state); } diff --git a/opm/polymer/IncompTpfaPolymer.cpp b/opm/polymer/IncompTpfaPolymer.cpp index 138b12499..61905d87b 100644 --- a/opm/polymer/IncompTpfaPolymer.cpp +++ b/opm/polymer/IncompTpfaPolymer.cpp @@ -20,6 +20,8 @@ #include +#include + #include #include @@ -99,8 +101,8 @@ namespace Opm PolymerState& state, WellState& well_state) { - c_ = &state.concentration(); - cmax_ = &state.maxconcentration(); + c_ = &state.getCellData( state.CONCENTRATION ); + cmax_ = &state.getCellData( state.CMAX) ; if (rock_comp_props_ != 0 && rock_comp_props_->isActive()) { solveRockComp(dt, state, well_state); } else { @@ -116,7 +118,7 @@ namespace Opm /// Compute per-solve dynamic properties. void IncompTpfaPolymer::computePerSolveDynamicData(const double /*dt*/, - const SimulatorState& state, + const SimulationDataContainer& state, const WellState& /*well_state*/) { // Computed here: diff --git a/opm/polymer/IncompTpfaPolymer.hpp b/opm/polymer/IncompTpfaPolymer.hpp index 8f61e5dc9..9feeb9bbe 100644 --- a/opm/polymer/IncompTpfaPolymer.hpp +++ b/opm/polymer/IncompTpfaPolymer.hpp @@ -37,6 +37,7 @@ namespace Opm class LinearSolverInterface; class PolymerState; class WellState; + class SimulationDataContainer; /// Encapsulating a tpfa pressure solver for the incompressible-fluid case with polymer. /// Supports gravity, wells controlled by bhp or reservoir rates, @@ -96,7 +97,7 @@ namespace Opm private: virtual void computePerSolveDynamicData(const double dt, - const SimulatorState& state, + const SimulationDataContainer& state, const WellState& well_state); private: // ------ Data that will remain unmodified after construction. ------ diff --git a/opm/polymer/PolymerBlackoilState.cpp b/opm/polymer/PolymerBlackoilState.cpp new file mode 100644 index 000000000..e89d16c29 --- /dev/null +++ b/opm/polymer/PolymerBlackoilState.cpp @@ -0,0 +1,42 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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 + +#include + +namespace Opm +{ + const std::string PolymerBlackoilState::CONCENTRATION = "CONCENTRATION"; + const std::string PolymerBlackoilState::CMAX = "CMAX"; + + PolymerBlackoilState::PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases) : + BlackoilState( number_of_cells , number_of_faces, num_phases) + { + registerCellData(CONCENTRATION , 1 , 0 ); + registerCellData(CMAX , 1 , 0 ); + } + + +} // namespace Opm + + + + diff --git a/opm/polymer/PolymerBlackoilState.hpp b/opm/polymer/PolymerBlackoilState.hpp index 6fb1c36e5..521c284f4 100644 --- a/opm/polymer/PolymerBlackoilState.hpp +++ b/opm/polymer/PolymerBlackoilState.hpp @@ -33,27 +33,10 @@ namespace Opm class PolymerBlackoilState : public BlackoilState { public: - void init(const UnstructuredGrid& g, int num_phases) - { - this->init(g.number_of_cells, g.number_of_faces, num_phases); - } + static const std::string CONCENTRATION; + static const std::string CMAX; - void init(int number_of_cells, int number_of_faces, int num_phases) - { - BlackoilState::init(number_of_cells, number_of_faces, num_phases); - concentration_.resize(number_of_cells, 0.0); - cmax_.resize(number_of_cells, 0.0); - } - - std::vector& concentration() { return concentration_; } - std::vector& maxconcentration() { return cmax_; } - - const std::vector& concentration() const { return concentration_; } - const std::vector& maxconcentration() const { return cmax_; } - - private: - std::vector concentration_; - std::vector cmax_; + PolymerBlackoilState(int number_of_cells, int number_of_faces, int num_phases); }; } // namespace Opm diff --git a/opm/polymer/PolymerState.cpp b/opm/polymer/PolymerState.cpp new file mode 100644 index 000000000..83c323d56 --- /dev/null +++ b/opm/polymer/PolymerState.cpp @@ -0,0 +1,42 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + + 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 + +#include + +namespace Opm +{ + const std::string PolymerState::CONCENTRATION = "CONCENTRATION"; + const std::string PolymerState::CMAX = "CMAX"; + + PolymerState::PolymerState(int number_of_cells, int number_of_faces, int num_phases) : + SimulationDataContainer( number_of_cells , number_of_faces , num_phases ) + { + registerCellData(CONCENTRATION , 1 , 0 ); + registerCellData(CMAX , 1 , 0 ); + } + + +} // namespace Opm + + + + diff --git a/opm/polymer/PolymerState.hpp b/opm/polymer/PolymerState.hpp index cf6cd5fd6..87f810a00 100644 --- a/opm/polymer/PolymerState.hpp +++ b/opm/polymer/PolymerState.hpp @@ -20,8 +20,8 @@ #ifndef OPM_POLYMERSTATE_HEADER_INCLUDED #define OPM_POLYMERSTATE_HEADER_INCLUDED +#include -#include #include #include @@ -29,34 +29,13 @@ namespace Opm { /// Simulator state for a two-phase simulator with polymer. - class PolymerState : public SimulatorState + class PolymerState : public SimulationDataContainer { public: + static const std::string CONCENTRATION; + static const std::string CMAX; - void init(int number_of_cells, int number_of_faces, int num_phases) - { - SimulatorState::init( number_of_cells , number_of_faces , num_phases ); - registerCellData("CONCENTRATION" , 1 , 0 ); - registerCellData("CMAX" , 1 , 0 ); - } - - const std::vector& concentration() const { - return getCellData("CONCENTRATION"); - } - - const std::vector& maxconcentration() const { - return getCellData("CMAX"); - } - - - std::vector& concentration() { - return getCellData("CONCENTRATION"); - } - - std::vector& maxconcentration() { - return getCellData("CMAX"); - } - + PolymerState(int number_of_cells, int number_of_faces, int num_phases); }; } // namespace Opm diff --git a/opm/polymer/SimulatorCompressiblePolymer.cpp b/opm/polymer/SimulatorCompressiblePolymer.cpp index a483cda58..b20e91492 100644 --- a/opm/polymer/SimulatorCompressiblePolymer.cpp +++ b/opm/polymer/SimulatorCompressiblePolymer.cpp @@ -268,7 +268,7 @@ namespace Opm total_timer.start(); double init_surfvol[2] = { 0.0 }; double inplace_surfvol[2] = { 0.0 }; - double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); + double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol()); double polymass_adsorbed = computePolymerAdsorbed(grid_, props_, poly_props_, state, rock_comp_props_); double init_polymass = polymass + polymass_adsorbed; double tot_injected[2] = { 0.0 }; @@ -301,7 +301,7 @@ namespace Opm if (check_well_controls_) { computeFractionalFlow(props_, poly_props_, allcells_, state.pressure(), state.temperature(), state.surfacevol(), state.saturation(), - state.concentration(), state.maxconcentration(), + state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ) , fractional_flows); wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); } @@ -397,7 +397,7 @@ namespace Opm state.pressure(), state.temperature(), &initial_porevol[0], &porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize, state.saturation(), state.surfacevol(), - state.concentration(), state.maxconcentration()); + state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX )); double substep_injected[2] = { 0.0 }; double substep_produced[2] = { 0.0 }; double substep_polyinj = 0.0; @@ -416,7 +416,7 @@ namespace Opm if (gravity_ != 0 && use_segregation_split_) { tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol(), - state.concentration(), state.maxconcentration()); + state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX )); } } transport_timer.stop(); @@ -426,7 +426,7 @@ namespace Opm // Report volume balances. Opm::computeSaturatedVol(porevol, state.surfacevol(), inplace_surfvol); - polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); + polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol()); polymass_adsorbed = Opm::computePolymerAdsorbed(grid_, props_, poly_props_, state, rock_comp_props_); tot_injected[0] += injected[0]; @@ -539,8 +539,8 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - dm["concentration"] = &state.concentration(); - dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ; + dm["cmax"] = &state.getCellData( state.CMAX ) ; dm["surfvol"] = &state.surfacevol(); std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); @@ -556,8 +556,8 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - dm["concentration"] = &state.concentration(); - dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ; + dm["cmax"] = &state.getCellData( state.CMAX ) ; dm["surfvol"] = &state.surfacevol(); std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); diff --git a/opm/polymer/SimulatorPolymer.cpp b/opm/polymer/SimulatorPolymer.cpp index 615bbd4eb..973b61e23 100644 --- a/opm/polymer/SimulatorPolymer.cpp +++ b/opm/polymer/SimulatorPolymer.cpp @@ -292,8 +292,8 @@ namespace Opm total_timer.start(); double init_satvol[2] = { 0.0 }; double satvol[2] = { 0.0 }; - double polymass = computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); - double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration()); + double polymass = computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol()); + double polymass_adsorbed = computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX )); double init_polymass = polymass + polymass_adsorbed; double injected[2] = { 0.0 }; double produced[2] = { 0.0 }; @@ -330,7 +330,7 @@ namespace Opm // Solve pressure. if (check_well_controls_) { computeFractionalFlow(props_, poly_props_, allcells_, - state.saturation(), state.concentration(), state.maxconcentration(), + state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX ), fractional_flows); wells_manager_.applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); } @@ -425,7 +425,7 @@ namespace Opm injected[0] = injected[1] = produced[0] = produced[1] = polyinj = polyprod = 0.0; for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) { tsolver_.solve(&state.faceflux()[0], &initial_porevol[0], &transport_src[0], &polymer_inflow_c[0], stepsize, - state.saturation(), state.concentration(), state.maxconcentration()); + state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX )); Opm::computeInjectedProduced(props_, poly_props_, state, transport_src, polymer_inflow_c, stepsize, @@ -438,7 +438,7 @@ namespace Opm polyprod += substep_polyprod; if (use_segregation_split_) { tsolver_.solveGravity(columns_, &porevol[0], stepsize, - state.saturation(), state.concentration(), state.maxconcentration()); + state.saturation(), state.getCellData( state.CONCENTRATION ), state.getCellData( state.CMAX )); } } transport_timer.stop(); @@ -448,8 +448,8 @@ namespace Opm // Report volume balances. Opm::computeSaturatedVol(porevol, state.saturation(), satvol); - polymass = Opm::computePolymerMass(porevol, state.saturation(), state.concentration(), poly_props_.deadPoreVol()); - polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.maxconcentration()); + polymass = Opm::computePolymerMass(porevol, state.saturation(), state.getCellData( state.CONCENTRATION ), poly_props_.deadPoreVol()); + polymass_adsorbed = Opm::computePolymerAdsorbed(props_, poly_props_, porevol, state.getCellData( state.CMAX )); tot_injected[0] += injected[0]; tot_injected[1] += injected[1]; tot_produced[0] += produced[0]; @@ -559,8 +559,8 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - dm["concentration"] = &state.concentration(); - dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ; + dm["cmax"] = &state.getCellData( state.CMAX ) ; std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); dm["velocity"] = &cell_velocity; @@ -575,8 +575,8 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - dm["concentration"] = &state.concentration(); - dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ; + dm["cmax"] = &state.getCellData( state.CMAX ) ; std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); dm["velocity"] = &cell_velocity; @@ -611,8 +611,8 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - dm["concentration"] = &state.concentration(); - dm["cmax"] = &state.maxconcentration(); + dm["concentration"] = &state.getCellData( state.CONCENTRATION ) ; + dm["cmax"] = &state.getCellData( state.CMAX ) ; std::vector cell_velocity; Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); dm["velocity"] = &cell_velocity; diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 83dead022..485c48e14 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -125,8 +125,10 @@ namespace Opm { WellState& well_state) { Base::prepareStep(dt, reservoir_state, well_state); + auto& max_concentration = reservoir_state.getCellData( reservoir_state.CMAX ); // Initial max concentration of this time step from PolymerBlackoilState. - cmax_ = Eigen::Map(reservoir_state.maxconcentration().data(), Opm::AutoDiffGrid::numCells(grid_)); + + cmax_ = Eigen::Map(max_concentration.data(), Opm::AutoDiffGrid::numCells(grid_)); } @@ -168,9 +170,10 @@ namespace Opm { // Initial polymer concentration. if (has_polymer_) { - assert (not x.concentration().empty()); - const int nc = x.concentration().size(); - const V c = Eigen::Map(&x.concentration()[0], nc); + const auto& concentration = x.getCellData( x.CONCENTRATION ); + assert (not concentration.empty()); + const int nc = concentration.size(); + const V c = Eigen::Map(concentration.data() , nc); // Concentration belongs after other reservoir vars but before well vars. auto concentration_pos = vars0.begin() + fluid_.numPhases(); assert(concentration_pos == vars0.end() - 2); @@ -255,12 +258,14 @@ namespace Opm { template void BlackoilPolymerModel::computeCmax(ReservoirState& state) { - const int nc = AutoDiffGrid::numCells(grid_); - V tmp = V::Zero(nc); - for (int i = 0; i < nc; ++i) { - tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]); - } - std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin()); + auto& max_concentration = state.getCellData( state.CMAX ); + const auto& concentration = state.getCellData( state.CONCENTRATION ); + std::transform( max_concentration.begin() , + max_concentration.end() , + concentration.begin() , + max_concentration.begin() , + [](double c_max , double c) { return std::max( c_max , c ); }); + } @@ -397,10 +402,13 @@ namespace Opm { // Call base version. Base::updateState(modified_dx, reservoir_state, well_state); - // Update concentration. - const V c_old = Eigen::Map(&reservoir_state.concentration()[0], nc, 1); - const V c = (c_old - dc).max(zero); - std::copy(&c[0], &c[0] + nc, reservoir_state.concentration().begin()); + { + auto& concentration = reservoir_state.getCellData( reservoir_state.CONCENTRATION ); + // Update concentration. + const V c_old = Eigen::Map(concentration.data(), nc, 1); + const V c = (c_old - dc).max(zero); + std::copy(&c[0], &c[0] + nc, concentration.begin()); + } } else { // Just forward call to base version. Base::updateState(dx, reservoir_state, well_state); diff --git a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp index 1ca4df5dd..6e9e195e7 100644 --- a/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp +++ b/opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp @@ -206,7 +206,7 @@ namespace { const std::vector& polymer_inflow = xw.polymerInflow(); // Initial max concentration of this time step from PolymerBlackoilState. - cmax_ = Eigen::Map(&x.maxconcentration()[0], Opm::AutoDiffGrid::numCells(grid_)); + cmax_ = Eigen::Map(&x.getCellData( x.CMAX )[0], Opm::AutoDiffGrid::numCells(grid_)); const SolutionState state = constantState(x, xw); computeAccum(state, 0); @@ -366,9 +366,18 @@ namespace { state.saturation[1] = ADB::constant(so, bpat); // Concentration - assert(not x.concentration().empty()); - const V c = Eigen::Map(&x.concentration()[0], nc); - state.concentration = ADB::constant(c); + { + auto& concentration = x.getCellData( x.CONCENTRATION ); + assert(concentration.empty()); + const V c = Eigen::Map(concentration.data(), nc); + // Do not understand: + //concentration = ADB::constant(c); + // Old code based on concentraton() method had the statement: + // + // state.concentration = ADB::constant(c) + // + // This looks like it was a method assignment - how did it even compile? + } // Well rates. assert (not xw.wellRates().empty()); @@ -413,9 +422,12 @@ namespace { vars0.push_back(sw0); // Initial concentration. - assert (not x.concentration().empty()); - const V c = Eigen::Map(&x.concentration()[0], nc); - vars0.push_back(c); + { + auto& concentration = x.getCellData( x.CONCENTRATION ); + assert (not concentration.empty()); + const V c = Eigen::Map(concentration.data() , nc); + vars0.push_back(c); + } // Initial well rates. assert (not xw.wellRates().empty()); @@ -511,11 +523,14 @@ namespace { { const int nc = grid_.number_of_cells; V tmp = V::Zero(nc); + const auto& concentration = state.getCellData( state.CONCENTRATION ); + auto& cmax = state.getCellData( state.CMAX ); + for (int i = 0; i < nc; ++i) { - tmp[i] = std::max(state.maxconcentration()[i], state.concentration()[i]); + tmp[i] = std::max(cmax[i], concentration[i]); } - std::copy(&tmp[0], &tmp[0] + nc, state.maxconcentration().begin()); + std::copy(&tmp[0], &tmp[0] + nc, cmax.begin()); } @@ -764,11 +779,14 @@ namespace { // Concentration updates. // const double dcmax = 0.3 * polymer_props_ad_.cMax(); // std::cout << "\n the max concentration: " << dcmax / 0.3 << std::endl; - const V c_old = Eigen::Map(&state.concentration()[0], nc, 1); -// const V dc_limited = sign(dc) * dc.abs().min(dcmax); -// const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02()); - const V c = (c_old - dc).max(zero); - std::copy(&c[0], &c[0] + nc, state.concentration().begin()); + { + auto& concentration = state.getCellData( state.CONCENTRATION ); + const V c_old = Eigen::Map(concentration.data() , nc, 1); + // const V dc_limited = sign(dc) * dc.abs().min(dcmax); + // const V c = (c_old - dc_limited).max(zero);//unaryExpr(Chop02()); + const V c = (c_old - dc).max(zero); + std::copy(&c[0], &c[0] + nc, concentration.begin()); + } // Qs update. // Since we need to update the wellrates, that are ordered by wells, diff --git a/opm/polymer/polymerUtilities.cpp b/opm/polymer/polymerUtilities.cpp index 7fa369fa3..971066f5f 100644 --- a/opm/polymer/polymerUtilities.cpp +++ b/opm/polymer/polymerUtilities.cpp @@ -211,8 +211,8 @@ namespace Opm OPM_THROW(std::runtime_error, "Sizes of state vectors do not match number of cells."); } const std::vector& s = state.saturation(); - const std::vector& c = state.concentration(); - const std::vector& cmax = state.maxconcentration(); + const std::vector& c = state.getCellData( state.CONCENTRATION ); + const std::vector& cmax = state.getCellData( state.CMAX ); std::fill(injected, injected + np, 0.0); std::fill(produced, produced + np, 0.0); polyinj = 0.0; @@ -282,8 +282,8 @@ namespace Opm const std::vector& temp = state.temperature(); const std::vector& s = state.saturation(); const std::vector& z = state.surfacevol(); - const std::vector& c = state.concentration(); - const std::vector& cmax = state.maxconcentration(); + const std::vector& c = state.getCellData( state.CONCENTRATION ); + const std::vector& cmax = state.getCellData( state.CMAX ); std::fill(injected, injected + np, 0.0); std::fill(produced, produced + np, 0.0); polyinj = 0.0; @@ -406,7 +406,7 @@ namespace Opm porosity.assign(props.porosity(), props.porosity() + num_cells); } double abs_mass = 0.0; - const std::vector& cmax = state.maxconcentration(); + const std::vector& cmax = state.getCellData( state.CMAX ); for (int cell = 0; cell < num_cells; ++cell) { double c_ads; polyprops.simpleAdsorption(cmax[cell], c_ads); diff --git a/tests/test_rateconverter.cpp b/tests/test_rateconverter.cpp index 6c9e28898..15ec67aee 100644 --- a/tests/test_rateconverter.cpp +++ b/tests/test_rateconverter.cpp @@ -32,6 +32,7 @@ #include +#include #include #include #include @@ -107,8 +108,7 @@ BOOST_FIXTURE_TEST_CASE(ThreePhase, TestFixture) Region reg{ 0 }; RCvrt cvrt(ad_props, reg); - Opm::BlackoilState x; - x.init(*grid.c_grid(), 3); + Opm::BlackoilState x( Opm::UgGridHelpers::numCells( *grid.c_grid()) , Opm::UgGridHelpers::numFaces( *grid.c_grid()) , 3); cvrt.defineState(x); diff --git a/tests/test_singlecellsolves.cpp b/tests/test_singlecellsolves.cpp index 25f298d85..1a9ba9800 100644 --- a/tests/test_singlecellsolves.cpp +++ b/tests/test_singlecellsolves.cpp @@ -71,7 +71,7 @@ try boost::scoped_ptr grid; boost::scoped_ptr props; - PolymerState state; + std::unique_ptr state; Opm::PolymerProperties poly_props; // bool check_well_controls = false; // int max_well_control_iterations = 0; @@ -81,23 +81,35 @@ try // Grid init. grid.reset(new GridManager(2, 1, 1, 1.0, 1.0, 1.0)); // Rock and fluid init. - props.reset(new IncompPropertiesBasic(param, grid->c_grid()->dimensions, grid->c_grid()->number_of_cells)); - // Init state variables (saturation and pressure). - initStateBasic(*grid->c_grid(), *props, param, 0.0, state); - // Init Polymer state - if (param.has("poly_init")) { - double poly_init = param.getDefault("poly_init", 0.0); - for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) { - double smin[2], smax[2]; - props->satRange(1, &cell, smin, smax); - if (state.saturation()[2*cell] > 0.5*(smin[0] + smax[0])) { - state.concentration()[cell] = poly_init; - state.maxconcentration()[cell] = poly_init; - } else { - state.saturation()[2*cell + 0] = 0.; - state.saturation()[2*cell + 1] = 1.; - state.concentration()[cell] = 0.; - state.maxconcentration()[cell] = 0.; + { + const UnstructuredGrid& ug_grid = *(grid->c_grid()); + state.reset( new PolymerState( UgGridHelpers::numCells( ug_grid ) , UgGridHelpers::numFaces( ug_grid ), 2)); + + props.reset(new IncompPropertiesBasic(param, ug_grid.dimensions, UgGridHelpers::numCells( ug_grid ))); + // Init state variables (saturation and pressure). + initStateBasic(*grid->c_grid(), *props, param, 0.0, *state); + // Init Polymer state + if (param.has("poly_init")) { + double poly_init = param.getDefault("poly_init", 0.0); + for (int cell = 0; cell < grid->c_grid()->number_of_cells; ++cell) { + double smin[2], smax[2]; + props->satRange(1, &cell, smin, smax); + + auto& saturation = state->saturation(); + auto& concentration = state->getCellData( state->CONCENTRATION ); + auto& max_concentration = state->getCellData( state->CMAX ); + + props->satRange(1, &cell, smin, smax); + if (saturation[2*cell] > 0.5*(smin[0] + smax[0])) { + concentration[cell] = poly_init; + max_concentration[cell] = poly_init; + } else { + saturation[2*cell + 0] = 0.; + saturation[2*cell + 1] = 1.; + concentration[cell] = 0.; + max_concentration[cell] = 0.; + } + } } } @@ -210,7 +222,7 @@ try if (face01 == -1) { OPM_THROW(std::runtime_error, "Could not find face adjacent to cells [0 1]"); } - state.faceflux()[face01] = src[0]; + state->faceflux()[face01] = src[0]; for (int sats = 0; sats < num_sats; ++sats) { const double s = double(sats)/double(num_sats - 1); const double ff = s; // Simplified a lot... @@ -220,20 +232,26 @@ try // std::cout << "(s, c) = (" << s << ", " << c << ")\n"; transport_src[0] = src[0]*ff; // Resetting the state for next run. - state.saturation()[0] = 0.0; - state.saturation()[1] = 0.0; - state.concentration()[0] = 0.0; - state.concentration()[1] = 0.0; - state.maxconcentration()[0] = 0.0; - state.maxconcentration()[1] = 0.0; - reorder_model.solve(&state.faceflux()[0], + auto& saturation = state->saturation(); + auto& concentration = state->getCellData( state->CONCENTRATION ); + auto& max_concentration = state->getCellData( state->CMAX ); + + + saturation[0] = 0.0; + saturation[1] = 0.0; + concentration[0] = 0.0; + concentration[1] = 0.0; + max_concentration[0] = 0.0; + max_concentration[1] = 0.0; + + reorder_model.solve(&state->faceflux()[0], &porevol[0], &transport_src[0], &polymer_inflow_c[0], dt, - state.saturation(), - state.concentration(), - state.maxconcentration()); + saturation, + concentration, + max_concentration); #ifdef PROFILING // Extract residual counts. From 98862f449c32650b071a9f2462c374cf94527901 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Mar 2016 10:40:48 +0100 Subject: [PATCH 16/24] Restore compilation (init -> constructor). Also suppress new warnings due to unsigned/signed comparisons. --- opm/autodiff/RedistributeDataHandles.hpp | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/opm/autodiff/RedistributeDataHandles.hpp b/opm/autodiff/RedistributeDataHandles.hpp index e401f00d0..e2c5724c9 100644 --- a/opm/autodiff/RedistributeDataHandles.hpp +++ b/opm/autodiff/RedistributeDataHandles.hpp @@ -232,7 +232,7 @@ public: { assert( T::codimension == 0); - for ( int i=0; i(size_arg); double val; - for ( int i=0; i::max()); BlackoilStateDataHandle state_handle(global_grid, grid, From 463866b523e1b14a640c61fe37bf912fa1b8aff2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 15 Mar 2016 14:03:55 +0100 Subject: [PATCH 17/24] Create state_ in all cases. --- opm/autodiff/FlowMain.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index 1f5c7f4f0..b8af7d6a4 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -483,6 +483,9 @@ namespace Opm initStateEquil(grid, props, deck_, eclipse_state_, gravity_[2], *state_); //state_.faceflux().resize(Opm::UgGridHelpers::numFaces(grid), 0.0); } else { + state_.reset( new ReservoirState( Opm::UgGridHelpers::numCells(grid), + Opm::UgGridHelpers::numFaces(grid), + props.numPhases())); initBlackoilStateFromDeck(Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid), Opm::UgGridHelpers::numFaces(grid), From 29053dc18de06db8d53e658fed66dbb1afa6d210 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Mon, 21 Mar 2016 10:58:50 +0100 Subject: [PATCH 18/24] SimulFulImplBOOutput: asynchronous serial output. --- .../SimulatorFullyImplicitBlackoilOutput.cpp | 146 ++++++++++++------ .../SimulatorFullyImplicitBlackoilOutput.hpp | 6 + 2 files changed, 101 insertions(+), 51 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index a57fcdb8a..f57eb5872 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -1,6 +1,6 @@ /* Copyright (c) 2014 SINTEF ICT, Applied Mathematics. - Copyright (c) 2015 IRIS AS + Copyright (c) 2015-2016 IRIS AS This file is part of the Open Porous Media project (OPM). @@ -19,6 +19,8 @@ */ #include "config.h" +#include + #include "SimulatorFullyImplicitBlackoilOutput.hpp" #include @@ -252,6 +254,34 @@ namespace Opm } } + struct WriterCall + { + BlackoilOutputWriter& writer_; + std::unique_ptr< SimulatorTimerInterface > timer_; + const SimulatorState state_; + const WellState wellState_; + const bool substep_; + + explicit WriterCall( BlackoilOutputWriter& writer, + const SimulatorTimerInterface& timer, + const SimulatorState& state, + const WellState& wellState, + bool substep) + : writer_( writer ), + timer_( timer.clone() ), + state_( state ), + wellState_( wellState ), + substep_( substep ) + {} + + // callback to writer's serial writeTimeStep method + void operator () () + { + writer_.writeTimeStepSerial( *timer_, state_, wellState_, substep_ ); + } + }; + + void BlackoilOutputWriter:: writeTimeStep(const SimulatorTimerInterface& timer, @@ -274,63 +304,77 @@ namespace Opm const SimulatorState& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState; const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState; - // output is only done on I/O rank + // serial output is only done on I/O rank if( isIORank ) { - // Matlab output - if( matlabWriter_ ) { - matlabWriter_->writeTimeStep( timer, state, wellState, substep ); - } - // ECL output - if ( eclWriter_ ) { - const auto initConfig = eclipseState_->getInitConfig(); - if (initConfig->getRestartInitiated() && ((initConfig->getRestartStep()) == (timer.currentStepNum()))) { - std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; - } else { - eclWriter_->writeTimeStep(timer, state, wellState, substep ); - } - } + // spawn write thread that calls eclWriter.writeTimeStepSerial + WriterCall call( *this, timer, state, wellState, substep ); + std::async( std::move( call ) ); + } + } - // write backup file - if( backupfile_ ) + void + BlackoilOutputWriter:: + writeTimeStepSerial(const SimulatorTimerInterface& timer, + const SimulatorState& state, + const WellState& wellState, + bool substep) + { + // Matlab output + if( matlabWriter_ ) { + matlabWriter_->writeTimeStep( timer, state, wellState, substep ); + } + + // ECL output + if ( eclWriter_ ) + { + const auto initConfig = eclipseState_->getInitConfig(); + if (initConfig->getRestartInitiated() && ((initConfig->getRestartStep()) == (timer.currentStepNum()))) { + std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; + } else { + eclWriter_->writeTimeStep(timer, state, wellState, substep ); + } + } + + // write backup file + if( backupfile_ ) + { + int reportStep = timer.reportStepNum(); + int currentTimeStep = timer.currentStepNum(); + if( (reportStep == currentTimeStep || // true for SimulatorTimer + currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep + timer.done() ) // true for AdaptiveSimulatorTimer at reportStep + && lastBackupReportStep_ != reportStep ) // only backup report step once { - int reportStep = timer.reportStepNum(); - int currentTimeStep = timer.currentStepNum(); - if( (reportStep == currentTimeStep || // true for SimulatorTimer - currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep - timer.done() ) // true for AdaptiveSimulatorTimer at reportStep - && lastBackupReportStep_ != reportStep ) // only backup report step once - { - // store report step - lastBackupReportStep_ = reportStep; - // write resport step number - backupfile_.write( (const char *) &reportStep, sizeof(int) ); + // store report step + lastBackupReportStep_ = reportStep; + // write resport step number + backupfile_.write( (const char *) &reportStep, sizeof(int) ); - try { - const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state); - backupfile_ << boState; + try { + const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state); + backupfile_ << boState; - const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState); - backupfile_ << boWellState; - } - catch ( const std::bad_cast& e ) - { - - } - - /* - const WellStateFullyImplicitBlackoil* boWellState = - dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState); - if( boWellState ) { - backupfile_ << (*boWellState); - } - else - OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed"); - */ - backupfile_ << std::flush; + const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState); + backupfile_ << boWellState; } - } // end backup - } // end isIORank + catch ( const std::bad_cast& e ) + { + + } + + /* + const WellStateFullyImplicitBlackoil* boWellState = + dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState); + if( boWellState ) { + backupfile_ << (*boWellState); + } + else + OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed"); + */ + backupfile_ << std::flush; + } + } // end backup } void diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index a4e8e879c..1cd42e9e4 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -218,6 +218,12 @@ namespace Opm const Opm::WellState& wellState, bool substep = false); + /** \copydoc Opm::OutputWriter::writeTimeStep */ + void writeTimeStepSerial(const SimulatorTimerInterface& timer, + const SimulatorState& reservoirState, + const Opm::WellState& wellState, + bool substep); + /** \brief return output directory */ const std::string& outputDirectory() const { return outputDir_; } From d5c375f297806b948f5f278664b355c1740f7b35 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Tue, 29 Mar 2016 10:48:24 +0200 Subject: [PATCH 19/24] SimulatorFullyImplicitOutput: added flag for async output. --- .../SimulatorFullyImplicitBlackoilOutput.cpp | 12 +++++++++--- .../SimulatorFullyImplicitBlackoilOutput.hpp | 4 +++- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index f57eb5872..372ed44b9 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -307,9 +307,15 @@ namespace Opm // serial output is only done on I/O rank if( isIORank ) { - // spawn write thread that calls eclWriter.writeTimeStepSerial - WriterCall call( *this, timer, state, wellState, substep ); - std::async( std::move( call ) ); + if( asyncOutput_ ) { + // spawn write thread that calls eclWriter.writeTimeStepSerial + WriterCall call( *this, timer, state, wellState, substep ); + std::async( std::move( call ) ); + } + else { + // just write the data to disk + writeTimeStepSerial( timer, state, wellState, substep ); + } } } diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index 1cd42e9e4..ad02b75c4 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -261,6 +261,7 @@ namespace Opm std::unique_ptr< OutputWriter > matlabWriter_; std::unique_ptr< EclipseWriter > eclWriter_; EclipseStateConstPtr eclipseState_; + const bool asyncOutput_ ; }; @@ -293,7 +294,8 @@ namespace Opm parallelOutput_->numCells(), parallelOutput_->globalCell() ) : 0 ), - eclipseState_(eclipseState) + eclipseState_(eclipseState), + asyncOutput_( output_ ? param.getDefault("async_output", bool( false ) ) : false ) { // For output. if (output_ && parallelOutput_->isIORank() ) { From 5044a07e467cc6e5014918d27a589862419cb054 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Tue, 29 Mar 2016 15:22:33 +0200 Subject: [PATCH 20/24] SimFullImplOutput: added future to ensure that write history is accounted for. --- .../SimulatorFullyImplicitBlackoilOutput.cpp | 25 ++++++++++++++----- .../SimulatorFullyImplicitBlackoilOutput.hpp | 2 ++ 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index 372ed44b9..4a393fc28 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -262,22 +262,33 @@ namespace Opm const WellState wellState_; const bool substep_; + std::future< bool > asyncWait_; + explicit WriterCall( BlackoilOutputWriter& writer, const SimulatorTimerInterface& timer, const SimulatorState& state, const WellState& wellState, - bool substep) + bool substep, + std::future< bool >&& asyncWait) : writer_( writer ), timer_( timer.clone() ), state_( state ), wellState_( wellState ), - substep_( substep ) + substep_( substep ), + asyncWait_( std::move( asyncWait ) ) {} // callback to writer's serial writeTimeStep method - void operator () () + bool operator () () { + // wait for previous time step thread to be finished + // to ensure that write history is correct. + if( asyncWait_.valid() ) { + asyncWait_.wait(); + } + writer_.writeTimeStepSerial( *timer_, state_, wellState_, substep_ ); + return true; } }; @@ -294,7 +305,7 @@ namespace Opm vtkWriter_->writeTimeStep( timer, localState, localWellState, false ); } - bool isIORank = true ; + bool isIORank = output_ ; if( parallelOutput_ && parallelOutput_->isParallel() ) { // collect all solutions to I/O rank @@ -308,9 +319,11 @@ namespace Opm if( isIORank ) { if( asyncOutput_ ) { + // spawn write thread that calls eclWriter.writeTimeStepSerial - WriterCall call( *this, timer, state, wellState, substep ); - std::async( std::move( call ) ); + // timer, state, and wellState are copied, previous async future is moved + WriterCall call( *this, timer, state, wellState, substep, std::move( asyncWait_ ) ); + asyncWait_ = std::move( std::async( std::move( call ) ) ); } else { // just write the data to disk diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index ad02b75c4..e0208bff2 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -44,6 +44,7 @@ #include #include #include +#include #include @@ -262,6 +263,7 @@ namespace Opm std::unique_ptr< EclipseWriter > eclWriter_; EclipseStateConstPtr eclipseState_; const bool asyncOutput_ ; + std::future< bool > asyncWait_; }; From 9df3b2fda974a3d7bf3c3ef787e73a62df8a75b3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Kvalsvik?= Date: Wed, 16 Mar 2016 14:40:54 +0100 Subject: [PATCH 21/24] Flow accepts base name for input Deck Enables flow to accept a basename for a case by appending a .DATA suffix should it not be provided. It already supported reading the basename from a .DATA extension file, but not opening said file by handing it to the parser. --- opm/autodiff/FlowMain.hpp | 35 ++++++++++++++++++++++++++++++----- 1 file changed, 30 insertions(+), 5 deletions(-) diff --git a/opm/autodiff/FlowMain.hpp b/opm/autodiff/FlowMain.hpp index b3095fd2b..1c39c5c28 100644 --- a/opm/autodiff/FlowMain.hpp +++ b/opm/autodiff/FlowMain.hpp @@ -94,6 +94,7 @@ #include #include #include +#include @@ -101,6 +102,33 @@ namespace Opm { + boost::filesystem::path simulationCaseName( const std::string& casename ) { + namespace fs = boost::filesystem; + + const auto exists = []( const fs::path& f ) -> bool { + if( !fs::exists( f ) ) return false; + + if( fs::is_regular_file( f ) ) return true; + + return fs::is_symlink( f ) + && fs::is_regular_file( fs::read_symlink( f ) ); + }; + + auto simcase = fs::path( casename ); + + if( exists( simcase ) ) { + return simcase; + } + + for( const auto& ext : { std::string("data"), std::string("DATA") } ) { + if( exists( simcase.replace_extension( ext ) ) ) { + return simcase; + } + } + + throw std::invalid_argument( "Cannot find input case " + casename ); + } + /// This class encapsulates the setup and running of /// a simulator based on an input deck. template @@ -258,10 +286,6 @@ namespace Opm } } - - - - // Read parameters, see if a deck was specified on the command line, and if // it was, insert it into parameters. // Writes to: @@ -281,7 +305,8 @@ namespace Opm std::cerr << "You can only specify a single input deck on the command line.\n"; return false; } else { - param_.insertParameter("deck_filename", param_.unhandledArguments()[0]); + const auto casename = simulationCaseName( param_.unhandledArguments()[ 0 ] ); + param_.insertParameter("deck_filename", casename.string() ); } } From 8d11be71e08b201d1652165e6f198c5a0ccd6553 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Wed, 30 Mar 2016 11:01:34 +0200 Subject: [PATCH 22/24] SimulatorState --> SimulationDataContainer. --- .../SimulatorFullyImplicitBlackoilOutput.cpp | 47 ++----------------- .../SimulatorFullyImplicitBlackoilOutput.hpp | 4 +- 2 files changed, 6 insertions(+), 45 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index f0055b268..d1f8e4034 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -259,7 +259,7 @@ namespace Opm { BlackoilOutputWriter& writer_; std::unique_ptr< SimulatorTimerInterface > timer_; - const SimulatorState state_; + const SimulationDataContainer state_; const WellState wellState_; const bool substep_; @@ -267,7 +267,7 @@ namespace Opm explicit WriterCall( BlackoilOutputWriter& writer, const SimulatorTimerInterface& timer, - const SimulatorState& state, + const SimulationDataContainer& state, const WellState& wellState, bool substep, std::future< bool >&& asyncWait) @@ -336,7 +336,7 @@ namespace Opm void BlackoilOutputWriter:: writeTimeStepSerial(const SimulatorTimerInterface& timer, - const SimulatorState& state, + const SimulationDataContainer& state, const WellState& wellState, bool substep) { @@ -372,54 +372,15 @@ namespace Opm backupfile_.write( (const char *) &reportStep, sizeof(int) ); try { - const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state); - backupfile_ << boState; + backupfile_ << state; const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState); backupfile_ << boWellState; } catch ( const std::bad_cast& e ) { - // store report step - lastBackupReportStep_ = reportStep; - // write resport step number - backupfile_.write( (const char *) &reportStep, sizeof(int) ); - - /* - try { - const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state); - backupfile_ << boState; - - const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState); - backupfile_ << boWellState; - } - catch ( const std::bad_cast& e ) - { - - } - */ - - /* - const WellStateFullyImplicitBlackoil* boWellState = - dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState); - if( boWellState ) { - backupfile_ << (*boWellState); - } - else - OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed"); - */ - backupfile_ << std::flush; } - /* - const WellStateFullyImplicitBlackoil* boWellState = - dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState); - if( boWellState ) { - backupfile_ << (*boWellState); - } - else - OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed"); - */ backupfile_ << std::flush; } } // end backup diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index e1ed6915c..c5a0bd908 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -86,7 +86,7 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - for (const auto& pair : state.cellData()) + for (const auto& pair : state.cellData()) { const std::string& name = pair.first; std::string key; @@ -223,7 +223,7 @@ namespace Opm /** \copydoc Opm::OutputWriter::writeTimeStep */ void writeTimeStepSerial(const SimulatorTimerInterface& timer, - const SimulatorState& reservoirState, + const SimulationDataContainer& reservoirState, const Opm::WellState& wellState, bool substep); From b39b9a85a2854b64f4ed34b9b17e0edab3761ff5 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 30 Mar 2016 15:41:03 +0200 Subject: [PATCH 23/24] making the interleaved solver works for blackoil polymer simulator. The CPR solver does not work yet. An error will be thrown if people specify to use CPR linear solver. --- opm/autodiff/FlowMainPolymer.hpp | 32 +++++++++++++++++-- .../fullyimplicit/BlackoilPolymerModel.hpp | 2 ++ .../BlackoilPolymerModel_impl.hpp | 21 ++++++++++++ 3 files changed, 53 insertions(+), 2 deletions(-) diff --git a/opm/autodiff/FlowMainPolymer.hpp b/opm/autodiff/FlowMainPolymer.hpp index b594b3c0b..1a5e647f8 100644 --- a/opm/autodiff/FlowMainPolymer.hpp +++ b/opm/autodiff/FlowMainPolymer.hpp @@ -38,6 +38,10 @@ namespace Opm { protected: using Base = FlowMainBase, Grid, Simulator>; + using Base::eclipse_state_; + using Base::param_; + using Base::fis_solver_; + using Base::parallel_information_; friend Base; // Set in setupGridAndProps() @@ -90,10 +94,34 @@ namespace Opm // Setup linear solver. // Writes to: // fis_solver_ + // Currently, the CPR solver is not ready for polymer solver yet void setupLinearSolver() { - OPM_MESSAGE("Caution: polymer solver currently only works with direct linear solver."); - Base::fis_solver_.reset(new NewtonIterationBlackoilSimple(Base::param_, Base::parallel_information_)); + const std::string cprSolver = "cpr"; + const std::string interleavedSolver = "interleaved"; + const std::string directSolver = "direct"; + const std::string flowDefaultSolver = interleavedSolver; + + std::shared_ptr simCfg = eclipse_state_->getSimulationConfig(); + std::string solver_approach = flowDefaultSolver; + + if (param_.has("solver_approach")) { + solver_approach = param_.template get("solver_approach"); + } else { + if (simCfg->useCPR()) { + solver_approach = cprSolver; + } + } + + if (solver_approach == cprSolver) { + OPM_THROW( std::runtime_error , "CPR solver is not ready for use with polymer solver yet."); + } else if (solver_approach == interleavedSolver) { + fis_solver_.reset(new NewtonIterationBlackoilInterleaved(param_, parallel_information_)); + } else if (solver_approach == directSolver) { + fis_solver_.reset(new NewtonIterationBlackoilSimple(param_, parallel_information_)); + } else { + OPM_THROW( std::runtime_error , "Internal error - solver approach " << solver_approach << " not recognized."); + } } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp index 53d452e21..da06207cf 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp @@ -234,6 +234,8 @@ namespace Opm { const SolutionState& state, WellState& xw); + void updateEquationsScaling(); + void computeMassFlux(const int actph , const V& transi, diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 83dead022..fe2a98056 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -107,6 +107,7 @@ namespace Opm { if (!active_[Water]) { OPM_THROW(std::logic_error, "Polymer must solved in water!\n"); } + residual_.matbalscale.resize(fluid_.numPhases() + 1, 1.1169); // use the same as the water phase // If deck has polymer, residual_ should contain polymer equation. rq_.resize(fluid_.numPhases() + 1); residual_.material_balance_eq.resize(fluid_.numPhases() + 1, ADB::null()); @@ -339,6 +340,26 @@ namespace Opm { residual_.material_balance_eq[poly_pos_] = pvdt_ * (rq_[poly_pos_].accum[1] - rq_[poly_pos_].accum[0]) + ops_.div*rq_[poly_pos_].mflux; } + + + if (param_.update_equations_scaling_) { + updateEquationsScaling(); + } + + } + + + + + + template + void BlackoilPolymerModel::updateEquationsScaling() + { + Base::updateEquationsScaling(); + if (has_polymer_) { + const int water_pos = fluid_.phaseUsage().phase_pos[Water]; + residual_.matbalscale[poly_pos_] = residual_.matbalscale[water_pos]; + } } From 13b86a782511c13cda3e2195b80f45fa481cdafe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 31 Mar 2016 09:43:01 +0200 Subject: [PATCH 24/24] Revert "Asynchronous output." --- .../SimulatorFullyImplicitBlackoilOutput.cpp | 156 ++++++------------ .../SimulatorFullyImplicitBlackoilOutput.hpp | 14 +- 2 files changed, 55 insertions(+), 115 deletions(-) diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp index d1f8e4034..d5f32b21c 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp @@ -1,6 +1,6 @@ /* Copyright (c) 2014 SINTEF ICT, Applied Mathematics. - Copyright (c) 2015-2016 IRIS AS + Copyright (c) 2015 IRIS AS This file is part of the Open Porous Media project (OPM). @@ -19,8 +19,6 @@ */ #include "config.h" -#include - #include "SimulatorFullyImplicitBlackoilOutput.hpp" #include @@ -255,45 +253,6 @@ namespace Opm } } - struct WriterCall - { - BlackoilOutputWriter& writer_; - std::unique_ptr< SimulatorTimerInterface > timer_; - const SimulationDataContainer state_; - const WellState wellState_; - const bool substep_; - - std::future< bool > asyncWait_; - - explicit WriterCall( BlackoilOutputWriter& writer, - const SimulatorTimerInterface& timer, - const SimulationDataContainer& state, - const WellState& wellState, - bool substep, - std::future< bool >&& asyncWait) - : writer_( writer ), - timer_( timer.clone() ), - state_( state ), - wellState_( wellState ), - substep_( substep ), - asyncWait_( std::move( asyncWait ) ) - {} - - // callback to writer's serial writeTimeStep method - bool operator () () - { - // wait for previous time step thread to be finished - // to ensure that write history is correct. - if( asyncWait_.valid() ) { - asyncWait_.wait(); - } - - writer_.writeTimeStepSerial( *timer_, state_, wellState_, substep_ ); - return true; - } - }; - - void BlackoilOutputWriter:: writeTimeStep(const SimulatorTimerInterface& timer, @@ -306,7 +265,7 @@ namespace Opm vtkWriter_->writeTimeStep( timer, localState, localWellState, false ); } - bool isIORank = output_ ; + bool isIORank = true ; if( parallelOutput_ && parallelOutput_->isParallel() ) { // collect all solutions to I/O rank @@ -316,74 +275,65 @@ namespace Opm const SimulationDataContainer& state = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalReservoirState() : localState; const WellState& wellState = (parallelOutput_ && parallelOutput_->isParallel() ) ? parallelOutput_->globalWellState() : localWellState; - // serial output is only done on I/O rank + // output is only done on I/O rank if( isIORank ) { - if( asyncOutput_ ) { - - // spawn write thread that calls eclWriter.writeTimeStepSerial - // timer, state, and wellState are copied, previous async future is moved - WriterCall call( *this, timer, state, wellState, substep, std::move( asyncWait_ ) ); - asyncWait_ = std::move( std::async( std::move( call ) ) ); + // Matlab output + if( matlabWriter_ ) { + matlabWriter_->writeTimeStep( timer, state, wellState, substep ); } - else { - // just write the data to disk - writeTimeStepSerial( timer, state, wellState, substep ); + // ECL output + if ( eclWriter_ ) { + const auto initConfig = eclipseState_->getInitConfig(); + if (initConfig->getRestartInitiated() && ((initConfig->getRestartStep()) == (timer.currentStepNum()))) { + std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; + } else { + eclWriter_->writeTimeStep(timer, state, wellState, substep ); + } } - } - } - void - BlackoilOutputWriter:: - writeTimeStepSerial(const SimulatorTimerInterface& timer, - const SimulationDataContainer& state, - const WellState& wellState, - bool substep) - { - // Matlab output - if( matlabWriter_ ) { - matlabWriter_->writeTimeStep( timer, state, wellState, substep ); - } - - // ECL output - if ( eclWriter_ ) - { - const auto initConfig = eclipseState_->getInitConfig(); - if (initConfig->getRestartInitiated() && ((initConfig->getRestartStep()) == (timer.currentStepNum()))) { - std::cout << "Skipping restart write in start of step " << timer.currentStepNum() << std::endl; - } else { - eclWriter_->writeTimeStep(timer, state, wellState, substep ); - } - } - - // write backup file - if( backupfile_ ) - { - int reportStep = timer.reportStepNum(); - int currentTimeStep = timer.currentStepNum(); - if( (reportStep == currentTimeStep || // true for SimulatorTimer - currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep - timer.done() ) // true for AdaptiveSimulatorTimer at reportStep - && lastBackupReportStep_ != reportStep ) // only backup report step once + // write backup file + if( backupfile_ ) { - // store report step - lastBackupReportStep_ = reportStep; - // write resport step number - backupfile_.write( (const char *) &reportStep, sizeof(int) ); - - try { - backupfile_ << state; - - const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState); - backupfile_ << boWellState; - } - catch ( const std::bad_cast& e ) + int reportStep = timer.reportStepNum(); + int currentTimeStep = timer.currentStepNum(); + if( (reportStep == currentTimeStep || // true for SimulatorTimer + currentTimeStep == 0 || // true for AdaptiveSimulatorTimer at reportStep + timer.done() ) // true for AdaptiveSimulatorTimer at reportStep + && lastBackupReportStep_ != reportStep ) // only backup report step once { - } + // store report step + lastBackupReportStep_ = reportStep; + // write resport step number + backupfile_.write( (const char *) &reportStep, sizeof(int) ); - backupfile_ << std::flush; - } - } // end backup + /* + try { + const BlackoilState& boState = dynamic_cast< const BlackoilState& > (state); + backupfile_ << boState; + + const WellStateFullyImplicitBlackoil& boWellState = static_cast< const WellStateFullyImplicitBlackoil& > (wellState); + backupfile_ << boWellState; + } + catch ( const std::bad_cast& e ) + { + + } + */ + + /* + const WellStateFullyImplicitBlackoil* boWellState = + dynamic_cast< const WellStateFullyImplicitBlackoil* > (&wellState); + if( boWellState ) { + backupfile_ << (*boWellState); + } + else + OPM_THROW(std::logic_error,"cast to WellStateFullyImplicitBlackoil failed"); + */ + backupfile_ << std::flush; + } + } // end backup + } // end isIORank } void diff --git a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp index c5a0bd908..4e51690e6 100644 --- a/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp +++ b/opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp @@ -43,7 +43,6 @@ #include #include #include -#include #include @@ -86,7 +85,7 @@ namespace Opm Opm::DataMap dm; dm["saturation"] = &state.saturation(); dm["pressure"] = &state.pressure(); - for (const auto& pair : state.cellData()) + for (const auto& pair : state.cellData()) { const std::string& name = pair.first; std::string key; @@ -221,12 +220,6 @@ namespace Opm const Opm::WellState& wellState, bool substep = false); - /** \copydoc Opm::OutputWriter::writeTimeStep */ - void writeTimeStepSerial(const SimulatorTimerInterface& timer, - const SimulationDataContainer& reservoirState, - const Opm::WellState& wellState, - bool substep); - /** \brief return output directory */ const std::string& outputDirectory() const { return outputDir_; } @@ -264,8 +257,6 @@ namespace Opm std::unique_ptr< OutputWriter > matlabWriter_; std::unique_ptr< EclipseWriter > eclWriter_; EclipseStateConstPtr eclipseState_; - const bool asyncOutput_ ; - std::future< bool > asyncWait_; }; @@ -298,8 +289,7 @@ namespace Opm parallelOutput_->numCells(), parallelOutput_->globalCell() ) : 0 ), - eclipseState_(eclipseState), - asyncOutput_( output_ ? param.getDefault("async_output", bool( false ) ) : false ) + eclipseState_(eclipseState) { // For output. if (output_ && parallelOutput_->isIORank() ) {