From 3bcfc905bdefebfadaef8a19209dc0cff8861bed Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 7 Apr 2016 15:41:08 +0200 Subject: [PATCH 01/13] makding StandardWells a template based on SolutionState and WellState may need to be adjusted later. --- CMakeLists_files.cmake | 5 +- opm/autodiff/BlackoilModelBase.hpp | 6 +-- opm/autodiff/StandardWells.hpp | 4 ++ ...andardWells.cpp => StandardWells_impl.hpp} | 53 ++++++++++++++----- 4 files changed, 52 insertions(+), 16 deletions(-) rename opm/autodiff/{StandardWells.cpp => StandardWells_impl.hpp} (63%) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 893bffa12..255fae7f0 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -26,6 +26,7 @@ # originally generated with the command: # find opm -name '*.c*' -printf '\t%p\n' | sort list (APPEND MAIN_SOURCE_FILES + opm/autodiff/BlackoilPropsAdInterface.cpp opm/autodiff/ExtractParallelGridInformationToISTL.cpp opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -47,7 +48,6 @@ list (APPEND MAIN_SOURCE_FILES opm/autodiff/VFPProdProperties.cpp opm/autodiff/VFPInjProperties.cpp opm/autodiff/WellMultiSegment.cpp - opm/autodiff/StandardWells.cpp opm/autodiff/BlackoilSolventState.cpp opm/autodiff/ThreadHandle.hpp opm/polymer/PolymerState.cpp @@ -68,6 +68,7 @@ list (APPEND MAIN_SOURCE_FILES ) + # originally generated with the command: # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort list (APPEND TEST_SOURCE_FILES @@ -130,6 +131,7 @@ list (APPEND PROGRAM_SOURCE_FILES # originally generated with the command: # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort list (APPEND PUBLIC_HEADER_FILES +<<<<<<< HEAD opm/autodiff/AdditionalObjectDeleter.hpp opm/autodiff/AutoDiffBlock.hpp opm/autodiff/AutoDiffHelpers.hpp @@ -195,6 +197,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/WellStateMultiSegment.hpp opm/autodiff/WellMultiSegment.hpp opm/autodiff/StandardWells.hpp + opm/autodiff/StandardWells_impl.hpp opm/polymer/CompressibleTpfaPolymer.hpp opm/polymer/GravityColumnSolverPolymer.hpp opm/polymer/GravityColumnSolverPolymer_impl.hpp diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index b8ceae236..1a95fdb31 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -276,7 +276,7 @@ namespace Opm { const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_; const RockCompressibility* rock_comp_props_; - StandardWells std_wells_; + StandardWells std_wells_; VFPProperties vfp_properties_; const NewtonIterationBlackoilInterface& linsolver_; // For each canonical phase -> true if active @@ -329,8 +329,8 @@ namespace Opm { } /// return the StandardWells object - StandardWells& stdWells() { return std_wells_; } - const StandardWells& stdWells() const { return std_wells_; } + StandardWells& stdWells() { return std_wells_; } + const StandardWells& stdWells() const { return std_wells_; } /// return the Well struct in the StandardWells const Wells& wells() const { return std_wells_.wells(); } diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 1c3d556b7..eda81114b 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -40,6 +40,7 @@ namespace Opm { typedef ADB::V Vector; /// Class for handling the standard well model. + template class StandardWells { protected: struct WellOps { @@ -81,4 +82,7 @@ namespace Opm { } // namespace Opm + +#include "StandardWells_impl.hpp" + #endif diff --git a/opm/autodiff/StandardWells.cpp b/opm/autodiff/StandardWells_impl.hpp similarity index 63% rename from opm/autodiff/StandardWells.cpp rename to opm/autodiff/StandardWells_impl.hpp index 6b49be6ab..b7b418c51 100644 --- a/opm/autodiff/StandardWells.cpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -27,7 +27,8 @@ namespace Opm { - StandardWells:: + template + StandardWells:: WellOps::WellOps(const Wells* wells) : w2p(), p2w(), @@ -65,7 +66,9 @@ namespace Opm - StandardWells::StandardWells(const Wells* wells_arg) + template + StandardWells:: + StandardWells(const Wells* wells_arg) : wells_(wells_arg) , wops_(wells_arg) , well_perforation_densities_(Vector()) @@ -77,7 +80,10 @@ namespace Opm - const Wells& StandardWells::wells() const + template + const Wells& + StandardWells:: + wells() const { assert(wells_ != 0); return *(wells_); @@ -87,7 +93,10 @@ namespace Opm - bool StandardWells::wellsActive() const + template + bool + StandardWells:: + wellsActive() const { return wells_active_; } @@ -96,7 +105,10 @@ namespace Opm - void StandardWells::setWellsActive(const bool wells_active) + template + void + StandardWells:: + setWellsActive(const bool wells_active) { wells_active_ = wells_active; } @@ -105,7 +117,10 @@ namespace Opm - bool StandardWells::localWellsActive() const + template + bool + StandardWells:: + localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; } @@ -114,8 +129,10 @@ namespace Opm - const StandardWells::WellOps& - StandardWells::wellOps() const + template + const typename StandardWells::WellOps& + StandardWells:: + wellOps() const { return wops_; } @@ -124,7 +141,10 @@ namespace Opm - Vector& StandardWells::wellPerforationDensities() + template + Vector& + StandardWells:: + wellPerforationDensities() { return well_perforation_densities_; } @@ -133,7 +153,10 @@ namespace Opm - const Vector& StandardWells::wellPerforationDensities() const + template + const Vector& + StandardWells:: + wellPerforationDensities() const { return well_perforation_densities_; } @@ -142,7 +165,10 @@ namespace Opm - Vector& StandardWells::wellPerforationPressureDiffs() + template + Vector& + StandardWells:: + wellPerforationPressureDiffs() { return well_perforation_pressure_diffs_; } @@ -151,7 +177,10 @@ namespace Opm - const Vector& StandardWells::wellPerforationPressureDiffs() const + template + const Vector& + StandardWells:: + wellPerforationPressureDiffs() const { return well_perforation_pressure_diffs_; } From 266f6f2df5f3cb232ea7e935a42ac4c6fb5a8d78 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Thu, 7 Apr 2016 16:17:25 +0200 Subject: [PATCH 02/13] adding computePropertiesForWellConnectionPressures to StandardWells tested with flow. Not sure how other flow siblings work. --- opm/autodiff/BlackoilModelBase_impl.hpp | 2 +- opm/autodiff/StandardWells.hpp | 20 ++++++ opm/autodiff/StandardWells_impl.hpp | 88 +++++++++++++++++++++++++ 3 files changed, 109 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 8711544b1..363b5a3d7 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -849,7 +849,7 @@ namespace detail { std::vector rsmax_perf; std::vector rvmax_perf; std::vector surf_dens_perf; - asImpl().computePropertiesForWellConnectionPressures(state, xw, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // 2. Compute densities std::vector cd = diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index eda81114b..e57f6dee0 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -31,6 +31,8 @@ #include #include +#include +#include namespace Opm { @@ -39,6 +41,13 @@ namespace Opm { typedef AutoDiffBlock ADB; typedef ADB::V Vector; + // copied from BlackoilModelBase + // should put to somewhere better + typedef Eigen::Array DataBlock; + /// Class for handling the standard well model. template class StandardWells { @@ -72,6 +81,17 @@ namespace Opm { Vector& wellPerforationPressureDiffs(); const Vector& wellPerforationPressureDiffs() const; + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); + + protected: bool wells_active_; const Wells* wells_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index b7b418c51..16478b5f2 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -23,6 +23,7 @@ + namespace Opm { @@ -185,4 +186,91 @@ namespace Opm return well_perforation_pressure_diffs_; } + + + + template + void + StandardWells:: + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + + // Compute the average pressure in each well block + const Vector perf_press = Eigen::Map(xw.perfPress().data(), nperf); + Vector avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + } + + const std::vector& well_cells = wellOps().well_cells; + + // Use cell values for the temperature as the wells don't knows its temperature yet. + const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); + std::vector perf_cond(nperf); + // const std::vector& pc = phaseCondition(); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = pc[well_cells[perf]]; + } + const PhaseUsage& pu = fluid.phaseUsage(); + DataBlock b(nperf, pu.num_phases); + if (pu.phase_used[BlackoilPhases::Aqua]) { + const Vector bw = fluid.bWat(avg_press_ad, perf_temp, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; + } + assert(active[Oil]); + const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); + if (pu.phase_used[BlackoilPhases::Liquid]) { + const ADB perf_rs = subset(state.rs, well_cells); + const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + // const V rssat = fluidRsSat(avg_press, perf_so, well_cells); + const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + } else { + rsmax_perf.assign(nperf, 0.0); + } + if (pu.phase_used[BlackoilPhases::Vapour]) { + const ADB perf_rv = subset(state.rv, well_cells); + const Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + // const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); + const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); + } else { + rvmax_perf.assign(nperf, 0.0); + } + // b is row major, so can just copy data. + b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + + // Surface density. + // The compute density segment wants the surface densities as + // an np * number of wells cells array + Vector rho = superset(fluid.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + rho += superset(fluid.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); + + } + + } From bca23d34c1afc10e62b44eb782fd778dc09f7c5c Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Apr 2016 10:20:11 +0200 Subject: [PATCH 03/13] making StandardWells non-templated we will focus on template functions instead. --- opm/autodiff/BlackoilModelBase.hpp | 6 ++-- opm/autodiff/StandardWells.hpp | 2 +- opm/autodiff/StandardWells_impl.hpp | 52 ++++++++--------------------- 3 files changed, 17 insertions(+), 43 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 1a95fdb31..b8ceae236 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -276,7 +276,7 @@ namespace Opm { const BlackoilPropsAdInterface& fluid_; const DerivedGeology& geo_; const RockCompressibility* rock_comp_props_; - StandardWells std_wells_; + StandardWells std_wells_; VFPProperties vfp_properties_; const NewtonIterationBlackoilInterface& linsolver_; // For each canonical phase -> true if active @@ -329,8 +329,8 @@ namespace Opm { } /// return the StandardWells object - StandardWells& stdWells() { return std_wells_; } - const StandardWells& stdWells() const { return std_wells_; } + StandardWells& stdWells() { return std_wells_; } + const StandardWells& stdWells() const { return std_wells_; } /// return the Well struct in the StandardWells const Wells& wells() const { return std_wells_.wells(); } diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index e57f6dee0..3c9723df0 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -49,7 +49,6 @@ namespace Opm { Eigen::RowMajor> DataBlock; /// Class for handling the standard well model. - template class StandardWells { protected: struct WellOps { @@ -81,6 +80,7 @@ namespace Opm { Vector& wellPerforationPressureDiffs(); const Vector& wellPerforationPressureDiffs() const; + template void computePropertiesForWellConnectionPressures(const SolutionState& state, const WellState& xw, const BlackoilPropsAdInterface& fluid, diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 16478b5f2..df58cd376 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -28,8 +28,7 @@ namespace Opm { - template - StandardWells:: + StandardWells:: WellOps::WellOps(const Wells* wells) : w2p(), p2w(), @@ -67,9 +66,7 @@ namespace Opm - template - StandardWells:: - StandardWells(const Wells* wells_arg) + StandardWells::StandardWells(const Wells* wells_arg) : wells_(wells_arg) , wops_(wells_arg) , well_perforation_densities_(Vector()) @@ -81,10 +78,7 @@ namespace Opm - template - const Wells& - StandardWells:: - wells() const + const Wells& StandardWells::wells() const { assert(wells_ != 0); return *(wells_); @@ -94,10 +88,7 @@ namespace Opm - template - bool - StandardWells:: - wellsActive() const + bool StandardWells::wellsActive() const { return wells_active_; } @@ -106,10 +97,7 @@ namespace Opm - template - void - StandardWells:: - setWellsActive(const bool wells_active) + void StandardWells::setWellsActive(const bool wells_active) { wells_active_ = wells_active; } @@ -118,10 +106,7 @@ namespace Opm - template - bool - StandardWells:: - localWellsActive() const + bool StandardWells::localWellsActive() const { return wells_ ? (wells_->number_of_wells > 0 ) : false; } @@ -130,10 +115,8 @@ namespace Opm - template - const typename StandardWells::WellOps& - StandardWells:: - wellOps() const + const StandardWells::WellOps& + StandardWells::wellOps() const { return wops_; } @@ -142,10 +125,7 @@ namespace Opm - template - Vector& - StandardWells:: - wellPerforationDensities() + Vector& StandardWells::wellPerforationDensities() { return well_perforation_densities_; } @@ -154,10 +134,8 @@ namespace Opm - template const Vector& - StandardWells:: - wellPerforationDensities() const + StandardWells::wellPerforationDensities() const { return well_perforation_densities_; } @@ -166,10 +144,8 @@ namespace Opm - template Vector& - StandardWells:: - wellPerforationPressureDiffs() + StandardWells::wellPerforationPressureDiffs() { return well_perforation_pressure_diffs_; } @@ -178,10 +154,8 @@ namespace Opm - template const Vector& - StandardWells:: - wellPerforationPressureDiffs() const + StandardWells::wellPerforationPressureDiffs() const { return well_perforation_pressure_diffs_; } @@ -191,7 +165,7 @@ namespace Opm template void - StandardWells:: + StandardWells:: computePropertiesForWellConnectionPressures(const SolutionState& state, const WellState& xw, const BlackoilPropsAdInterface& fluid, From 5d99fac207180cc433ad33db4d2bf6e672188bed Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Apr 2016 12:50:51 +0200 Subject: [PATCH 04/13] adding StandardWellsSolvent for Solvent model. --- CMakeLists_files.cmake | 5 +- opm/autodiff/BlackoilModelBase.hpp | 7 - opm/autodiff/BlackoilModelBase_impl.hpp | 73 --------- opm/autodiff/BlackoilSolventModel.hpp | 15 +- opm/autodiff/BlackoilSolventModel_impl.hpp | 129 +-------------- opm/autodiff/StandardWellsSolvent.hpp | 61 +++++++ opm/autodiff/StandardWellsSolvent_impl.hpp | 181 +++++++++++++++++++++ 7 files changed, 253 insertions(+), 218 deletions(-) create mode 100644 opm/autodiff/StandardWellsSolvent.hpp create mode 100644 opm/autodiff/StandardWellsSolvent_impl.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 255fae7f0..fe95bcbc1 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -131,7 +131,7 @@ list (APPEND PROGRAM_SOURCE_FILES # originally generated with the command: # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort list (APPEND PUBLIC_HEADER_FILES -<<<<<<< HEAD + opm/autodiff/AdditionalObjectDeleter.hpp opm/autodiff/AutoDiffBlock.hpp opm/autodiff/AutoDiffHelpers.hpp @@ -197,7 +197,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/WellStateMultiSegment.hpp opm/autodiff/WellMultiSegment.hpp opm/autodiff/StandardWells.hpp - opm/autodiff/StandardWells_impl.hpp + opm/autodiff/StandardWellsSolvent.hpp opm/polymer/CompressibleTpfaPolymer.hpp opm/polymer/GravityColumnSolverPolymer.hpp opm/polymer/GravityColumnSolverPolymer_impl.hpp @@ -227,3 +227,4 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/SimulatorIncompTwophase.hpp ) + diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index b8ceae236..6e27f4c20 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -383,13 +383,6 @@ namespace Opm { void computeWellConnectionPressures(const SolutionState& state, const WellState& xw); - void computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf); - void assembleMassBalanceEq(const SolutionState& state); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 363b5a3d7..323a346bd 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -758,81 +758,8 @@ namespace detail { } } - template - void BlackoilModelBase::computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf) - { - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; - - // Compute the average pressure in each well block - const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); - V avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - avg_press[perf] = p_avg; - } - } - - const std::vector& well_cells = stdWells().wellOps().well_cells; - - // Use cell values for the temperature as the wells don't knows its temperature yet. - const ADB perf_temp = subset(state.temperature, well_cells); - - // Compute b, rsmax, rvmax values for perforations. - // Evaluate the properties using average well block pressures - // and cell values for rs, rv, phase condition and temperature. - const ADB avg_press_ad = ADB::constant(avg_press); - std::vector perf_cond(nperf); - const std::vector& pc = phaseCondition(); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = pc[well_cells[perf]]; - } - const PhaseUsage& pu = fluid_.phaseUsage(); - DataBlock b(nperf, pu.num_phases); - if (pu.phase_used[BlackoilPhases::Aqua]) { - const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; - } - assert(active_[Oil]); - const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); - if (pu.phase_used[BlackoilPhases::Liquid]) { - const ADB perf_rs = subset(state.rs, well_cells); - const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; - const V rssat = fluidRsSat(avg_press, perf_so, well_cells); - rsmax_perf.assign(rssat.data(), rssat.data() + nperf); - } else { - rsmax_perf.assign(nperf, 0.0); - } - if (pu.phase_used[BlackoilPhases::Vapour]) { - const ADB perf_rv = subset(state.rv, well_cells); - const V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; - const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); - rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); - } else { - rvmax_perf.assign(nperf, 0.0); - } - // b is row major, so can just copy data. - b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); - // Surface density. - // The compute density segment wants the surface densities as - // an np * number of wells cells array - V rho = superset(fluid_.surfaceDensity(0 , well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); - for (int phase = 1; phase < pu.num_phases; ++phase) { - rho += superset(fluid_.surfaceDensity(phase , well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - surf_dens_perf.assign(rho.data(), rho.data() + nperf * pu.num_phases); - } template diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 505f7a3a6..02b7dd3ac 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -25,6 +25,7 @@ #include #include #include +#include namespace Opm { @@ -103,6 +104,9 @@ namespace Opm { const bool is_miscible_; std::vector mu_eff_; std::vector b_eff_; + StandardWellsSolvent std_wells_; + const StandardWellsSolvent& stdWells() const { return std_wells_; } + StandardWellsSolvent& stdWells() { return std_wells_; } // Need to declare Base members we want to use here. @@ -130,7 +134,7 @@ namespace Opm { // --------- Protected methods --------- // Need to declare Base members we want to use here. - using Base::stdWells; + // using Base::stdWells; using Base::wells; using Base::variableState; using Base::computeGasPressure; @@ -148,7 +152,7 @@ namespace Opm { using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; - using Base::computePropertiesForWellConnectionPressures; + // using Base::computePropertiesForWellConnectionPressures; std::vector computeRelPerm(const SolutionState& state) const; @@ -202,13 +206,6 @@ namespace Opm { const SolutionState& state, WellState& xw); - void computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf); - void updateEquationsScaling(); void diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index c292c0626..6145267df 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -89,7 +89,8 @@ namespace Opm { has_solvent_(has_solvent), solvent_pos_(detail::solventPos(fluid.phaseUsage())), solvent_props_(solvent_props), - is_miscible_(is_miscible) + is_miscible_(is_miscible), + std_wells_(wells_arg, solvent_props) { if (has_solvent_) { @@ -381,132 +382,6 @@ namespace Opm { } } - template - void BlackoilSolventModel::computePropertiesForWellConnectionPressures(const SolutionState& state, - const WellState& xw, - std::vector& b_perf, - std::vector& rsmax_perf, - std::vector& rvmax_perf, - std::vector& surf_dens_perf) - { - using namespace Opm::AutoDiffGrid; - // 1. Compute properties required by computeConnectionPressureDelta(). - // Note that some of the complexity of this part is due to the function - // taking std::vector arguments, and not Eigen objects. - const int nperf = wells().well_connpos[wells().number_of_wells]; - const int nw = wells().number_of_wells; - const std::vector well_cells(wells().well_cells, wells().well_cells + nperf); - - // Compute the average pressure in each well block - const V perf_press = Eigen::Map(xw.perfPress().data(), nperf); - V avg_press = perf_press*0; - for (int w = 0; w < nw; ++w) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; - const double p_avg = (perf_press[perf] + p_above)/2; - avg_press[perf] = p_avg; - } - } - - // Use cell values for the temperature as the wells don't knows its temperature yet. - const ADB perf_temp = subset(state.temperature, well_cells); - - // Compute b, rsmax, rvmax values for perforations. - // Evaluate the properties using average well block pressures - // and cell values for rs, rv, phase condition and temperature. - const ADB avg_press_ad = ADB::constant(avg_press); - std::vector perf_cond(nperf); - const std::vector& pc = phaseCondition(); - for (int perf = 0; perf < nperf; ++perf) { - perf_cond[perf] = pc[well_cells[perf]]; - } - - const PhaseUsage& pu = fluid_.phaseUsage(); - DataBlock b(nperf, pu.num_phases); - - const V bw = fluid_.bWat(avg_press_ad, perf_temp, well_cells).value(); - if (pu.phase_used[BlackoilPhases::Aqua]) { - b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; - } - - assert(active_[Oil]); - assert(active_[Gas]); - const ADB perf_rv = subset(state.rv, well_cells); - const ADB perf_rs = subset(state.rs, well_cells); - const V perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); - if (pu.phase_used[BlackoilPhases::Liquid]) { - const V bo = fluid_.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); - //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value(); - b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; - const V rssat = fluidRsSat(avg_press, perf_so, well_cells); - rsmax_perf.assign(rssat.data(), rssat.data() + nperf); - } else { - rsmax_perf.assign(0.0, nperf); - } - V surf_dens_copy = superset(fluid_.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); - for (int phase = 1; phase < pu.num_phases; ++phase) { - if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) { - continue; // the gas surface density is added after the solvent is accounted for. - } - surf_dens_copy += superset(fluid_.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); - } - - if (pu.phase_used[BlackoilPhases::Vapour]) { - // Unclear wether the effective or the pure values should be used for the wells - // the current usage of unmodified properties values gives best match. - //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value(); - V bg = fluid_.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); - V rhog = fluid_.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells); - if (has_solvent_) { - - const V bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value(); - //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value(); - - // A weighted sum of the b-factors of gas and solvent are used. - const int nc = Opm::AutoDiffGrid::numCells(grid_); - - const ADB zero = ADB::constant(V::Zero(nc)); - const ADB& ss = state.solvent_saturation; - const ADB& sg = (active_[ Gas ] - ? state.saturation[ pu.phase_pos[ Gas ] ] - : zero); - - Selector zero_selector(ss.value() + sg.value(), Selector::Zero); - V F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value(); - - V injectedSolventFraction = Eigen::Map(&xw.solventFraction()[0], nperf); - - V isProducer = V::Zero(nperf); - V ones = V::Constant(nperf,1.0); - for (int w = 0; w < nw; ++w) { - if(wells().type[w] == PRODUCER) { - for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { - isProducer[perf] = 1; - } - } - } - - F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction; - - bg = bg * (ones - F_solvent); - bg = bg + F_solvent * bs; - - const V& rhos = solvent_props_.solventSurfaceDensity(well_cells); - rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos); - } - b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; - surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases); - - const V rvsat = fluidRvSat(avg_press, perf_so, well_cells); - rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); - } else { - rvmax_perf.assign(0.0, nperf); - } - - // b and surf_dens_perf is row major, so can just copy data. - b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); - surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases); - } diff --git a/opm/autodiff/StandardWellsSolvent.hpp b/opm/autodiff/StandardWellsSolvent.hpp new file mode 100644 index 000000000..9906eae3b --- /dev/null +++ b/opm/autodiff/StandardWellsSolvent.hpp @@ -0,0 +1,61 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED +#define OPM_STANDARDWELLSSOLVENT_HEADER_INCLUDED + +#include +#include + +namespace Opm { + + + /// Class for handling the standard well model for solvent model + class StandardWellsSolvent : public StandardWells + { + public: + + using Base = StandardWells; + + // --------- Public methods --------- + explicit StandardWellsSolvent(const Wells* wells, const SolventPropsAdFromDeck& solvent_props); + + template + void computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf); + protected: + const SolventPropsAdFromDeck& solvent_props_; + + }; + + +} // namespace Opm + +#include "StandardWellsSolvent_impl.hpp" + +#endif diff --git a/opm/autodiff/StandardWellsSolvent_impl.hpp b/opm/autodiff/StandardWellsSolvent_impl.hpp new file mode 100644 index 000000000..38b538675 --- /dev/null +++ b/opm/autodiff/StandardWellsSolvent_impl.hpp @@ -0,0 +1,181 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#include + + + + +namespace Opm +{ + + + + + + + StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, + const SolventPropsAdFromDeck& solvent_props) + : Base(wells_arg) + , solvent_props_(solvent_props) + { + } + + + + + + template + void + StandardWellsSolvent:: + computePropertiesForWellConnectionPressures(const SolutionState& state, + const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + const std::vector& pc, + std::vector& b_perf, + std::vector& rsmax_perf, + std::vector& rvmax_perf, + std::vector& surf_dens_perf) + { + // 1. Compute properties required by computeConnectionPressureDelta(). + // Note that some of the complexity of this part is due to the function + // taking std::vector arguments, and not Eigen objects. + const int nperf = wells().well_connpos[wells().number_of_wells]; + const int nw = wells().number_of_wells; + + // Compute the average pressure in each well block + const Vector perf_press = Eigen::Map(xw.perfPress().data(), nperf); + Vector avg_press = perf_press*0; + for (int w = 0; w < nw; ++w) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + const double p_above = perf == wells().well_connpos[w] ? state.bhp.value()[w] : perf_press[perf - 1]; + const double p_avg = (perf_press[perf] + p_above)/2; + avg_press[perf] = p_avg; + } + } + + const std::vector& well_cells = wellOps().well_cells; + + // Use cell values for the temperature as the wells don't knows its temperature yet. + const ADB perf_temp = subset(state.temperature, well_cells); + + // Compute b, rsmax, rvmax values for perforations. + // Evaluate the properties using average well block pressures + // and cell values for rs, rv, phase condition and temperature. + const ADB avg_press_ad = ADB::constant(avg_press); + std::vector perf_cond(nperf); + for (int perf = 0; perf < nperf; ++perf) { + perf_cond[perf] = pc[well_cells[perf]]; + } + + const PhaseUsage& pu = fluid.phaseUsage(); + DataBlock b(nperf, pu.num_phases); + + const Vector bw = fluid.bWat(avg_press_ad, perf_temp, well_cells).value(); + if (pu.phase_used[BlackoilPhases::Aqua]) { + b.col(pu.phase_pos[BlackoilPhases::Aqua]) = bw; + } + + assert(active[Oil]); + assert(active[Gas]); + const ADB perf_rv = subset(state.rv, well_cells); + const ADB perf_rs = subset(state.rs, well_cells); + const Vector perf_so = subset(state.saturation[pu.phase_pos[Oil]].value(), well_cells); + if (pu.phase_used[BlackoilPhases::Liquid]) { + const Vector bo = fluid.bOil(avg_press_ad, perf_temp, perf_rs, perf_cond, well_cells).value(); + //const V bo_eff = subset(rq_[pu.phase_pos[Oil] ].b , well_cells).value(); + b.col(pu.phase_pos[BlackoilPhases::Liquid]) = bo; + // const Vector rssat = fluidRsSat(avg_press, perf_so, well_cells); + const Vector rssat = fluid.rsSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rsmax_perf.assign(rssat.data(), rssat.data() + nperf); + } else { + rsmax_perf.assign(0.0, nperf); + } + V surf_dens_copy = superset(fluid.surfaceDensity(0, well_cells), Span(nperf, pu.num_phases, 0), nperf*pu.num_phases); + for (int phase = 1; phase < pu.num_phases; ++phase) { + if ( phase == pu.phase_pos[BlackoilPhases::Vapour]) { + continue; // the gas surface density is added after the solvent is accounted for. + } + surf_dens_copy += superset(fluid.surfaceDensity(phase, well_cells), Span(nperf, pu.num_phases, phase), nperf*pu.num_phases); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + // Unclear wether the effective or the pure values should be used for the wells + // the current usage of unmodified properties values gives best match. + //V bg_eff = subset(rq_[pu.phase_pos[Gas]].b,well_cells).value(); + Vector bg = fluid.bGas(avg_press_ad, perf_temp, perf_rv, perf_cond, well_cells).value(); + Vector rhog = fluid.surfaceDensity(pu.phase_pos[BlackoilPhases::Vapour], well_cells); + // to handle solvent related + { + + const Vector bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value(); + //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value(); + + // number of cells + const int nc = state.pressure.size(); + + const ADB zero = ADB::constant(Vector::Zero(nc)); + const ADB& ss = state.solvent_saturation; + const ADB& sg = (active[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : zero); + + Selector zero_selector(ss.value() + sg.value(), Selector::Zero); + Vector F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells).value(); + + Vector injectedSolventFraction = Eigen::Map(&xw.solventFraction()[0], nperf); + + Vector isProducer = Vector::Zero(nperf); + Vector ones = Vector::Constant(nperf,1.0); + for (int w = 0; w < nw; ++w) { + if(wells().type[w] == PRODUCER) { + for (int perf = wells().well_connpos[w]; perf < wells().well_connpos[w+1]; ++perf) { + isProducer[perf] = 1; + } + } + } + + F_solvent = isProducer * F_solvent + (ones - isProducer) * injectedSolventFraction; + + bg = bg * (ones - F_solvent); + bg = bg + F_solvent * bs; + + const Vector& rhos = solvent_props_.solventSurfaceDensity(well_cells); + rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos); + } + b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg; + surf_dens_copy += superset(rhog, Span(nperf, pu.num_phases, pu.phase_pos[BlackoilPhases::Vapour]), nperf*pu.num_phases); + + // const Vector rvsat = fluidRvSat(avg_press, perf_so, well_cells); + const Vector rvsat = fluid.rvSat(ADB::constant(avg_press), ADB::constant(perf_so), well_cells).value(); + rvmax_perf.assign(rvsat.data(), rvsat.data() + nperf); + } else { + rvmax_perf.assign(0.0, nperf); + } + + // b and surf_dens_perf is row major, so can just copy data. + b_perf.assign(b.data(), b.data() + nperf * pu.num_phases); + surf_dens_perf.assign(surf_dens_copy.data(), surf_dens_copy.data() + nperf * pu.num_phases); + } + + +} From c8b66821d516705ce81a4643c311e7c7f87cd893 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Apr 2016 16:05:56 +0200 Subject: [PATCH 05/13] adding computeWellConnectionDensitesPressures to StandardWells class. --- opm/autodiff/BlackoilModelBase_impl.hpp | 42 +++++++++---------------- opm/autodiff/StandardWells.hpp | 10 ++++++ opm/autodiff/StandardWells_impl.hpp | 30 ++++++++++++++++++ 3 files changed, 54 insertions(+), 28 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 323a346bd..43bc57e59 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -778,31 +778,17 @@ namespace detail { std::vector surf_dens_perf; asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - // 2. Compute densities - std::vector cd = - WellDensitySegmented::computeConnectionDensities( - wells(), xw, fluid_.phaseUsage(), - b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); - - const int nperf = wells().well_connpos[wells().number_of_wells]; - const std::vector& well_cells = stdWells().wellOps().well_cells; - // Extract well connection depths. - const V depth = cellCentroidsZToEigen(grid_); - const V pdepth = subset(depth, well_cells); - std::vector perf_depth(pdepth.data(), pdepth.data() + nperf); + const Vector depth = cellCentroidsZToEigen(grid_); + const Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells); + const int nperf = wells().well_connpos[wells().number_of_wells]; + const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); // Gravity - double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); + const double grav = detail::getGravity(geo_.gravity(), dimensions(grid_)); - // 3. Compute pressure deltas - std::vector cdp = - WellDensitySegmented::computeConnectionPressureDelta( - wells(), perf_depth, cd, grav); + asImpl().stdWells().computeWellConnectionDensitesPressures(xw, fluid_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf, depth_perf, grav); - // 4. Store the results - stdWells().wellPerforationDensities() = Eigen::Map(cd.data(), nperf); - stdWells().wellPerforationPressureDiffs() = Eigen::Map(cdp.data(), nperf); } @@ -1057,7 +1043,7 @@ namespace detail { const std::vector& well_cells = stdWells().wellOps().well_cells; // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = stdWells().wellPerforationPressureDiffs(); + const V& cdp = asImpl().stdWells().wellPerforationPressureDiffs(); // Extract needed quantities for the perforation cells const ADB& p_perfcells = subset(state.pressure, well_cells); const ADB& rv_perfcells = subset(state.rv, well_cells); @@ -1213,7 +1199,7 @@ namespace detail { xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); // Update the perforation pressures. - const V& cdp = stdWells().wellPerforationPressureDiffs(); + const V& cdp = asImpl().stdWells().wellPerforationPressureDiffs(); const V perfpressure = (stdWells().wellOps().w2p * state.bhp.value().matrix()).array() + cdp; xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); } @@ -1498,14 +1484,14 @@ namespace detail { if (well_type == INJECTOR) { double dp = detail::computeHydrostaticCorrection( wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); + asImpl().stdWells().wellPerforationDensities(), gravity); xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; } else if (well_type == PRODUCER) { double dp = detail::computeHydrostaticCorrection( wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); + asImpl().stdWells().wellPerforationDensities(), gravity); xw.bhp()[w] = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; } @@ -1737,7 +1723,7 @@ namespace detail { case THP: { const int perf = wells().well_connpos[w]; - rho_v[w] = stdWells().wellPerforationDensities()[perf]; + rho_v[w] = asImpl().stdWells().wellPerforationDensities()[perf]; const int table_id = well_controls_iget_vfp(wc, current); const double target = well_controls_iget_target(wc, current); @@ -1800,7 +1786,7 @@ namespace detail { //Perform hydrostatic correction to computed targets double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - const ADB::V dp_v = detail::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, stdWells().wellPerforationDensities(), gravity); + const ADB::V dp_v = detail::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, asImpl().stdWells().wellPerforationDensities(), gravity); const ADB dp = ADB::constant(dp_v); const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw); @@ -2241,14 +2227,14 @@ namespace detail { if (well_type == INJECTOR) { double dp = detail::computeHydrostaticCorrection( wells(), w, vfp_properties_.getInj()->getTable(table_id)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); + asImpl().stdWells().wellPerforationDensities(), gravity); well_state.thp()[w] = vfp_properties_.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp); } else if (well_type == PRODUCER) { double dp = detail::computeHydrostaticCorrection( wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(), - stdWells().wellPerforationDensities(), gravity); + asImpl().stdWells().wellPerforationDensities(), gravity); well_state.thp()[w] = vfp_properties_.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); } diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 3c9723df0..85222b7bc 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -91,6 +91,16 @@ namespace Opm { std::vector& rvmax_perf, std::vector& surf_dens_perf); + template + void computeWellConnectionDensitesPressures(const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav); + protected: bool wells_active_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index df58cd376..dd5034172 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -20,6 +20,7 @@ #include +#include @@ -246,5 +247,34 @@ namespace Opm } + template + void + StandardWells:: + computeWellConnectionDensitesPressures(const WellState& xw, + const BlackoilPropsAdInterface& fluid, + const std::vector& b_perf, + const std::vector& rsmax_perf, + const std::vector& rvmax_perf, + const std::vector& surf_dens_perf, + const std::vector& depth_perf, + const double grav) + { + // Compute densities + std::vector cd = + WellDensitySegmented::computeConnectionDensities( + wells(), xw, fluid.phaseUsage(), + b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); + + const int nperf = wells().well_connpos[wells().number_of_wells]; + + // Compute pressure deltas + std::vector cdp = + WellDensitySegmented::computeConnectionPressureDelta( + wells(), depth_perf, cd, grav); + + // Store the results + well_perforation_densities_ = Eigen::Map(cd.data(), nperf); + well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); + } } From e9e1b9fda829ff02bced74e93c57d718ec952bad Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Apr 2016 16:48:31 +0200 Subject: [PATCH 06/13] adding extractWellPerfProperties to StandardWells It causes problem for the flow_multisegment. So the version in BlackoilModelBase is kept for now. On the other hand, it is few of functions that will be both required by the standard wells and multisegment wells. Some decision will be made later on how to put this function. --- opm/autodiff/BlackoilModelBase_impl.hpp | 2 +- opm/autodiff/StandardWells.hpp | 6 +++++ opm/autodiff/StandardWells_impl.hpp | 33 +++++++++++++++++++++++++ 3 files changed, 40 insertions(+), 1 deletion(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 43bc57e59..44a9c544d 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -851,7 +851,7 @@ namespace detail { std::vector mob_perfcells; std::vector b_perfcells; - asImpl().extractWellPerfProperties(state, mob_perfcells, b_perfcells); + asImpl().stdWells().extractWellPerfProperties(rq_, fluid_.numPhases(), mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 85222b7bc..3d2d85e8e 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -101,6 +101,12 @@ namespace Opm { const std::vector& depth_perf, const double grav); + template + void extractWellPerfProperties(const std::vector& rq, + const int np, + std::vector& mob_perfcells, + std::vector& b_perfcells) const; + protected: bool wells_active_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index dd5034172..55de376f8 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -247,6 +247,10 @@ namespace Opm } + + + + template void StandardWells:: @@ -277,4 +281,33 @@ namespace Opm well_perforation_pressure_diffs_ = Eigen::Map(cdp.data(), nperf); } + + + + + template + void + StandardWells:: + extractWellPerfProperties(const std::vector& rq, + const int np, + std::vector& mob_perfcells, + std::vector& b_perfcells) const + { + // If we have wells, extract the mobilities and b-factors for + // the well-perforated cells. + if ( !localWellsActive() ) { + mob_perfcells.clear(); + b_perfcells.clear(); + return; + } else { + const std::vector& well_cells = wellOps().well_cells; + mob_perfcells.resize(np, ADB::null()); + b_perfcells.resize(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + mob_perfcells[phase] = subset(rq[phase].mob, well_cells); + b_perfcells[phase] = subset(rq[phase].b, well_cells); + } + } + } + } From c398a6e424b7d78f936adadfeac43683feaf6f92 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Fri, 8 Apr 2016 17:26:07 +0200 Subject: [PATCH 07/13] puting computeWellFlux in StandardWells --- opm/autodiff/BlackoilModelBase.hpp | 7 - opm/autodiff/BlackoilModelBase_impl.hpp | 156 +----------------- opm/autodiff/StandardWells.hpp | 10 ++ opm/autodiff/StandardWells_impl.hpp | 154 +++++++++++++++++ .../BlackoilPolymerModel_impl.hpp | 4 +- 5 files changed, 168 insertions(+), 163 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 6e27f4c20..0e8784a81 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -397,13 +397,6 @@ namespace Opm { SolutionState& state, WellState& well_state); - void - computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s) const; - void updatePerfPhaseRatesAndPressures(const std::vector& cq_s, const SolutionState& state, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 44a9c544d..079475ee0 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -858,7 +858,7 @@ namespace detail { } V aliveWells; std::vector cq_s; - asImpl().computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + asImpl().stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); asImpl().addWellFluxEq(cq_s, state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); @@ -1024,158 +1024,6 @@ namespace detail { - - template - void - BlackoilModelBase::computeWellFlux(const SolutionState& state, - const std::vector& mob_perfcells, - const std::vector& b_perfcells, - V& aliveWells, - std::vector& cq_s) const - { - if( ! localWellsActive() ) return ; - - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - V Tw = Eigen::Map(wells().WI, nperf); - const std::vector& well_cells = stdWells().wellOps().well_cells; - - // pressure diffs computed already (once per step, not changing per iteration) - const V& cdp = asImpl().stdWells().wellPerforationPressureDiffs(); - // Extract needed quantities for the perforation cells - const ADB& p_perfcells = subset(state.pressure, well_cells); - const ADB& rv_perfcells = subset(state.rv, well_cells); - const ADB& rs_perfcells = subset(state.rs, well_cells); - - // Perforation pressure - const ADB perfpressure = (stdWells().wellOps().w2p * state.bhp) + cdp; - - // Pressure drawdown (also used to determine direction of flow) - const ADB drawdown = p_perfcells - perfpressure; - - // Compute vectors with zero and ones that - // selects the wanted quantities. - - // selects injection perforations - V selectInjectingPerforations = V::Zero(nperf); - // selects producing perforations - V selectProducingPerforations = V::Zero(nperf); - for (int c = 0; c < nperf; ++c){ - if (drawdown.value()[c] < 0) - selectInjectingPerforations[c] = 1; - else - selectProducingPerforations[c] = 1; - } - - // Handle cross flow - const V numInjectingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); - const V numProducingPerforations = (stdWells().wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); - for (int w = 0; w < nw; ++w) { - if (!wells().allow_cf[w]) { - for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { - // Crossflow is not allowed; reverse flow is prevented. - // At least one of the perforation must be open in order to have a meeningful - // equation to solve. For the special case where all perforations have reverse flow, - // and the target rate is non-zero all of the perforations are keept open. - if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { - selectProducingPerforations[perf] = 0.0; - } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ - selectInjectingPerforations[perf] = 0.0; - } - } - } - } - - // HANDLE FLOW INTO WELLBORE - // compute phase volumetric rates at standard conditions - std::vector cq_ps(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); - cq_ps[phase] = b_perfcells[phase] * cq_p; - } - if (active_[Oil] && active_[Gas]) { - const int oilpos = pu.phase_pos[Oil]; - const int gaspos = pu.phase_pos[Gas]; - const ADB cq_psOil = cq_ps[oilpos]; - const ADB cq_psGas = cq_ps[gaspos]; - cq_ps[gaspos] += rs_perfcells * cq_psOil; - cq_ps[oilpos] += rv_perfcells * cq_psGas; - } - - // HANDLE FLOW OUT FROM WELLBORE - // Using total mobilities - ADB total_mob = mob_perfcells[0]; - for (int phase = 1; phase < np; ++phase) { - total_mob += mob_perfcells[phase]; - } - // injection perforations total volume rates - const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); - - // compute wellbore mixture for injecting perforations - // The wellbore mixture depends on the inflow from the reservoar - // and the well injection rates. - - // compute avg. and total wellbore phase volumetric rates at standard conds - const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); - std::vector wbq(np, ADB::null()); - ADB wbqt = ADB::constant(V::Zero(nw)); - for (int phase = 0; phase < np; ++phase) { - const ADB& q_ps = stdWells().wellOps().p2w * cq_ps[phase]; - const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); - Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); - const int pos = pu.phase_pos[phase]; - wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(V::Zero(nw)))) - q_ps; - wbqt += wbq[phase]; - } - // compute wellbore mixture at standard conditions. - Selector notDeadWells_selector(wbqt.value(), Selector::Zero); - std::vector cmix_s(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - const int pos = pu.phase_pos[phase]; - cmix_s[phase] = stdWells().wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); - } - - // compute volume ratio between connection at standard conditions - ADB volumeRatio = ADB::constant(V::Zero(nperf)); - const ADB d = V::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; - for (int phase = 0; phase < np; ++phase) { - ADB tmp = cmix_s[phase]; - - if (phase == Oil && active_[Gas]) { - const int gaspos = pu.phase_pos[Gas]; - tmp -= rv_perfcells * cmix_s[gaspos] / d; - } - if (phase == Gas && active_[Oil]) { - const int oilpos = pu.phase_pos[Oil]; - tmp -= rs_perfcells * cmix_s[oilpos] / d; - } - volumeRatio += tmp / b_perfcells[phase]; - } - - // injecting connections total volumerates at standard conditions - ADB cqt_is = cqt_i/volumeRatio; - - // connection phase volumerates at standard conditions - cq_s.resize(np, ADB::null()); - for (int phase = 0; phase < np; ++phase) { - cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; - } - - // check for dead wells (used in the well controll equations) - aliveWells = V::Constant(nw, 1.0); - for (int w = 0; w < nw; ++w) { - if (wbqt.value()[w] == 0) { - aliveWells[w] = 0.0; - } - } - } - - - - - template void BlackoilModelBase::updatePerfPhaseRatesAndPressures(const std::vector& cq_s, const SolutionState& state, @@ -1587,7 +1435,7 @@ namespace detail { SolutionState wellSolutionState = state0; asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState); - asImpl().computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); + asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); asImpl().addWellFluxEq(cq_s, wellSolutionState); asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 3d2d85e8e..632ed183d 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -107,6 +107,16 @@ namespace Opm { std::vector& mob_perfcells, std::vector& b_perfcells) const; + template + void computeWellFlux(const SolutionState& state, + const Opm::PhaseUsage& phase_usage, + const std::vector& active, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + std::vector& cq_s) const; + + protected: bool wells_active_; diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 55de376f8..c154bf130 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -310,4 +310,158 @@ namespace Opm } } + + + + + + template + void + StandardWells:: + computeWellFlux(const SolutionState& state, + const Opm::PhaseUsage& pu, + const std::vector& active, + const std::vector& mob_perfcells, + const std::vector& b_perfcells, + Vector& aliveWells, + std::vector& cq_s) const + { + if( ! localWellsActive() ) return ; + + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + Vector Tw = Eigen::Map(wells().WI, nperf); + const std::vector& well_cells = wellOps().well_cells; + + // pressure diffs computed already (once per step, not changing per iteration) + const Vector& cdp = wellPerforationPressureDiffs(); + // Extract needed quantities for the perforation cells + const ADB& p_perfcells = subset(state.pressure, well_cells); + const ADB& rv_perfcells = subset(state.rv, well_cells); + const ADB& rs_perfcells = subset(state.rs, well_cells); + + // Perforation pressure + const ADB perfpressure = (wellOps().w2p * state.bhp) + cdp; + + // Pressure drawdown (also used to determine direction of flow) + const ADB drawdown = p_perfcells - perfpressure; + + // Compute vectors with zero and ones that + // selects the wanted quantities. + + // selects injection perforations + Vector selectInjectingPerforations = Vector::Zero(nperf); + // selects producing perforations + Vector selectProducingPerforations = Vector::Zero(nperf); + for (int c = 0; c < nperf; ++c){ + if (drawdown.value()[c] < 0) + selectInjectingPerforations[c] = 1; + else + selectProducingPerforations[c] = 1; + } + + // Handle cross flow + const Vector numInjectingPerforations = (wellOps().p2w * ADB::constant(selectInjectingPerforations)).value(); + const Vector numProducingPerforations = (wellOps().p2w * ADB::constant(selectProducingPerforations)).value(); + for (int w = 0; w < nw; ++w) { + if (!wells().allow_cf[w]) { + for (int perf = wells().well_connpos[w] ; perf < wells().well_connpos[w+1]; ++perf) { + // Crossflow is not allowed; reverse flow is prevented. + // At least one of the perforation must be open in order to have a meeningful + // equation to solve. For the special case where all perforations have reverse flow, + // and the target rate is non-zero all of the perforations are keept open. + if (wells().type[w] == INJECTOR && numInjectingPerforations[w] > 0) { + selectProducingPerforations[perf] = 0.0; + } else if (wells().type[w] == PRODUCER && numProducingPerforations[w] > 0 ){ + selectInjectingPerforations[perf] = 0.0; + } + } + } + } + + // HANDLE FLOW INTO WELLBORE + // compute phase volumetric rates at standard conditions + std::vector cq_ps(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const ADB cq_p = -(selectProducingPerforations * Tw) * (mob_perfcells[phase] * drawdown); + cq_ps[phase] = b_perfcells[phase] * cq_p; + } + if (active[Oil] && active[Gas]) { + const int oilpos = pu.phase_pos[Oil]; + const int gaspos = pu.phase_pos[Gas]; + const ADB cq_psOil = cq_ps[oilpos]; + const ADB cq_psGas = cq_ps[gaspos]; + cq_ps[gaspos] += rs_perfcells * cq_psOil; + cq_ps[oilpos] += rv_perfcells * cq_psGas; + } + + // HANDLE FLOW OUT FROM WELLBORE + // Using total mobilities + ADB total_mob = mob_perfcells[0]; + for (int phase = 1; phase < np; ++phase) { + total_mob += mob_perfcells[phase]; + } + // injection perforations total volume rates + const ADB cqt_i = -(selectInjectingPerforations * Tw) * (total_mob * drawdown); + + // compute wellbore mixture for injecting perforations + // The wellbore mixture depends on the inflow from the reservoar + // and the well injection rates. + + // compute avg. and total wellbore phase volumetric rates at standard conds + const DataBlock compi = Eigen::Map(wells().comp_frac, nw, np); + std::vector wbq(np, ADB::null()); + ADB wbqt = ADB::constant(Vector::Zero(nw)); + for (int phase = 0; phase < np; ++phase) { + const ADB& q_ps = wellOps().p2w * cq_ps[phase]; + const ADB& q_s = subset(state.qs, Span(nw, 1, phase*nw)); + Selector injectingPhase_selector(q_s.value(), Selector::GreaterZero); + const int pos = pu.phase_pos[phase]; + wbq[phase] = (compi.col(pos) * injectingPhase_selector.select(q_s,ADB::constant(Vector::Zero(nw)))) - q_ps; + wbqt += wbq[phase]; + } + // compute wellbore mixture at standard conditions. + Selector notDeadWells_selector(wbqt.value(), Selector::Zero); + std::vector cmix_s(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + const int pos = pu.phase_pos[phase]; + cmix_s[phase] = wellOps().w2p * notDeadWells_selector.select(ADB::constant(compi.col(pos)), wbq[phase]/wbqt); + } + + // compute volume ratio between connection at standard conditions + ADB volumeRatio = ADB::constant(Vector::Zero(nperf)); + const ADB d = Vector::Constant(nperf,1.0) - rv_perfcells * rs_perfcells; + for (int phase = 0; phase < np; ++phase) { + ADB tmp = cmix_s[phase]; + + if (phase == Oil && active[Gas]) { + const int gaspos = pu.phase_pos[Gas]; + tmp -= rv_perfcells * cmix_s[gaspos] / d; + } + if (phase == Gas && active[Oil]) { + const int oilpos = pu.phase_pos[Oil]; + tmp -= rs_perfcells * cmix_s[oilpos] / d; + } + volumeRatio += tmp / b_perfcells[phase]; + } + + // injecting connections total volumerates at standard conditions + ADB cqt_is = cqt_i/volumeRatio; + + // connection phase volumerates at standard conditions + cq_s.resize(np, ADB::null()); + for (int phase = 0; phase < np; ++phase) { + cq_s[phase] = cq_ps[phase] + cmix_s[phase]*cqt_is; + } + + // check for dead wells (used in the well controll equations) + aliveWells = Vector::Constant(nw, 1.0); + for (int w = 0; w < nw; ++w) { + if (wbqt.value()[w] == 0) { + aliveWells[w] = 0.0; + } + } + } + } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index c50cf38d8..53a81315e 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -550,7 +550,7 @@ namespace Opm { Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); } - Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); if (has_plyshlog_) { std::vector water_vel_wells; @@ -569,7 +569,7 @@ namespace Opm { mob_perfcells[water_pos] = mob_perfcells[water_pos] / shear_mult_wells_adb; } - Base::computeWellFlux(state, mob_perfcells, b_perfcells, aliveWells, cq_s); + stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); Base::addWellFluxEq(cq_s, state); addWellContributionToMassBalanceEq(cq_s, state, well_state); From 376c940f0988dfa1610830994051cfb27e20275b Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 11 Apr 2016 10:46:21 +0200 Subject: [PATCH 08/13] moving updatePerfPhaseRatesAndPressures to StandardWells --- opm/autodiff/BlackoilModelBase.hpp | 5 --- opm/autodiff/BlackoilModelBase_impl.hpp | 36 ++----------------- opm/autodiff/StandardWells.hpp | 5 +++ opm/autodiff/StandardWells_impl.hpp | 34 ++++++++++++++++++ .../BlackoilPolymerModel_impl.hpp | 2 +- 5 files changed, 42 insertions(+), 40 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 0e8784a81..f9cf5d8fb 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -397,11 +397,6 @@ namespace Opm { SolutionState& state, WellState& well_state); - void - updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - const SolutionState& state, - WellState& xw) const; - void addWellFluxEq(const std::vector& cq_s, const SolutionState& state); diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 079475ee0..ae9dca1cd 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -859,7 +859,7 @@ namespace detail { V aliveWells; std::vector cq_s; asImpl().stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); - asImpl().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); + asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); asImpl().addWellFluxEq(cq_s, state); asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state); asImpl().addWellControlEq(state, well_state, aliveWells); @@ -1024,38 +1024,6 @@ namespace detail { - template - void BlackoilModelBase::updatePerfPhaseRatesAndPressures(const std::vector& cq_s, - const SolutionState& state, - WellState& xw) const - { - if ( !asImpl().localWellsActive() ) - { - // If there are no wells in the subdomain of the proces then - // cq_s has zero size and will cause a segmentation fault below. - return; - } - - // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const int nperf = wells().well_connpos[nw]; - V cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); - for (int phase = 1; phase < np; ++phase) { - cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); - } - xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); - - // Update the perforation pressures. - const V& cdp = asImpl().stdWells().wellPerforationPressureDiffs(); - const V perfpressure = (stdWells().wellOps().w2p * state.bhp.value().matrix()).array() + cdp; - xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); - } - - - - - template void BlackoilModelBase::addWellFluxEq(const std::vector& cq_s, const SolutionState& state) @@ -1436,7 +1404,7 @@ namespace detail { SolutionState wellSolutionState = state0; asImpl().variableStateExtractWellsVars(indices, vars, wellSolutionState); asImpl().stdWells().computeWellFlux(wellSolutionState, fluid_.phaseUsage(), active_, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); - asImpl().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); + asImpl().stdWells().updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); asImpl().addWellFluxEq(cq_s, wellSolutionState); asImpl().addWellControlEq(wellSolutionState, well_state, aliveWells); converged = getWellConvergence(it); diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 632ed183d..005cd2fe3 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -116,6 +116,11 @@ namespace Opm { Vector& aliveWells, std::vector& cq_s) const; + template + void updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const; + protected: diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index c154bf130..1794478a3 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -464,4 +464,38 @@ namespace Opm } } + + + + + template + void + StandardWells:: + updatePerfPhaseRatesAndPressures(const std::vector& cq_s, + const SolutionState& state, + WellState& xw) const + { + if ( !localWellsActive() ) + { + // If there are no wells in the subdomain of the proces then + // cq_s has zero size and will cause a segmentation fault below. + return; + } + + // Update the perforation phase rates (used to calculate the pressure drop in the wellbore). + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + const int nperf = wells().well_connpos[nw]; + Vector cq = superset(cq_s[0].value(), Span(nperf, np, 0), nperf*np); + for (int phase = 1; phase < np; ++phase) { + cq += superset(cq_s[phase].value(), Span(nperf, np, phase), nperf*np); + } + xw.perfPhaseRates().assign(cq.data(), cq.data() + nperf*np); + + // Update the perforation pressures. + const Vector& cdp = wellPerforationPressureDiffs(); + const Vector perfpressure = (wellOps().w2p * state.bhp.value().matrix()).array() + cdp; + xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); + } + } diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp index 53a81315e..03c1b47ee 100644 --- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp +++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp @@ -570,7 +570,7 @@ namespace Opm { } stdWells().computeWellFlux(state, fluid_.phaseUsage(), active_, mob_perfcells, b_perfcells, aliveWells, cq_s); - Base::updatePerfPhaseRatesAndPressures(cq_s, state, well_state); + stdWells().updatePerfPhaseRatesAndPressures(cq_s, state, well_state); Base::addWellFluxEq(cq_s, state); addWellContributionToMassBalanceEq(cq_s, state, well_state); addWellControlEq(state, well_state, aliveWells); From 5da57973fe1c7c1c938aac9f2bdfae4c17c5aa44 Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Mon, 11 Apr 2016 15:45:03 +0200 Subject: [PATCH 09/13] adding updateWellState and updateWellControls to StandardWells --- CMakeLists_files.cmake | 1 + opm/autodiff/BlackoilModelBase.hpp | 5 - opm/autodiff/BlackoilModelBase_impl.hpp | 401 +----------------- .../BlackoilMultiSegmentModel_impl.hpp | 2 +- opm/autodiff/BlackoilSolventModel.hpp | 2 +- opm/autodiff/StandardWells.hpp | 16 + opm/autodiff/StandardWells_impl.hpp | 259 +++++++++++ opm/autodiff/WellHelpers.hpp | 169 ++++++++ .../fullyimplicit/BlackoilPolymerModel.hpp | 3 +- .../BlackoilPolymerModel_impl.hpp | 4 +- 10 files changed, 466 insertions(+), 396 deletions(-) create mode 100644 opm/autodiff/WellHelpers.hpp diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index fe95bcbc1..bd7868398 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -196,6 +196,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/autodiff/VFPInjProperties.hpp opm/autodiff/WellStateMultiSegment.hpp opm/autodiff/WellMultiSegment.hpp + opm/autodiff/WellHelpers.hpp opm/autodiff/StandardWells.hpp opm/autodiff/StandardWellsSolvent.hpp opm/polymer/CompressibleTpfaPolymer.hpp diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index f9cf5d8fb..ccdbb45f2 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -411,11 +411,6 @@ namespace Opm { const WellState& xw, const V& aliveWells); - void updateWellControls(WellState& xw) const; - - void updateWellState(const V& dwells, - WellState& well_state); - bool getWellConvergence(const int iteration); bool isVFPActive() const; diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index ae9dca1cd..f19e33077 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include #include #include @@ -816,7 +817,10 @@ namespace detail { // Possibly switch well controls and updating well state to // get reasonable initial conditions for the wells - asImpl().updateWellControls(well_state); + // asImpl().updateWellControls(well_state); + // asImpl().stdWells().updateWellControls(well_state); + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + asImpl().stdWells().updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); // Create the primary variables. SolutionState state = asImpl().variableState(reservoir_state, well_state); @@ -1050,138 +1054,6 @@ namespace detail { - namespace detail - { - inline - double rateToCompare(const std::vector& well_phase_flow_rate, - const int well, - const int num_phases, - const double* distr) - { - double rate = 0.0; - for (int phase = 0; phase < num_phases; ++phase) { - // Important: well_phase_flow_rate is ordered with all phase rates for first - // well first, then all phase rates for second well etc. - rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase]; - } - return rate; - } - - inline - bool constraintBroken(const std::vector& bhp, - const std::vector& thp, - const std::vector& well_phase_flow_rate, - const int well, - const int num_phases, - const WellType& well_type, - const WellControls* wc, - const int ctrl_index) - { - const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index); - const double target = well_controls_iget_target(wc, ctrl_index); - const double* distr = well_controls_iget_distr(wc, ctrl_index); - - bool broken = false; - - switch (well_type) { - case INJECTOR: - { - switch (ctrl_type) { - case BHP: - broken = bhp[well] > target; - break; - - case THP: - broken = thp[well] > target; - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - broken = rateToCompare(well_phase_flow_rate, - well, num_phases, distr) > target; - break; - } - } - break; - - case PRODUCER: - { - switch (ctrl_type) { - case BHP: - broken = bhp[well] < target; - break; - - case THP: - broken = thp[well] < target; - break; - - case RESERVOIR_RATE: // Intentional fall-through - case SURFACE_RATE: - // Note that the rates compared below are negative, - // so breaking the constraints means: too high flow rate - // (as for injection). - broken = rateToCompare(well_phase_flow_rate, - well, num_phases, distr) < target; - break; - } - } - break; - - default: - OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells."); - } - - return broken; - } - } // namespace detail - - - namespace detail { - /** - * Simple hydrostatic correction for VFP table - * @param wells - wells struct - * @param w Well number - * @param vfp_table VFP table - * @param well_perforation_densities Densities at well perforations - * @param gravity Gravitational constant (e.g., 9.81...) - */ - inline - double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth, - const ADB::V& well_perforation_densities, const double gravity) { - if ( wells.well_connpos[w] == wells.well_connpos[w+1] ) - { - // This is a well with no perforations. - // If this is the last well we would subscript over the - // bounds below. - // we assume well_perforation_densities to be 0 - return 0; - } - const double well_ref_depth = wells.depth_ref[w]; - const double dh = vfp_ref_depth - well_ref_depth; - const int perf = wells.well_connpos[w]; - const double rho = well_perforation_densities[perf]; - const double dp = rho*gravity*dh; - - return dp; - } - - inline - ADB::V computeHydrostaticCorrection(const Wells& wells, const ADB::V vfp_ref_depth, - const ADB::V& well_perforation_densities, const double gravity) { - const int nw = wells.number_of_wells; - ADB::V retval = ADB::V::Zero(nw); - -#pragma omp parallel for schedule(static) - for (int i=0; i bool BlackoilModelBase::isVFPActive() const { @@ -1214,156 +1086,6 @@ namespace detail { } - template - void BlackoilModelBase::updateWellControls(WellState& xw) const - { - if( ! localWellsActive() ) return ; - - std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; - // Find, for each well, if any constraints are broken. If so, - // switch control to first broken constraint. - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); -#pragma omp parallel for schedule(dynamic) - for (int w = 0; w < nw; ++w) { - const WellControls* wc = wells().ctrls[w]; - // The current control in the well state overrides - // the current control set in the Wells struct, which - // is instead treated as a default. - int current = xw.currentControls()[w]; - // Loop over all controls except the current one, and also - // skip any RESERVOIR_RATE controls, since we cannot - // handle those. - const int nwc = well_controls_get_num(wc); - int ctrl_index = 0; - for (; ctrl_index < nwc; ++ctrl_index) { - if (ctrl_index == current) { - // This is the currently used control, so it is - // used as an equation. So this is not used as an - // inequality constraint, and therefore skipped. - continue; - } - if (detail::constraintBroken( - xw.bhp(), xw.thp(), xw.wellRates(), - w, np, wells().type[w], wc, ctrl_index)) { - // ctrl_index will be the index of the broken constraint after the loop. - break; - } - } - if (ctrl_index != nwc) { - // Constraint number ctrl_index was broken, switch to it. - if (terminal_output_) - { - std::cout << "Switching control mode for well " << wells().name[w] - << " from " << modestring[well_controls_iget_type(wc, current)] - << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; - } - xw.currentControls()[w] = ctrl_index; - current = xw.currentControls()[w]; - } - - //Get gravity for THP hydrostatic corrrection - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - - // Updating well state and primary variables. - // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE - const double target = well_controls_iget_target(wc, current); - const double* distr = well_controls_iget_distr(wc, current); - switch (well_controls_iget_type(wc, current)) { - case BHP: - xw.bhp()[w] = target; - break; - - case THP: { - double aqua = 0.0; - double liquid = 0.0; - double vapour = 0.0; - - if (active_[ Water ]) { - aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; - } - if (active_[ Oil ]) { - liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; - } - if (active_[ Gas ]) { - vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; - } - - 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); - - //Set *BHP* target by calculating bhp from THP - const WellType& well_type = wells().type[w]; - - if (well_type == INJECTOR) { - double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(), - asImpl().stdWells().wellPerforationDensities(), gravity); - - xw.bhp()[w] = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; - } - else if (well_type == PRODUCER) { - double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(), - asImpl().stdWells().wellPerforationDensities(), gravity); - - xw.bhp()[w] = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; - } - else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - break; - } - - case RESERVOIR_RATE: - // No direct change to any observable quantity at - // surface condition. In this case, use existing - // flow rates as initial conditions as reservoir - // rate acts only in aggregate. - break; - - case SURFACE_RATE: - // 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) { - - // 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) { - numPhasesWithTargetsUnderThisControl += 1; - } - } - for (int phase = 0; phase < np; ++phase) { - if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { - xw.wellRates()[np*w + phase] = target * distr[phase]; - } - } - } else { - OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); - } - - - break; - } - - } - } - - - template @@ -1429,8 +1151,10 @@ namespace detail { ADB::V total_residual_v = total_residual.value(); const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); assert(dx.size() == total_residual_v.size()); - asImpl().updateWellState(dx.array(), well_state); - asImpl().updateWellControls(well_state); + // asImpl().updateWellState(dx.array(), well_state); + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + asImpl().stdWells().updateWellState(dx.array(), gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state); + asImpl().stdWells(). updateWellControls(fluid_.phaseUsage(), gravity, vfp_properties_, terminal_output_, active_, well_state); } } while (it < 15); @@ -1602,7 +1326,7 @@ namespace detail { //Perform hydrostatic correction to computed targets double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - const ADB::V dp_v = detail::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, asImpl().stdWells().wellPerforationDensities(), gravity); + const ADB::V dp_v = wellhelpers::computeHydrostaticCorrection(wells(), vfp_ref_depth_v, asImpl().stdWells().wellPerforationDensities(), gravity); const ADB dp = ADB::constant(dp_v); const ADB dp_inj = superset(subset(dp, thp_inj_elems), thp_inj_elems, nw); const ADB dp_prod = superset(subset(dp, thp_prod_elems), thp_prod_elems, nw); @@ -1960,7 +1684,10 @@ namespace detail { } - asImpl().updateWellState(dwells,well_state); + // TODO: gravity should be stored as a member + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + // asImpl().updateWellState(dwells,well_state); + asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state); // Update phase conditions used for property calculations. updatePhaseCondFromPrimalVariable(); @@ -1970,106 +1697,6 @@ namespace detail { - template - void - BlackoilModelBase::updateWellState(const V& dwells, - WellState& well_state) - { - - if( localWellsActive() ) - { - const int np = wells().number_of_phases; - const int nw = wells().number_of_wells; - - // Extract parts of dwells corresponding to each part. - int varstart = 0; - const V dqs = subset(dwells, Span(np*nw, 1, varstart)); - varstart += dqs.size(); - const V dbhp = subset(dwells, Span(nw, 1, varstart)); - varstart += dbhp.size(); - assert(varstart == dwells.size()); - const double dpmaxrel = dpMaxRel(); - - - // Qs update. - // Since we need to update the wellrates, that are ordered by wells, - // from dqs which are ordered by phase, the simplest is to compute - // dwr, which is the data from dqs but ordered by wells. - const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); - const V dwr = Eigen::Map(wwr.data(), nw*np); - const V wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); - const V wr = wr_old - dwr; - std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); - - // Bhp update. - const V bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); - const V dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); - const V bhp = bhp_old - dbhp_limited; - std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); - - //Get gravity for THP hydrostatic correction - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - - // Thp update - const Opm::PhaseUsage& pu = fluid_.phaseUsage(); - //Loop over all wells -#pragma omp parallel for schedule(static) - for (int w=0; wgetTable(table_id)->getDatumDepth(), - asImpl().stdWells().wellPerforationDensities(), gravity); - - well_state.thp()[w] = vfp_properties_.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp); - } - else if (well_type == PRODUCER) { - double dp = detail::computeHydrostaticCorrection( - wells(), w, vfp_properties_.getProd()->getTable(table_id)->getDatumDepth(), - asImpl().stdWells().wellPerforationDensities(), gravity); - - well_state.thp()[w] = vfp_properties_.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); - } - else { - OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); - } - - //Assume only one THP control specified for each well - break; - } - } - } - } - } - - - - - template std::vector BlackoilModelBase::computeRelPerm(const SolutionState& state) const diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index 84d63fa8d..6892771b8 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -964,7 +964,7 @@ namespace Opm { // inequality constraint, and therefore skipped. continue; } - if (detail::constraintBroken( + if (wellhelpers::constraintBroken( xw.bhp(), xw.thp(), xw.wellRates(), w, np, wellsMultiSegment()[w]->wellType(), wc, ctrl_index)) { // ctrl_index will be the index of the broken constraint after the loop. diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp index 02b7dd3ac..c7881a3a4 100644 --- a/opm/autodiff/BlackoilSolventModel.hpp +++ b/opm/autodiff/BlackoilSolventModel.hpp @@ -149,7 +149,7 @@ namespace Opm { using Base::dsMax; using Base::drMaxRel; using Base::maxResidualAllowed; - using Base::updateWellControls; + // using Base::updateWellControls; using Base::computeWellConnectionPressures; using Base::addWellControlEq; // using Base::computePropertiesForWellConnectionPressures; diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 005cd2fe3..492eb72e9 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -121,6 +121,22 @@ namespace Opm { const SolutionState& state, WellState& xw) const; + template + void updateWellState(const Vector& dwells, + const double gravity, + const double dpmaxrel, + const Opm::PhaseUsage& pu, + const std::vector& active, + const VFPProperties& vfp_properties, + WellState& well_state); + + template + void updateWellControls(const Opm::PhaseUsage& pu, + const double gravity, + const VFPProperties& vfp_properties, + const bool terminal_output, + const std::vector& active, + WellState& xw) const; protected: diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 1794478a3..40d4b3a3e 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -22,6 +22,10 @@ #include #include +#include +#include +#include + @@ -498,4 +502,259 @@ namespace Opm xw.perfPress().assign(perfpressure.data(), perfpressure.data() + nperf); } + + + + + template + void + StandardWells:: + updateWellState(const Vector& dwells, + const double gravity, + const double dpmaxrel, + const Opm::PhaseUsage& pu, + const std::vector& active, + const VFPProperties& vfp_properties, + WellState& well_state) + { + if( localWellsActive() ) + { + // TODO: these parameter should be stored in the StandardWells class + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; + + // Extract parts of dwells corresponding to each part. + int varstart = 0; + const Vector dqs = subset(dwells, Span(np*nw, 1, varstart)); + varstart += dqs.size(); + const Vector dbhp = subset(dwells, Span(nw, 1, varstart)); + varstart += dbhp.size(); + assert(varstart == dwells.size()); + + + // Qs update. + // Since we need to update the wellrates, that are ordered by wells, + // from dqs which are ordered by phase, the simplest is to compute + // dwr, which is the data from dqs but ordered by wells. + const DataBlock wwr = Eigen::Map(dqs.data(), np, nw).transpose(); + const Vector dwr = Eigen::Map(wwr.data(), nw*np); + const Vector wr_old = Eigen::Map(&well_state.wellRates()[0], nw*np); + const Vector wr = wr_old - dwr; + std::copy(&wr[0], &wr[0] + wr.size(), well_state.wellRates().begin()); + + // Bhp update. + const Vector bhp_old = Eigen::Map(&well_state.bhp()[0], nw, 1); + const Vector dbhp_limited = sign(dbhp) * dbhp.abs().min(bhp_old.abs()*dpmaxrel); + const Vector bhp = bhp_old - dbhp_limited; + std::copy(&bhp[0], &bhp[0] + bhp.size(), well_state.bhp().begin()); + + //Loop over all wells +#pragma omp parallel for schedule(static) + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + const int nwc = well_controls_get_num(wc); + //Loop over all controls until we find a THP control + //that specifies what we need... + //Will only update THP for wells with THP control + for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) { + if (well_controls_iget_type(wc, ctrl_index) == THP) { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + if (active[ Water ]) { + aqua = wr[w*np + pu.phase_pos[ Water ] ]; + } + if (active[ Oil ]) { + liquid = wr[w*np + pu.phase_pos[ Oil ] ]; + } + if (active[ Gas ]) { + vapour = wr[w*np + pu.phase_pos[ Gas ] ]; + } + + double alq = well_controls_iget_alq(wc, ctrl_index); + int table_id = well_controls_iget_vfp(wc, ctrl_index); + + const WellType& well_type = wells().type[w]; + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getInj()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity); + + well_state.thp()[w] = vfp_properties.getInj()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp); + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getProd()->getTable(table_id)->getDatumDepth(), + wellPerforationDensities(), gravity); + + well_state.thp()[w] = vfp_properties.getProd()->thp(table_id, aqua, liquid, vapour, bhp[w] + dp, alq); + } + else { + OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well"); + } + + //Assume only one THP control specified for each well + break; + } + } + } + } + } + + + + + + template + void + StandardWells:: + updateWellControls(const Opm::PhaseUsage& pu, + const double gravity, + const VFPProperties& vfp_properties, + const bool terminal_output, + const std::vector& active, + WellState& xw) const + { + if( !localWellsActive() ) return ; + + std::string modestring[4] = { "BHP", "THP", "RESERVOIR_RATE", "SURFACE_RATE" }; + // Find, for each well, if any constraints are broken. If so, + // switch control to first broken constraint. + const int np = wells().number_of_phases; + const int nw = wells().number_of_wells; +#pragma omp parallel for schedule(dynamic) + for (int w = 0; w < nw; ++w) { + const WellControls* wc = wells().ctrls[w]; + // The current control in the well state overrides + // the current control set in the Wells struct, which + // is instead treated as a default. + int current = xw.currentControls()[w]; + // Loop over all controls except the current one, and also + // skip any RESERVOIR_RATE controls, since we cannot + // handle those. + const int nwc = well_controls_get_num(wc); + int ctrl_index = 0; + for (; ctrl_index < nwc; ++ctrl_index) { + if (ctrl_index == current) { + // This is the currently used control, so it is + // used as an equation. So this is not used as an + // inequality constraint, and therefore skipped. + continue; + } + if (wellhelpers::constraintBroken( + xw.bhp(), xw.thp(), xw.wellRates(), + w, np, wells().type[w], wc, ctrl_index)) { + // ctrl_index will be the index of the broken constraint after the loop. + break; + } + } + if (ctrl_index != nwc) { + // Constraint number ctrl_index was broken, switch to it. + if (terminal_output) + { + std::cout << "Switching control mode for well " << wells().name[w] + << " from " << modestring[well_controls_iget_type(wc, current)] + << " to " << modestring[well_controls_iget_type(wc, ctrl_index)] << std::endl; + } + xw.currentControls()[w] = ctrl_index; + current = xw.currentControls()[w]; + } + + // Updating well state and primary variables. + // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE + const double target = well_controls_iget_target(wc, current); + const double* distr = well_controls_iget_distr(wc, current); + switch (well_controls_iget_type(wc, current)) { + case BHP: + xw.bhp()[w] = target; + break; + + case THP: { + double aqua = 0.0; + double liquid = 0.0; + double vapour = 0.0; + + if (active[ Water ]) { + aqua = xw.wellRates()[w*np + pu.phase_pos[ Water ] ]; + } + if (active[ Oil ]) { + liquid = xw.wellRates()[w*np + pu.phase_pos[ Oil ] ]; + } + if (active[ Gas ]) { + vapour = xw.wellRates()[w*np + pu.phase_pos[ Gas ] ]; + } + + 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); + + //Set *BHP* target by calculating bhp from THP + const WellType& well_type = wells().type[w]; + + if (well_type == INJECTOR) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getInj()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity); + + xw.bhp()[w] = vfp_properties.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp; + } + else if (well_type == PRODUCER) { + double dp = wellhelpers::computeHydrostaticCorrection( + wells(), w, vfp_properties.getProd()->getTable(vfp)->getDatumDepth(), + wellPerforationDensities(), gravity); + + xw.bhp()[w] = vfp_properties.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp; + } + else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + break; + } + + case RESERVOIR_RATE: + // No direct change to any observable quantity at + // surface condition. In this case, use existing + // flow rates as initial conditions as reservoir + // rate acts only in aggregate. + break; + + case SURFACE_RATE: + // 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) { + + // 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) { + numPhasesWithTargetsUnderThisControl += 1; + } + } + for (int phase = 0; phase < np; ++phase) { + if (distr[phase] > 0.0 && numPhasesWithTargetsUnderThisControl < 2 ) { + xw.wellRates()[np*w + phase] = target * distr[phase]; + } + } + } else { + OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well"); + } + + + break; + } + } + + } + } diff --git a/opm/autodiff/WellHelpers.hpp b/opm/autodiff/WellHelpers.hpp new file mode 100644 index 000000000..a580f2c53 --- /dev/null +++ b/opm/autodiff/WellHelpers.hpp @@ -0,0 +1,169 @@ +/* + Copyright 2016 SINTEF ICT, Applied Mathematics. + Copyright 2016 Statoil ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + +#ifndef OPM_WELLHELPERS_HEADER_INCLUDED +#define OPM_WELLHELPERS_HEADER_INCLUDED + +#include +#include +// #include + +#include + +namespace Opm { + + + // --------- Types --------- + typedef AutoDiffBlock ADB; + typedef ADB::V Vector; + + + namespace wellhelpers + { + inline + double rateToCompare(const std::vector& well_phase_flow_rate, + const int well, + const int num_phases, + const double* distr) + { + double rate = 0.0; + for (int phase = 0; phase < num_phases; ++phase) { + // Important: well_phase_flow_rate is ordered with all phase rates for first + // well first, then all phase rates for second well etc. + rate += well_phase_flow_rate[well*num_phases + phase] * distr[phase]; + } + return rate; + } + + inline + bool constraintBroken(const std::vector& bhp, + const std::vector& thp, + const std::vector& well_phase_flow_rate, + const int well, + const int num_phases, + const WellType& well_type, + const WellControls* wc, + const int ctrl_index) + { + const WellControlType ctrl_type = well_controls_iget_type(wc, ctrl_index); + const double target = well_controls_iget_target(wc, ctrl_index); + const double* distr = well_controls_iget_distr(wc, ctrl_index); + + bool broken = false; + + switch (well_type) { + case INJECTOR: + { + switch (ctrl_type) { + case BHP: + broken = bhp[well] > target; + break; + + case THP: + broken = thp[well] > target; + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + broken = rateToCompare(well_phase_flow_rate, + well, num_phases, distr) > target; + break; + } + } + break; + + case PRODUCER: + { + switch (ctrl_type) { + case BHP: + broken = bhp[well] < target; + break; + + case THP: + broken = thp[well] < target; + break; + + case RESERVOIR_RATE: // Intentional fall-through + case SURFACE_RATE: + // Note that the rates compared below are negative, + // so breaking the constraints means: too high flow rate + // (as for injection). + broken = rateToCompare(well_phase_flow_rate, + well, num_phases, distr) < target; + break; + } + } + break; + + default: + OPM_THROW(std::logic_error, "Can only handle INJECTOR and PRODUCER wells."); + } + + return broken; + } + + /** + * Simple hydrostatic correction for VFP table + * @param wells - wells struct + * @param w Well number + * @param vfp_table VFP table + * @param well_perforation_densities Densities at well perforations + * @param gravity Gravitational constant (e.g., 9.81...) + */ + inline + double computeHydrostaticCorrection(const Wells& wells, const int w, double vfp_ref_depth, + const Vector& well_perforation_densities, const double gravity) { + if ( wells.well_connpos[w] == wells.well_connpos[w+1] ) + { + // This is a well with no perforations. + // If this is the last well we would subscript over the + // bounds below. + // we assume well_perforation_densities to be 0 + return 0; + } + const double well_ref_depth = wells.depth_ref[w]; + const double dh = vfp_ref_depth - well_ref_depth; + const int perf = wells.well_connpos[w]; + const double rho = well_perforation_densities[perf]; + const double dp = rho*gravity*dh; + + return dp; + } + + inline + Vector computeHydrostaticCorrection(const Wells& wells, const Vector vfp_ref_depth, + const Vector& well_perforation_densities, const double gravity) { + const int nw = wells.number_of_wells; + Vector retval = Vector::Zero(nw); + +#pragma omp parallel for schedule(static) + for (int i=0; i Date: Mon, 11 Apr 2016 17:10:15 +0200 Subject: [PATCH 10/13] adding extractWellPerfProperties to StandardWellsSolvent to fix the flow_solvent running. --- opm/autodiff/BlackoilModelBase.hpp | 2 + opm/autodiff/BlackoilModelBase_impl.hpp | 2 +- opm/autodiff/BlackoilSolventModel_impl.hpp | 2 +- opm/autodiff/StandardWells.hpp | 7 ++- opm/autodiff/StandardWellsSolvent.hpp | 13 +++++- opm/autodiff/StandardWellsSolvent_impl.hpp | 52 +++++++++++++++++++++- opm/autodiff/StandardWells_impl.hpp | 7 ++- 7 files changed, 77 insertions(+), 8 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index ccdbb45f2..2c432feb9 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -386,6 +386,8 @@ namespace Opm { void assembleMassBalanceEq(const SolutionState& state); + // TODO: only kept for now due to flow_multisegment + // will be removed soon void extractWellPerfProperties(const SolutionState& state, std::vector& mob_perfcells, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index f19e33077..9f158d243 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -855,7 +855,7 @@ namespace detail { std::vector mob_perfcells; std::vector b_perfcells; - asImpl().stdWells().extractWellPerfProperties(rq_, fluid_.numPhases(), mob_perfcells, b_perfcells); + asImpl().stdWells().extractWellPerfProperties(state, rq_, fluid_.numPhases(), fluid_, active_, mob_perfcells, b_perfcells); if (param_.solve_welleq_initially_ && initial_assembly) { // solve the well equations as a pre-processing step asImpl().solveWellEq(mob_perfcells, b_perfcells, state, well_state); diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 6145267df..bb86bc578 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -90,7 +90,7 @@ namespace Opm { solvent_pos_(detail::solventPos(fluid.phaseUsage())), solvent_props_(solvent_props), is_miscible_(is_miscible), - std_wells_(wells_arg, solvent_props) + std_wells_(wells_arg, solvent_props, solvent_pos_) { if (has_solvent_) { diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 492eb72e9..6adba9d16 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -101,9 +101,12 @@ namespace Opm { const std::vector& depth_perf, const double grav); - template - void extractWellPerfProperties(const std::vector& rq, + template + void extractWellPerfProperties(const SolutionState& state, + const std::vector& rq, const int np, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, std::vector& mob_perfcells, std::vector& b_perfcells) const; diff --git a/opm/autodiff/StandardWellsSolvent.hpp b/opm/autodiff/StandardWellsSolvent.hpp index 9906eae3b..083b30ab7 100644 --- a/opm/autodiff/StandardWellsSolvent.hpp +++ b/opm/autodiff/StandardWellsSolvent.hpp @@ -36,7 +36,7 @@ namespace Opm { using Base = StandardWells; // --------- Public methods --------- - explicit StandardWellsSolvent(const Wells* wells, const SolventPropsAdFromDeck& solvent_props); + explicit StandardWellsSolvent(const Wells* wells, const SolventPropsAdFromDeck& solvent_props, const int solvent_pos); template void computePropertiesForWellConnectionPressures(const SolutionState& state, @@ -48,8 +48,19 @@ namespace Opm { std::vector& rsmax_perf, std::vector& rvmax_perf, std::vector& surf_dens_perf); + + // TODO: fluid and active may be can put in the member list + template + void extractWellPerfProperties(const SolutionState& state, + const std::vector& rq, + const int np, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + std::vector& mob_perfcells, + std::vector& b_perfcells) const; protected: const SolventPropsAdFromDeck& solvent_props_; + const int solvent_pos_; }; diff --git a/opm/autodiff/StandardWellsSolvent_impl.hpp b/opm/autodiff/StandardWellsSolvent_impl.hpp index 38b538675..ee69a342d 100644 --- a/opm/autodiff/StandardWellsSolvent_impl.hpp +++ b/opm/autodiff/StandardWellsSolvent_impl.hpp @@ -33,9 +33,11 @@ namespace Opm StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg, - const SolventPropsAdFromDeck& solvent_props) + const SolventPropsAdFromDeck& solvent_props, + const int solvent_pos) : Base(wells_arg) , solvent_props_(solvent_props) + , solvent_pos_(solvent_pos) { } @@ -178,4 +180,52 @@ namespace Opm } + + + + + template + void + StandardWellsSolvent:: + extractWellPerfProperties(const SolutionState& state, + const std::vector& rq, + const int np, + const BlackoilPropsAdInterface& fluid, + const std::vector& active, + std::vector& mob_perfcells, + std::vector& b_perfcells) const + { + Base::extractWellPerfProperties(state, rq, np, fluid, active, mob_perfcells, b_perfcells); + // handle the solvent related + { + int gas_pos = fluid.phaseUsage().phase_pos[Gas]; + const std::vector& well_cells = wellOps().well_cells; + const int nperf = well_cells.size(); + // Gas and solvent is combinded and solved together + // The input in the well equation is then the + // total gas phase = hydro carbon gas + solvent gas + + // The total mobility is the sum of the solvent and gas mobiliy + mob_perfcells[gas_pos] += subset(rq[solvent_pos_].mob, well_cells); + + // A weighted sum of the b-factors of gas and solvent are used. + const int nc = rq[solvent_pos_].mob.size(); + + const Opm::PhaseUsage& pu = fluid.phaseUsage(); + const ADB zero = ADB::constant(Vector::Zero(nc)); + const ADB& ss = state.solvent_saturation; + const ADB& sg = (active[ Gas ] + ? state.saturation[ pu.phase_pos[ Gas ] ] + : zero); + + Selector zero_selector(ss.value() + sg.value(), Selector::Zero); + ADB F_solvent = subset(zero_selector.select(ss, ss / (ss + sg)),well_cells); + Vector ones = Vector::Constant(nperf,1.0); + + b_perfcells[gas_pos] = (ones - F_solvent) * b_perfcells[gas_pos]; + b_perfcells[gas_pos] += (F_solvent * subset(rq[solvent_pos_].b, well_cells)); + } + } + + } diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 40d4b3a3e..2466b0166 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -289,11 +289,14 @@ namespace Opm - template + template void StandardWells:: - extractWellPerfProperties(const std::vector& rq, + extractWellPerfProperties(const SolutionState& /* state */, + const std::vector& rq, const int np, + const BlackoilPropsAdInterface& /* fluid */, + const std::vector& /* active */, std::vector& mob_perfcells, std::vector& b_perfcells) const { From 14d774e08f2429f6b6c0bba952ba4ff6c0b924fb Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 12 Apr 2016 11:02:04 +0200 Subject: [PATCH 11/13] removing blank lines in CMakeLists_files Not sure how the lines happened, probably from the rebase process. --- CMakeLists_files.cmake | 4 ---- 1 file changed, 4 deletions(-) diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index bd7868398..692905f2a 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -26,7 +26,6 @@ # originally generated with the command: # find opm -name '*.c*' -printf '\t%p\n' | sort list (APPEND MAIN_SOURCE_FILES - opm/autodiff/BlackoilPropsAdInterface.cpp opm/autodiff/ExtractParallelGridInformationToISTL.cpp opm/autodiff/NewtonIterationBlackoilCPR.cpp @@ -68,7 +67,6 @@ list (APPEND MAIN_SOURCE_FILES ) - # originally generated with the command: # find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort list (APPEND TEST_SOURCE_FILES @@ -131,7 +129,6 @@ list (APPEND PROGRAM_SOURCE_FILES # originally generated with the command: # find opm -name '*.h*' -a ! -name '*-pch.hpp' -printf '\t%p\n' | sort list (APPEND PUBLIC_HEADER_FILES - opm/autodiff/AdditionalObjectDeleter.hpp opm/autodiff/AutoDiffBlock.hpp opm/autodiff/AutoDiffHelpers.hpp @@ -228,4 +225,3 @@ list (APPEND PUBLIC_HEADER_FILES opm/simulators/SimulatorIncompTwophase.hpp ) - From 75b73a893e322f0ab371c9d26f380aeb6cba01df Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 12 Apr 2016 11:35:03 +0200 Subject: [PATCH 12/13] putting Vector and ADB to be inside StandardWells and wellhelpers --- opm/autodiff/BlackoilModelBase_impl.hpp | 4 ++-- opm/autodiff/StandardWells.hpp | 24 ++++++++++++------------ opm/autodiff/StandardWells_impl.hpp | 8 ++++---- opm/autodiff/WellHelpers.hpp | 7 ++++--- 4 files changed, 22 insertions(+), 21 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 9f158d243..c4e77413d 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -780,8 +780,8 @@ namespace detail { asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf); // Extract well connection depths. - const Vector depth = cellCentroidsZToEigen(grid_); - const Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells); + const StandardWells::Vector depth = cellCentroidsZToEigen(grid_); + const StandardWells::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells); const int nperf = wells().well_connpos[wells().number_of_wells]; const std::vector depth_perf(pdepth.data(), pdepth.data() + nperf); diff --git a/opm/autodiff/StandardWells.hpp b/opm/autodiff/StandardWells.hpp index 6adba9d16..33a15812f 100644 --- a/opm/autodiff/StandardWells.hpp +++ b/opm/autodiff/StandardWells.hpp @@ -37,16 +37,6 @@ namespace Opm { - // --------- Types --------- - typedef AutoDiffBlock ADB; - typedef ADB::V Vector; - - // copied from BlackoilModelBase - // should put to somewhere better - typedef Eigen::Array DataBlock; /// Class for handling the standard well model. class StandardWells { @@ -59,6 +49,16 @@ namespace Opm { }; public: + // --------- Types --------- + using ADB = AutoDiffBlock; + using Vector = ADB::V; + + // copied from BlackoilModelBase + // should put to somewhere better + using DataBlock = Eigen::Array; // --------- Public methods --------- explicit StandardWells(const Wells* wells); @@ -73,11 +73,11 @@ namespace Opm { const WellOps& wellOps() const; /// Density of each well perforation - Vector& wellPerforationDensities(); + Vector& wellPerforationDensities(); // mutable version kept for BlackoilMultisegmentModel const Vector& wellPerforationDensities() const; /// Diff to bhp for each well perforation. - Vector& wellPerforationPressureDiffs(); + Vector& wellPerforationPressureDiffs(); // mutable version kept for BlackoilMultisegmentModel const Vector& wellPerforationPressureDiffs() const; template diff --git a/opm/autodiff/StandardWells_impl.hpp b/opm/autodiff/StandardWells_impl.hpp index 2466b0166..b086984f5 100644 --- a/opm/autodiff/StandardWells_impl.hpp +++ b/opm/autodiff/StandardWells_impl.hpp @@ -130,7 +130,7 @@ namespace Opm - Vector& StandardWells::wellPerforationDensities() + StandardWells::Vector& StandardWells::wellPerforationDensities() { return well_perforation_densities_; } @@ -139,7 +139,7 @@ namespace Opm - const Vector& + const StandardWells::Vector& StandardWells::wellPerforationDensities() const { return well_perforation_densities_; @@ -149,7 +149,7 @@ namespace Opm - Vector& + StandardWells::Vector& StandardWells::wellPerforationPressureDiffs() { return well_perforation_pressure_diffs_; @@ -159,7 +159,7 @@ namespace Opm - const Vector& + const StandardWells::Vector& StandardWells::wellPerforationPressureDiffs() const { return well_perforation_pressure_diffs_; diff --git a/opm/autodiff/WellHelpers.hpp b/opm/autodiff/WellHelpers.hpp index a580f2c53..b92bb0682 100644 --- a/opm/autodiff/WellHelpers.hpp +++ b/opm/autodiff/WellHelpers.hpp @@ -31,9 +31,6 @@ namespace Opm { - // --------- Types --------- - typedef AutoDiffBlock ADB; - typedef ADB::V Vector; namespace wellhelpers @@ -120,6 +117,10 @@ namespace Opm { return broken; } + + // --------- Types --------- + using Vector = AutoDiffBlock::V; + /** * Simple hydrostatic correction for VFP table * @param wells - wells struct From adb140a78836649265cff9944805470e31c2ec0d Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Tue, 12 Apr 2016 14:00:22 +0200 Subject: [PATCH 13/13] fixing the running of flow_multisegment --- opm/autodiff/BlackoilModelBase.hpp | 6 + opm/autodiff/BlackoilModelBase_impl.hpp | 24 +++- opm/autodiff/BlackoilMultiSegmentModel.hpp | 8 ++ .../BlackoilMultiSegmentModel_impl.hpp | 108 +++++++++++++++++- 4 files changed, 142 insertions(+), 4 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp index 2c432feb9..41af2dc3c 100644 --- a/opm/autodiff/BlackoilModelBase.hpp +++ b/opm/autodiff/BlackoilModelBase.hpp @@ -393,6 +393,12 @@ namespace Opm { std::vector& mob_perfcells, std::vector& b_perfcells) const; + // TODO: only kept for now due to flow_multisegment + // will be removed soon + void updateWellState(const V& dwells, + WellState& well_state); + + bool solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index c4e77413d..7dc6dc2b8 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1685,9 +1685,9 @@ namespace detail { // TODO: gravity should be stored as a member - const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); - // asImpl().updateWellState(dwells,well_state); - asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state); + // const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + // asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), active_, vfp_properties_, well_state); + asImpl().updateWellState(dwells,well_state); // Update phase conditions used for property calculations. updatePhaseCondFromPrimalVariable(); @@ -2527,6 +2527,24 @@ namespace detail { + + // TODO: only kept for now due to flow_multisegment + // will be removed soon + template + void + BlackoilModelBase::updateWellState(const V& dwells, + WellState& well_state) + { + const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_)); + asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(), + active_, vfp_properties_, well_state); + + } + + + + + } // namespace Opm #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp index b59e709c2..1e2c68c31 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp @@ -222,6 +222,7 @@ namespace Opm { using Base::convergenceReduction; using Base::maxResidualAllowed; using Base::variableState; + using Base::variableWellStateIndices; using Base::asImpl; const std::vector& wellsMultiSegment() const { return wells_multisegment_; } @@ -239,6 +240,13 @@ namespace Opm { void computeWellConnectionPressures(const SolutionState& state, const WellState& xw); + /// added to fixing the flow_multisegment running + bool + baseSolveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state); + bool solveWellEq(const std::vector& mob_perfcells, const std::vector& b_perfcells, diff --git a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp index 6892771b8..c87c2c7ac 100644 --- a/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp +++ b/opm/autodiff/BlackoilMultiSegmentModel_impl.hpp @@ -1033,7 +1033,7 @@ namespace Opm { SolutionState& state, WellState& well_state) { - const bool converged = Base::solveWellEq(mob_perfcells, b_perfcells, state, well_state); + const bool converged = baseSolveWellEq(mob_perfcells, b_perfcells, state, well_state); if (converged) { // We must now update the state.segp and state.segqs members, @@ -1547,6 +1547,112 @@ namespace Opm { } + + + /// added to fixing the flow_multisegment running + template + bool + BlackoilMultiSegmentModel::baseSolveWellEq(const std::vector& mob_perfcells, + const std::vector& b_perfcells, + SolutionState& state, + WellState& well_state) { + V aliveWells; + const int np = wells().number_of_phases; + std::vector cq_s(np, ADB::null()); + std::vector indices = variableWellStateIndices(); + SolutionState state0 = state; + WellState well_state0 = well_state; + makeConstantState(state0); + + std::vector mob_perfcells_const(np, ADB::null()); + std::vector b_perfcells_const(np, ADB::null()); + + if ( Base::localWellsActive() ){ + // If there are non well in the sudomain of the process + // thene mob_perfcells_const and b_perfcells_const would be empty + for (int phase = 0; phase < np; ++phase) { + mob_perfcells_const[phase] = ADB::constant(mob_perfcells[phase].value()); + b_perfcells_const[phase] = ADB::constant(b_perfcells[phase].value()); + } + } + + int it = 0; + bool converged; + do { + // bhp and Q for the wells + std::vector vars0; + vars0.reserve(2); + variableWellStateInitials(well_state, vars0); + std::vector vars = ADB::variables(vars0); + + SolutionState wellSolutionState = state0; + variableStateExtractWellsVars(indices, vars, wellSolutionState); + computeWellFlux(wellSolutionState, mob_perfcells_const, b_perfcells_const, aliveWells, cq_s); + updatePerfPhaseRatesAndPressures(cq_s, wellSolutionState, well_state); + addWellFluxEq(cq_s, wellSolutionState); + addWellControlEq(wellSolutionState, well_state, aliveWells); + converged = Base::getWellConvergence(it); + + if (converged) { + break; + } + + ++it; + if( Base::localWellsActive() ) + { + std::vector eqs; + eqs.reserve(2); + eqs.push_back(residual_.well_flux_eq); + eqs.push_back(residual_.well_eq); + ADB total_residual = vertcatCollapseJacs(eqs); + const std::vector& Jn = total_residual.derivative(); + typedef Eigen::SparseMatrix Sp; + Sp Jn0; + Jn[0].toSparse(Jn0); + const Eigen::SparseLU< Sp > solver(Jn0); + ADB::V total_residual_v = total_residual.value(); + const Eigen::VectorXd& dx = solver.solve(total_residual_v.matrix()); + assert(dx.size() == total_residual_v.size()); + // asImpl().updateWellState(dx.array(), well_state); + updateWellState(dx.array(), well_state); + updateWellControls(well_state); + } + } while (it < 15); + + if (converged) { + if ( terminal_output_ ) { + std::cout << "well converged iter: " << it << std::endl; + } + const int nw = wells().number_of_wells; + { + // We will set the bhp primary variable to the new ones, + // but we do not change the derivatives here. + ADB::V new_bhp = Eigen::Map(well_state.bhp().data(), nw); + // Avoiding the copy below would require a value setter method + // in AutoDiffBlock. + std::vector old_derivs = state.bhp.derivative(); + state.bhp = ADB::function(std::move(new_bhp), std::move(old_derivs)); + } + { + // Need to reshuffle well rates, from phase running fastest + // to wells running fastest. + // The transpose() below switches the ordering. + const DataBlock wrates = Eigen::Map(well_state.wellRates().data(), nw, np).transpose(); + ADB::V new_qs = Eigen::Map(wrates.data(), nw*np); + std::vector old_derivs = state.qs.derivative(); + state.qs = ADB::function(std::move(new_qs), std::move(old_derivs)); + } + computeWellConnectionPressures(state, well_state); + } + + if (!converged) { + well_state = well_state0; + } + + return converged; + } + + } // namespace Opm #endif // OPM_BLACKOILMODELBASE_IMPL_HEADER_INCLUDED