From d45543b8fb27e3731e2eec99e2e6f2a2a150ffe1 Mon Sep 17 00:00:00 2001 From: babrodtk Date: Fri, 7 Aug 2015 15:27:35 +0200 Subject: [PATCH] Proper integration of derivatives for THP --- opm/autodiff/BlackoilModelBase_impl.hpp | 56 +++-- opm/autodiff/VFPProperties.cpp | 271 ++++++++++++++++++++---- opm/autodiff/VFPProperties.hpp | 11 +- tests/test_vfpproperties.cpp | 19 +- 4 files changed, 291 insertions(+), 66 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 2bc709419..6d677e028 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -50,6 +50,7 @@ #include #include #include +#include //#include // A debugging utility. @@ -1373,11 +1374,26 @@ namespace detail { const ADB& liquid = subset(state.qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); const ADB& vapour = subset(state.qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); - V bhp_targets = V::Zero(nw); - V rate_targets = V::Zero(nw); - M rate_distr(nw, np*nw); + //1. Calculate THP targets + std::vector table_id(nw, -1); + ADB::V thp_v = ADB::V::Zero(nw); + ADB::V alq_v = ADB::V::Zero(nw); + + //Target vars + ADB::V bhp_targets = ADB::V::Zero(nw); + ADB::V rate_targets = ADB::V::Zero(nw); + ADB::M rate_distr(nw, np*nw); + + //Selection variables + std::vector bhp_elems; + std::vector thp_elems; + std::vector rate_elems; + + //Run through all wells to calculate BHP/RATE targets + //and gather info about current control for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; + auto wc = wells().ctrls[w]; + // The current control in the well state overrides // the current control set in the Wells struct, which // is instead treated as a default. @@ -1386,18 +1402,20 @@ namespace detail { switch (well_controls_iget_type(wc, current)) { case BHP: { - bhp_targets (w) = well_controls_iget_target(wc, current); + bhp_elems.push_back(w); + bhp_targets(w) = well_controls_iget_target(wc, current); rate_targets(w) = -1e100; } break; case THP: { - const int vfp = well_controls_iget_vfp(wc, current); - const double& thp = well_controls_iget_target(wc, current); - const double& alq = well_controls_iget_alq(wc, current); + table_id[w] = well_controls_iget_vfp(wc, current); + thp_v[w] = well_controls_iget_target(wc, current); + alq_v[w] = well_controls_iget_alq(wc, current); - bhp_targets (w) = vfp_properties_->getProd()->bhp(vfp, aqua.value()[w], liquid.value()[w], vapour.value()[w], thp, alq).value; + thp_elems.push_back(w); + bhp_targets(w) = -1e100; rate_targets(w) = -1e100; } break; @@ -1405,6 +1423,7 @@ namespace detail { case RESERVOIR_RATE: // Intentional fall-through case SURFACE_RATE: { + rate_elems.push_back(w); // RESERVOIR and SURFACE rates look the same, from a // high-level point of view, in the system of // simultaneous linear equations. @@ -1416,19 +1435,28 @@ namespace detail { rate_distr.insert(w, p*nw + w) = distr[p]; } - bhp_targets (w) = -1.0e100; + bhp_targets(w) = -1.0e100; rate_targets(w) = well_controls_iget_target(wc, current); } break; } } + + //Calculate BHP target from THP + const ADB thp = ADB::constant(thp_v); + const ADB alq = ADB::constant(alq_v); + const ADB thp_targets = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq); + + //Calculate residuals + const ADB thp_residual = state.bhp - thp_targets; const ADB bhp_residual = state.bhp - bhp_targets; const ADB rate_residual = rate_distr * state.qs - rate_targets; - //wells - // Choose bhp residual for positive bhp targets. - Selector bhp_selector(bhp_targets); - residual_.well_eq = bhp_selector.select(bhp_residual, rate_residual); + //Select the right residual for each well + residual_.well_eq = superset(subset(bhp_residual, bhp_elems), bhp_elems, bhp_residual.size()) + + superset(subset(thp_residual, thp_elems), thp_elems, thp_residual.size()) + + superset(subset(rate_residual, rate_elems), rate_elems, rate_residual.size()); + // For wells that are dead (not flowing), and therefore not communicating // with the reservoir, we set the equation to be equal to the well's total // flow. This will be a solution only if the target rate is also zero. diff --git a/opm/autodiff/VFPProperties.cpp b/opm/autodiff/VFPProperties.cpp index 66a65ad39..9f336467d 100644 --- a/opm/autodiff/VFPProperties.cpp +++ b/opm/autodiff/VFPProperties.cpp @@ -25,6 +25,8 @@ #include #include +#include +#include namespace Opm { @@ -99,7 +101,7 @@ void VFPProdProperties::init(const std::map& prod_tables) { } } -VFPProdProperties::ADB VFPProdProperties::bhp(int table_id, +VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, const Wells& wells, const ADB& qs, const ADB& thp, @@ -117,52 +119,231 @@ VFPProdProperties::ADB VFPProdProperties::bhp(int table_id, return bhp(table_id, w, o, g, thp, alq); } -VFPProdProperties::ADB VFPProdProperties::bhp(int table_id, +namespace detail { + /** + * Returns the type variable for FLO/GFR/WFR + */ + template + TYPE getType(const VFPProdTable* table); + + template <> + VFPProdTable::FLO_TYPE getType(const VFPProdTable* table) { + return table->getFloType(); + } + + template <> + VFPProdTable::WFR_TYPE getType(const VFPProdTable* table) { + return table->getWFRType(); + } + + template <> + VFPProdTable::GFR_TYPE getType(const VFPProdTable* table) { + return table->getGFRType(); + } + + /** + * Returns the actual ADB for the type of FLO/GFR/WFR type + */ + template + VFPProdProperties::ADB getValue( + const VFPProdProperties::ADB& aqua, + const VFPProdProperties::ADB& liquid, + const VFPProdProperties::ADB& vapour, TYPE type); + + template <> + VFPProdProperties::ADB getValue( + const VFPProdProperties::ADB& aqua, + const VFPProdProperties::ADB& liquid, + const VFPProdProperties::ADB& vapour, + VFPProdTable::FLO_TYPE type) { + return VFPProdProperties::getFlo(aqua, liquid, vapour, type); + } + + template <> + VFPProdProperties::ADB getValue( + const VFPProdProperties::ADB& aqua, + const VFPProdProperties::ADB& liquid, + const VFPProdProperties::ADB& vapour, + VFPProdTable::WFR_TYPE type) { + return VFPProdProperties::getWFR(aqua, liquid, vapour, type); + } + + template <> + VFPProdProperties::ADB getValue( + const VFPProdProperties::ADB& aqua, + const VFPProdProperties::ADB& liquid, + const VFPProdProperties::ADB& vapour, + VFPProdTable::GFR_TYPE type) { + return VFPProdProperties::getGFR(aqua, liquid, vapour, type); + } + + /** + * Given m wells and n types of VFP variables (e.g., FLO = {FLO_OIL, FLO_LIQ} + * this function combines the n types of ADB objects, so that each of the + * m wells gets the right ADB. + */ + template + VFPProdProperties::ADB gather_vars(const std::vector& well_tables, + const VFPProdProperties::ADB& aqua, + const VFPProdProperties::ADB& liquid, + const VFPProdProperties::ADB& vapour) { + + typedef VFPProdProperties::ADB ADB; + + const int num_wells = static_cast(well_tables.size()); + assert(aqua.size() == num_wells); + assert(liquid.size() == num_wells); + assert(vapour.size() == num_wells); + + //Caching variable for flo/wfr/gfr + std::map map; + + //Indexing variable used when combining the different ADB types + std::map > elems; + + //Compute all of the different ADB types, + //and record which wells use which types + for (int i=0; i(table); + + //"Caching" of flo_type etc: Only calculate used types + //Create type if it does not exist + if (map.find(type) == map.end()) { + map.insert(std::pair( + type, + detail::getValue(aqua, liquid, vapour, type) + )); + } + + //Add the index for assembly later in gather_vars + elems[type].push_back(i); + } + } + + //Loop over all types of ADB variables, and combine them + //so that each well gets the proper variable + ADB retval = ADB::constant(ADB::V::Zero(num_wells)); + for (const auto& entry : elems) { + const auto& key = entry.first; + const auto& value = entry.second; + + //Get the ADB for this type of variable + assert(map.find(key) != map.end()); + const ADB& values = map.find(key)->second; + + //Get indices to all elements that should use this ADB + const std::vector& elems = value; + + //Add these elements to retval + retval = retval + superset(subset(values, elems), elems, values.size()); + } + + return retval; + } + + void extendBlockPattern(const VFPProdProperties::ADB& x, std::vector& block_pattern) { + std::vector x_block_pattern = x.blockPattern(); + + if (x_block_pattern.empty()) { + return; + } + else { + if (block_pattern.empty()) { + block_pattern = x_block_pattern; + return; + } + else { + if (x_block_pattern != block_pattern) { + OPM_THROW(std::logic_error, "Block patterns do not match"); + } + } + } + } + + std::vector commonBlockPattern( + const VFPProdProperties::ADB& x1, + const VFPProdProperties::ADB& x2, + const VFPProdProperties::ADB& x3, + const VFPProdProperties::ADB& x4, + const VFPProdProperties::ADB& x5) { + std::vector block_pattern; + + extendBlockPattern(x1, block_pattern); + extendBlockPattern(x2, block_pattern); + extendBlockPattern(x3, block_pattern); + extendBlockPattern(x4, block_pattern); + extendBlockPattern(x5, block_pattern); + + return block_pattern; + } + +} //Namespace + +VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, const ADB& aqua, const ADB& liquid, const ADB& vapour, const ADB& thp, const ADB& alq) const { - const VFPProdTable* table = getProdTable(table_id); const int nw = thp.size(); - assert(aqua.size() == nw); - assert(liquid.size() == nw); - assert(vapour.size() == nw); - assert(thp.size() == nw); - assert(alq.size() == nw); + std::vector block_pattern = detail::commonBlockPattern(aqua, liquid, vapour, thp, alq); + + assert(static_cast(table_id.size()) == nw); + assert(aqua.size() == nw); + assert(liquid.size() == nw); + assert(vapour.size() == nw); + assert(thp.size() == nw); + assert(alq.size() == nw); //Allocate data for bhp's and partial derivatives - ADB::V value, dthp, dwfr, dgfr, dalq, dflo; - value.resize(nw); - dthp.resize(nw); - dwfr.resize(nw); - dgfr.resize(nw); - dalq.resize(nw); - dflo.resize(nw); + ADB::V value = ADB::V::Zero(nw); + ADB::V dthp = ADB::V::Zero(nw); + ADB::V dwfr = ADB::V::Zero(nw); + ADB::V dgfr = ADB::V::Zero(nw); + ADB::V dalq = ADB::V::Zero(nw); + ADB::V dflo = ADB::V::Zero(nw); - //Find interpolation variables - ADB flo = getFlo(aqua, liquid, vapour, table->getFloType()); - ADB wfr = getWFR(aqua, liquid, vapour, table->getWFRType()); - ADB gfr = getGFR(aqua, liquid, vapour, table->getGFRType()); + //Get the table for each well + std::vector well_tables(nw, NULL); + for (int i=0; i= 0) { + well_tables[i] = getProdTable(table_id[i]); + } + } + + //Get the right FLO/GFR/WFR variable for each well as a single ADB + const ADB flo = detail::gather_vars(well_tables, aqua, liquid, vapour); + const ADB wfr = detail::gather_vars(well_tables, aqua, liquid, vapour); + const ADB gfr = detail::gather_vars(well_tables, aqua, liquid, vapour); //Compute the BHP for each well independently for (int i=0; igetFloAxis()); - auto thp_i = find_interp_data(thp.value()[i], table->getTHPAxis()); - auto wfr_i = find_interp_data(wfr.value()[i], table->getWFRAxis()); - auto gfr_i = find_interp_data(gfr.value()[i], table->getGFRAxis()); - auto alq_i = find_interp_data(alq.value()[i], table->getALQAxis()); + const VFPProdTable* table = well_tables[i]; + if (table != NULL) { + //First, find the values to interpolate between + auto flo_i = find_interp_data(flo.value()[i], table->getFloAxis()); + auto thp_i = find_interp_data(thp.value()[i], table->getTHPAxis()); + auto wfr_i = find_interp_data(wfr.value()[i], table->getWFRAxis()); + auto gfr_i = find_interp_data(gfr.value()[i], table->getGFRAxis()); + auto alq_i = find_interp_data(alq.value()[i], table->getALQAxis()); - adb_like bhp_val = interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); + adb_like bhp_val = interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); - value[i] = bhp_val.value; - dthp[i] = bhp_val.dthp; - dwfr[i] = bhp_val.dwfr; - dgfr[i] = bhp_val.dgfr; - dalq[i] = bhp_val.dalq; - dflo[i] = bhp_val.dflo; + value[i] = bhp_val.value; + dthp[i] = bhp_val.dthp; + dwfr[i] = bhp_val.dwfr; + dgfr[i] = bhp_val.dgfr; + dalq[i] = bhp_val.dalq; + dflo[i] = bhp_val.dflo; + } + else { + value[i] = -1e100; //Signal that this value has not been calculated properly, due to "missing" table + } } //Create diagonal matrices from ADB::Vs @@ -172,17 +353,29 @@ VFPProdProperties::ADB VFPProdProperties::bhp(int table_id, ADB::M dalq_diag = spdiag(dalq); ADB::M dflo_diag = spdiag(dflo); - //Calculate the jacobians - const int num_blocks = aqua.numBlocks(); + //Calculate the Jacobians + const int num_blocks = block_pattern.size(); std::vector jacs(num_blocks); for (int block = 0; block < num_blocks; ++block) { //Could have used fastSparseProduct and temporary variables //but may not save too much on that. - jacs[block] = dthp_diag * thp.derivative()[block] + - dwfr_diag * wfr.derivative()[block] + - dgfr_diag * gfr.derivative()[block] + - dalq_diag * alq.derivative()[block] + - dflo_diag * flo.derivative()[block]; + jacs[block] = ADB::M(nw, block_pattern[block]); + + if (!thp.derivative().empty()) { + jacs[block] += dthp_diag * thp.derivative()[block]; + } + if (!wfr.derivative().empty()) { + jacs[block] += dwfr_diag * wfr.derivative()[block]; + } + if (!gfr.derivative().empty()) { + jacs[block] += dgfr_diag * gfr.derivative()[block]; + } + if (!alq.derivative().empty()) { + jacs[block] += dalq_diag * alq.derivative()[block]; + } + if (!flo.derivative().empty()) { + jacs[block] += dflo_diag * flo.derivative()[block]; + } } ADB retval = ADB::function(std::move(value), std::move(jacs)); diff --git a/opm/autodiff/VFPProperties.hpp b/opm/autodiff/VFPProperties.hpp index b416f0648..986546b5c 100644 --- a/opm/autodiff/VFPProperties.hpp +++ b/opm/autodiff/VFPProperties.hpp @@ -139,7 +139,7 @@ public: * @return The bottom hole pressure, interpolated/extrapolated linearly using * the above parameters from the values in the input table. */ - ADB bhp(int table_id, + ADB bhp(const std::vector& table_id, const Wells& wells, const ADB& qs, const ADB& thp, @@ -147,7 +147,10 @@ public: /** * Linear interpolation of bhp as a function of the input parameters given as ADBs - * @param table Table number to use + * Each entry corresponds typically to one well. + * @param table_id Table number to use. A negative entry (e.g., -1) + * will indicate that no table is used, and the corresponding + * BHP will be calculated as a constant -1e100. * @param aqua Water phase * @param liquid Oil phase * @param vapour Gas phase @@ -158,7 +161,7 @@ public: * the above parameters from the values in the input table, for each entry in the * input ADB objects. */ - ADB bhp(int table_id, + ADB bhp(const std::vector& table_id, const ADB& aqua, const ADB& liquid, const ADB& vapour, @@ -203,8 +206,6 @@ public: const double& bhp, const double& alq) const; - //FIXME: ARB: Implement inj_bhp to match the prod_bhp's, but for injection wells. - /** * Computes the flo parameter according to the flo_type_ * @return Production rate of oil, gas or liquid. diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 6847e5369..9aa61bdcd 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -57,7 +57,8 @@ struct TrivialFixture { typedef Opm::VFPProdProperties::ADB ADB; typedef Opm::VFPProdProperties::adb_like adb_like; - TrivialFixture() : thp_axis{0.0, 1.0}, + TrivialFixture() : table_ids(1, 1), + thp_axis{0.0, 1.0}, wfr_axis{0.0, 0.5, 1.0}, gfr_axis{0.0, 0.25, 0.5, 0.75, 1}, alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}, @@ -145,7 +146,7 @@ struct TrivialFixture { inline void initProperties() { //Initialize table - table.init(1, + table.init(table_ids[0], 1000.0, Opm::VFPProdTable::FLO_OIL, Opm::VFPProdTable::WFR_WOR, @@ -175,6 +176,7 @@ struct TrivialFixture { std::shared_ptr properties; Opm::VFPProdTable table; + std::vector table_ids; private: const std::vector thp_axis; @@ -232,9 +234,9 @@ BOOST_AUTO_TEST_CASE(GetTable) //Check that different versions of the prod_bph function work - ADB a = properties->bhp(1, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb); - adb_like b = properties->bhp(1, aqua_d, liquid_d, vapour_d, thp_d, alq_d); - ADB c = properties->bhp(1, *wells, qs_adb, thp_adb, alq_adb); + ADB a = properties->bhp(table_ids, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb); + adb_like b = properties->bhp(table_ids[0], aqua_d, liquid_d, vapour_d, thp_d, alq_d); + ADB c = properties->bhp(table_ids, *wells, qs_adb, thp_adb, alq_adb); //Check that results are actually equal double d = a.value()[0]; @@ -244,7 +246,8 @@ BOOST_AUTO_TEST_CASE(GetTable) BOOST_CHECK_EQUAL(d, f); //Table 2 does not exist. - BOOST_CHECK_THROW(properties->bhp(2, *wells, qs_adb, thp_adb, alq_adb), std::invalid_argument); + std::vector table_ids_wrong(1, 2); + BOOST_CHECK_THROW(properties->bhp(table_ids_wrong, *wells, qs_adb, thp_adb, alq_adb), std::invalid_argument); } /** @@ -510,7 +513,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB) ADB adb_u = ADB::constant(adb_v_u); ADB adb_liquid = ADB::constant(adb_v_liquid); - ADB bhp = properties->bhp(1, adb_aqua, adb_liquid, adb_vapour, adb_x, adb_u); + ADB bhp = properties->bhp(table_ids, adb_aqua, adb_liquid, adb_vapour, adb_x, adb_u); ADB::V bhp_val = bhp.value(); double value = 0.0; @@ -596,7 +599,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs) ADB alq = ADB::constant(alq_v); //Call the bhp function - ADB::V bhp = properties->bhp(1, *wells, qs, thp, alq).value(); + ADB::V bhp = properties->bhp(table_ids, *wells, qs, thp, alq).value(); //Calculate reference //First, find the three phases