From 34edf3a5b8819c573f626bf6403350aeb72aee57 Mon Sep 17 00:00:00 2001 From: babrodtk Date: Thu, 6 Aug 2015 09:25:41 +0200 Subject: [PATCH] Changed API of VFPProperties to take ADBs --- opm/autodiff/BlackoilModelBase_impl.hpp | 2 +- opm/autodiff/VFPProperties.cpp | 154 ++++++++++++++---------- opm/autodiff/VFPProperties.hpp | 37 ++++-- tests/test_vfpproperties.cpp | 123 +++++++++++-------- 4 files changed, 191 insertions(+), 125 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index bbe4090fa..213f13e9c 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1379,7 +1379,7 @@ namespace detail { const double& thp = well_controls_iget_target(wc, current); const double& alq = well_controls_iget_alq(wc, current); - bhp_targets (w) = vfp_properties_->getProd()->bhp(vfp, aqua[w], liquid[w], vapour[w], thp, alq); + bhp_targets (w) = vfp_properties_->getProd()->bhp(vfp, aqua[w], liquid[w], vapour[w], thp, alq).value; rate_targets(w) = -1e100; } break; diff --git a/opm/autodiff/VFPProperties.cpp b/opm/autodiff/VFPProperties.cpp index 3ed814bb0..cf1e056d1 100644 --- a/opm/autodiff/VFPProperties.cpp +++ b/opm/autodiff/VFPProperties.cpp @@ -99,30 +99,30 @@ void VFPProdProperties::init(const std::map& prod_tables) { } } -VFPProdProperties::ADB::V VFPProdProperties::bhp(int table_id, +VFPProdProperties::ADB VFPProdProperties::bhp(int table_id, const Wells& wells, - const ADB::V& qs, - const ADB::V& thp, - const ADB::V& alq) const { + const ADB& qs, + const ADB& thp, + const ADB& alq) const { const int np = wells.number_of_phases; const int nw = wells.number_of_wells; //Short-hands for water / oil / gas phases //TODO enable support for two-phase. assert(np == 3); - const ADB::V& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); - const ADB::V& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); - const ADB::V& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); + const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); + const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); + const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); return bhp(table_id, w, o, g, thp, alq); } -VFPProdProperties::ADB::V VFPProdProperties::bhp(int table_id, - const ADB::V& aqua, - const ADB::V& liquid, - const ADB::V& vapour, - const ADB::V& thp, - const ADB::V& alq) const { +VFPProdProperties::ADB VFPProdProperties::bhp(int table_id, + const ADB& aqua, + const ADB& liquid, + const ADB& vapour, + const ADB& thp, + const ADB& alq) const { const int nw = thp.size(); assert(aqua.size() == nw); @@ -131,16 +131,53 @@ VFPProdProperties::ADB::V VFPProdProperties::bhp(int table_id, 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); + //Compute the BHP for each well independently - ADB::V bhp_vals; - bhp_vals.resize(nw); for (int i=0; i jacs(num_blocks); + for (int block = 0; block < num_blocks; ++block) { + fastSparseProduct(dthp_diag, aqua.derivative()[block], jacs[block]); + ADB::M temp; + fastSparseProduct(dwfr_diag, liquid.derivative()[block], temp); + jacs[block] += temp; + } + + ADB retval = ADB::function(std::move(value), std::move(jacs)); + return retval; } -double VFPProdProperties::bhp(int table_id, +VFPProdProperties::adb_like VFPProdProperties::bhp(int table_id, const double& aqua, const double& liquid, const double& vapour, @@ -161,7 +198,8 @@ double VFPProdProperties::bhp(int table_id, auto alq_i = find_interp_data(alq, table->getALQAxis()); //Then perform the interpolation itself - return interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); + adb_like retval = interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); + return retval; } @@ -200,7 +238,7 @@ double VFPProdProperties::thp(int table_id, std::vector bhp_array(nthp); for (int i=0; i ADB; + /** + * An "ADB-like" structure with a single value and a set of derivatives + */ + struct adb_like { + adb_like() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {}; + double value; + double dthp; + double dwfr; + double dgfr; + double dalq; + double dflo; + }; + /** * Empty constructor */ @@ -125,11 +138,11 @@ public: * @return The bottom hole pressure, interpolated/extrapolated linearly using * the above parameters from the values in the input table. */ - ADB::V bhp(int table_id, + ADB bhp(int table_id, const Wells& wells, - const ADB::V& qs, - const ADB::V& thp, - const ADB::V& alq) const; + const ADB& qs, + const ADB& thp, + const ADB& alq) const; /** * Linear interpolation of bhp as a function of the input parameters given as ADBs @@ -144,12 +157,12 @@ public: * the above parameters from the values in the input table, for each entry in the * input ADB objects. */ - ADB::V bhp(int table_id, - const ADB::V& aqua, - const ADB::V& liquid, - const ADB::V& vapour, - const ADB::V& thp, - const ADB::V& alq) const; + ADB bhp(int table_id, + const ADB& aqua, + const ADB& liquid, + const ADB& vapour, + const ADB& thp, + const ADB& alq) const; /** * Linear interpolation of bhp as a function of the input parameters @@ -163,7 +176,7 @@ public: * @return The bottom hole pressure, interpolated/extrapolated linearly using * the above parameters from the values in the input table. */ - double bhp(int table_id, + adb_like bhp(int table_id, const double& aqua, const double& liquid, const double& vapour, @@ -294,7 +307,7 @@ private: /** * Helper function which interpolates data using the indices etc. given in the inputs. */ - static double interpolate(const VFPProdTable::array_type& array, + static adb_like interpolate(const VFPProdTable::array_type& array, const InterpData& flo_i, const InterpData& thp_i, const InterpData& wfr_i, diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index b8b40c7eb..70de3f560 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -55,6 +55,7 @@ const double sad_tol = 1.0e-8; */ struct TrivialFixture { typedef Opm::VFPProdProperties::ADB ADB; + typedef Opm::VFPProdProperties::adb_like adb_like; TrivialFixture() : thp_axis{0.0, 1.0}, wfr_axis{0.0, 0.5, 1.0}, @@ -80,7 +81,7 @@ struct TrivialFixture { /** * Fills our interpolation data with zeros */ - void fillData(double value) { + inline void fillData(double value) { for (int i=0; i(nx-1); for (int j=0; j properties; Opm::VFPProdTable table; @@ -208,25 +220,28 @@ BOOST_AUTO_TEST_CASE(GetTable) double thp_d = 4.5; double alq_d = 5.5; - ADB::V aqua_adb = ADB::V::Constant(1, aqua_d); - ADB::V liquid_adb = ADB::V::Constant(1, liquid_d); - ADB::V vapour_adb = ADB::V::Constant(1, vapour_d); - ADB::V thp_adb = ADB::V::Constant(1, thp_d); - ADB::V alq_adb = ADB::V::Constant(1, alq_d); - - ADB::V qs_adb(3); - qs_adb << aqua_adb, liquid_adb, vapour_adb; + ADB aqua_adb = createConstantScalarADB(aqua_d); + ADB liquid_adb = createConstantScalarADB(liquid_d); + ADB vapour_adb = createConstantScalarADB(vapour_d); + ADB thp_adb = createConstantScalarADB(thp_d); + ADB alq_adb = createConstantScalarADB(alq_d); + ADB::V qs_adb_v(3); + qs_adb_v << aqua_adb.value(), liquid_adb.value(), vapour_adb.value(); + ADB qs_adb = ADB::constant(qs_adb_v); //Check that different versions of the prod_bph function work - ADB::V a = properties->bhp(1, aqua_adb, liquid_adb, vapour_adb, thp_adb, alq_adb); - double b = properties->bhp(1, aqua_d, liquid_d, vapour_d, thp_d, alq_d); - ADB::V c = properties->bhp(1, *wells, qs_adb, thp_adb, alq_adb); + 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); //Check that results are actually equal - BOOST_CHECK_EQUAL(a[0], b); - BOOST_CHECK_EQUAL(a[0], c[0]); + double d = a.value()[0]; + double e = b.value; + double f = c.value()[0]; + BOOST_CHECK_EQUAL(d, e); + 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); @@ -256,7 +271,7 @@ BOOST_AUTO_TEST_CASE(InterpolateZero) const double v = m / static_cast(n-1); //Note order of arguments! - sum += properties->bhp(1, v, x, y, z, u); + sum += properties->bhp(1, v, x, y, z, u).value; } } } @@ -291,7 +306,7 @@ BOOST_AUTO_TEST_CASE(InterpolateOne) const double v = m / static_cast(n-1); //Note order of arguments! - const double value = properties->bhp(1, v, x, y, z, u); + const double value = properties->bhp(1, v, x, y, z, u).value; sum += value; } @@ -340,7 +355,7 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane) reference_sum += reference; //Note order of arguments! id, aqua, liquid, vapour, thp , alq - double value = properties->bhp(1, aqua, liquid, vapour, x, u); + double value = properties->bhp(1, aqua, liquid, vapour, x, u).value; sum += value; double abs_diff = std::abs(value - reference); @@ -396,7 +411,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane) reference_sum += reference; //Note order of arguments! id, aqua, liquid, vapour, thp , alq - double value = properties->bhp(1, aqua, liquid, vapour, x, u); + double value = properties->bhp(1, aqua, liquid, vapour, x, u).value; sum += value; double abs_diff = std::abs(value - reference); @@ -426,14 +441,6 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB) fillDataPlane(); initProperties(); - //Temporary variables used to represent independent wells - const int num_wells = 5; - ADB::V x_v(num_wells); - ADB::V aqua_v(num_wells); - ADB::V vapour_v(num_wells); - ADB::V u_v(num_wells); - ADB::V liquid_v(num_wells); - //Check linear extrapolation (i.e., using values of x, y, etc. outside our interpolant domain) double sum = 0.0; double reference_sum = 0.0; @@ -452,15 +459,30 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB) for (int m=1; m<=n+o; ++m) { const double liquid = m / static_cast(n); + //Temporary variables used to represent independent wells + const int num_wells = 5; + ADB::V adb_v_x(num_wells); + ADB::V adb_v_aqua(num_wells); + ADB::V adb_v_vapour(num_wells); + ADB::V adb_v_u(num_wells); + ADB::V adb_v_liquid(num_wells); + for (unsigned int w=0; wbhp(1, aqua_v, liquid_v, vapour_v, x_v, u_v); + ADB adb_x = ADB::constant(adb_v_x); + ADB adb_aqua = ADB::constant(adb_v_aqua); + ADB adb_vapour = ADB::constant(adb_v_vapour); + 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::V bhp_val = bhp.value(); double value = 0.0; double reference = 0.0; @@ -522,27 +544,30 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs) } //Create some artificial flow values for our wells between 0 and 1 - ADB::V qs(nphases*nwells); + ADB::V qs_v(nphases*nwells); for (int j=0; j(nwells*nphases-1.0); + qs_v[j*nwells+i] = (j*nwells+i) / static_cast(nwells*nphases-1.0); } } + ADB qs = ADB::constant(qs_v); //Create the THP for each well - ADB::V thp(nwells); + ADB::V thp_v(nwells); for (int i=0; i(nwells-1.0); + thp_v[i] = (i) / static_cast(nwells-1.0); } + ADB thp = ADB::constant(thp_v); //Create the ALQ for each well - ADB::V alq(nwells); + ADB::V alq_v(nwells); for (int i=0; ibhp(1, *wells, qs, thp, alq); + ADB::V bhp = properties->bhp(1, *wells, qs, thp, alq).value(); //Calculate reference //First, find the three phases @@ -550,9 +575,9 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs) std::vector oil(nwells); std::vector gas(nwells); for (int i=0; i