From 738ec64ac80140d4dfa2dd794d63138a40c7f25c Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 09:48:49 +0200 Subject: [PATCH 01/22] Added necessary input parameters for well checking in spu_2p --- examples/spu_2p.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 2c819b32..bcab7c78 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -297,6 +297,8 @@ main(int argc, char** argv) boost::scoped_ptr rock_comp; Opm::SimulatorTimer simtimer; Opm::TwophaseState state; + bool check_well_controls = false; + int max_well_control_iterations = 0; double gravity[3] = { 0.0 }; if (use_deck) { std::string deck_filename = param.get("deck_filename"); @@ -309,6 +311,8 @@ main(int argc, char** argv) props.reset(new Opm::IncompPropertiesFromDeck(deck, global_cell)); // Wells init. wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); + check_well_controls = param.getDefault("check_well_controls", false); + max_well_control_iterations = param.getDefault("max_well_control_iterations", 0); // Timer init. if (deck.hasField("TSTEP")) { simtimer.init(deck); From 338f5907a640821878fb619c71f1ae53b373d1b8 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 10:14:38 +0200 Subject: [PATCH 02/22] Added well controls in spu2p --- examples/spu_2p.cpp | 81 +++++++++++++++++++++++++++++---------------- 1 file changed, 53 insertions(+), 28 deletions(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index bcab7c78..41115469 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -516,6 +516,8 @@ main(int argc, char** argv) Opm::WellReport wellreport; std::vector well_bhp; std::vector well_perfrates; + std::vector fractional_flows; + std::vector well_resflows; if (wells->c_wells()) { const int nw = wells->c_wells()->number_of_wells; well_bhp.resize(nw, 0.0); @@ -539,38 +541,61 @@ main(int argc, char** argv) if (wells->c_wells()) { Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), gravity[2], true, wdp); } - pressure_timer.start(); - if (rock_comp->isActive()) { - rc.resize(num_cells); - std::vector initial_pressure = state.pressure(); - std::vector prev_pressure; - for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { - prev_pressure = state.pressure(); - for (int cell = 0; cell < num_cells; ++cell) { - rc[cell] = rock_comp->rockComp(state.pressure()[cell]); + if (check_well_controls) { + computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows); + } + bool well_control_passed = !check_well_controls; + int well_control_iteration = 0; + do { + pressure_timer.start(); + if (rock_comp->isActive()) { + rc.resize(num_cells); + std::vector initial_pressure = state.pressure(); + std::vector prev_pressure; + for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { + prev_pressure = state.pressure(); + for (int cell = 0; cell < num_cells; ++cell) { + rc[cell] = rock_comp->rockComp(state.pressure()[cell]); + } + state.pressure() = initial_pressure; + psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(), + state.pressure(), state.faceflux(), well_bhp, well_perfrates); + double max_change = 0.0; + for (int cell = 0; cell < num_cells; ++cell) { + max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell])); + } + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + if (max_change < nl_pressure_tolerance) { + break; + } } - state.pressure() = initial_pressure; - psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(), - state.pressure(), state.faceflux(), well_bhp, well_perfrates); - double max_change = 0.0; - for (int cell = 0; cell < num_cells; ++cell) { - max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell])); + computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + } else { + psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), + well_bhp, well_perfrates); + } + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + + if (check_well_controls) { + Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), + fractional_flows, + well_perfrates, + well_resflows); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells->conditionsMet(well_bhp, well_resflows, well_resflows); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations) { + THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); } - std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; - if (max_change < nl_pressure_tolerance) { - break; + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; } } - computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); - } else { - psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), - well_bhp, well_perfrates); - } - pressure_timer.stop(); - double pt = pressure_timer.secsSinceStart(); - std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; - ptime += pt; - + } while (!well_control_passed); // Process transport sources (to include bdy terms and well flows). Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, wells->c_wells(), well_perfrates, reorder_src); From 8e49914bbde4c5bec181a2971cc1183ca557cd9a Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 10:17:58 +0200 Subject: [PATCH 03/22] Added friendly printout --- examples/spu_2p.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 41115469..15ece9a3 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -312,7 +312,7 @@ main(int argc, char** argv) // Wells init. wells.reset(new Opm::WellsManager(deck, *grid->c_grid(), props->permeability())); check_well_controls = param.getDefault("check_well_controls", false); - max_well_control_iterations = param.getDefault("max_well_control_iterations", 0); + max_well_control_iterations = param.getDefault("max_well_control_iterations", 10); // Timer init. if (deck.hasField("TSTEP")) { simtimer.init(deck); @@ -593,6 +593,8 @@ main(int argc, char** argv) } if (!well_control_passed) { std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; } } } while (!well_control_passed); From 5cf205f225839fc0dd2494c53154890c80bc9be4 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 11:28:56 +0200 Subject: [PATCH 04/22] Added failure checking for solve method. --- opm/core/pressure/IncompTpfa.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/opm/core/pressure/IncompTpfa.cpp b/opm/core/pressure/IncompTpfa.cpp index cd0a07b6..95973d3e 100644 --- a/opm/core/pressure/IncompTpfa.cpp +++ b/opm/core/pressure/IncompTpfa.cpp @@ -210,13 +210,16 @@ namespace Opm F.bc = bcs; F.totmob = &totmob[0]; F.wdp = &wdp[0]; - + int ok = true; if (rock_comp.empty()) { - ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_); + ok = ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_); } else { - ifs_tpfa_assemble_comprock(gg, &F, &trans_[0], &gpress_omegaweighted_[0], + ok = ifs_tpfa_assemble_comprock(gg, &F, &trans_[0], &gpress_omegaweighted_[0], &porevol[0], &rock_comp[0], dt, &pressure[0], h_); } + if (!ok) { + THROW("Failed assembling pressure system."); + } linsolver_.solve(h_->A, h_->b, h_->x); From c72c4eed7c94c9007c52fe62f794ecdf99d85133 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 11:29:13 +0200 Subject: [PATCH 05/22] Corrected some basic logic tests. --- opm/core/WellsGroup.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 931087a5..7b56e8ca 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -453,7 +453,7 @@ namespace Opm // as that would check if we're under group control, something we're not. const double children_guide_rate = children_[i]->productionGuideRate(true); children_[i]->applyProdGroupControl(prod_mode, - (children_guide_rate / my_guide_rate) * getTarget(prod_mode), + (children_guide_rate / my_guide_rate) * getTarget(prod_mode), false); } break; @@ -691,7 +691,7 @@ namespace Opm { // Not changing if we're not forced to change if (!forced - && (injSpec().control_mode_ != InjectionSpecification::GRUP || injSpec().control_mode_ != InjectionSpecification::NONE)) { + && (injSpec().control_mode_ != InjectionSpecification::GRUP && injSpec().control_mode_ != InjectionSpecification::NONE)) { return; } if (!wells_->type[self_index_] == INJECTOR) { @@ -734,7 +734,8 @@ namespace Opm { // Not changing if we're not forced to change if (!forced && (prodSpec().control_mode_ != ProductionSpecification::GRUP - || prodSpec().control_mode_ != ProductionSpecification::NONE)) { + && prodSpec().control_mode_ != ProductionSpecification::NONE)) { + std::cout << "Returning" << std::endl; return; } if (!wells_->type[self_index_] == PRODUCER) { @@ -771,6 +772,7 @@ namespace Opm distr[phase_pos[BlackoilPhases::Vapour]] = 1.0; break; case ProductionSpecification::LRAT: + std::cout << "applying rate" << std::endl; wct = SURFACE_RATE; if (!phase_used[BlackoilPhases::Liquid]) { THROW("Oil phase not active and LRAT control specified."); From 965b6d54d8f499634af606c01b3f789db9077a8d Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 12:35:26 +0200 Subject: [PATCH 06/22] Added summation method to find total produced rates in a given group. --- opm/core/WellsGroup.cpp | 33 +++++++++++++++++++++++++++++++++ opm/core/WellsGroup.hpp | 28 ++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 7b56e8ca..b145793d 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -528,7 +528,21 @@ namespace Opm return sum; } + /// Gets the total production flow of the given phase. + /// \param[in] phase_flows A vector containing rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] phase The phase for which to sum up. + double WellsGroup::getTotalProductionFlow(const std::vector& phase_flows, + const BlackoilPhases::PhaseIndex phase) + { + double sum = 0.0; + for (size_t i = 0; i < children_.size(); ++i) { + sum += children_[i]->getTotalProductionFlow(phase_flows, phase); + } + return sum; + } // ============== WellNode members ============ @@ -728,6 +742,25 @@ namespace Opm } + /// Gets the total production flow of the given phase. + /// \param[in] phase_flows A vector containing rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] phase The phase for which to sum up. + + double WellNode::getTotalProductionFlow(const std::vector& phase_flows, + const BlackoilPhases::PhaseIndex phase) + { + if (type() == INJECTOR) { + return 0.0; + } + return phase_flows[self_index_*phaseUsage().num_phases + phaseUsage().phase_pos[phase]]; + } + + WellType WellNode::type() const { + return wells_->type[self_index_]; + } + void WellNode::applyProdGroupControl(const ProductionSpecification::ControlMode control_mode, const double target, const bool forced) diff --git a/opm/core/WellsGroup.hpp b/opm/core/WellsGroup.hpp index 6fe5508b..d1b98aed 100644 --- a/opm/core/WellsGroup.hpp +++ b/opm/core/WellsGroup.hpp @@ -179,6 +179,15 @@ namespace Opm /// wells under group control virtual double injectionGuideRate(bool only_group) = 0; + /// Gets the total production flow of the given phase. + /// \param[in] phase_flows A vector containing rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] phase The phase for which to sum up. + virtual double getTotalProductionFlow(const std::vector& phase_flows, + const BlackoilPhases::PhaseIndex phase) = 0; + + protected: /// Calculates the correct rate for the given ProductionSpecification::ControlMode double rateByMode(const double* res_rates, @@ -258,6 +267,14 @@ namespace Opm /// \param[in] only_group If true, will only accumelate guide rates for /// wells under group control virtual double injectionGuideRate(bool only_group); + + /// Gets the total production flow of the given phase. + /// \param[in] phase_flows A vector containing rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] phase The phase for which to sum up. + virtual double getTotalProductionFlow(const std::vector& phase_flows, + const BlackoilPhases::PhaseIndex phase); private: std::vector > children_; @@ -327,6 +344,17 @@ namespace Opm /// \param[in] only_group If true, will only accumelate guide rates for /// wells under group control virtual double injectionGuideRate(bool only_group); + + /// Gets the total production flow of the given phase. + /// \param[in] phase_flows A vector containing rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] phase The phase for which to sum up. + virtual double getTotalProductionFlow(const std::vector& phase_flows, + const BlackoilPhases::PhaseIndex phase); + + /// Returns the type of the well. + WellType type() const; private: Wells* wells_; From 2394a474fcf3a63049b24809dce713f05e130391 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 12:57:36 +0200 Subject: [PATCH 07/22] Added simple reinjection to group control (explicit) --- opm/core/WellsGroup.cpp | 53 +++++++++++++++++++++++++++++++++++++++++ opm/core/WellsGroup.hpp | 36 ++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index b145793d..a23fbcbe 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -544,6 +544,45 @@ namespace Opm return sum; } + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + void WellsGroup::applyExplicitReinjectionControls(const std::vector&, + const std::vector& well_surfacerates_phase) + { + if (injSpec().control_mode_ == InjectionSpecification::REIN) { + BlackoilPhases::PhaseIndex phase; + switch (injSpec().injector_type_) { + case InjectionSpecification::WATER: + phase = BlackoilPhases::Aqua; + break; + case InjectionSpecification::GAS: + phase = BlackoilPhases::Vapour; + break; + case InjectionSpecification::OIL: + phase = BlackoilPhases::Liquid; + break; + } + const double total_produced = getTotalProductionFlow(well_surfacerates_phase, phase); + const double my_guide_rate = injectionGuideRate(true); + for (size_t i = 0; i < children_.size(); ++i) { + // Apply for all children. + // Note, we do _not_ want to call the applyProdGroupControl in this object, + // as that would check if we're under group control, something we're not. + const double children_guide_rate = children_[i]->injectionGuideRate(true); + children_[i]->applyInjGroupControl(InjectionSpecification::RATE, + (children_guide_rate / my_guide_rate) * total_produced, + false); + } + } + } + // ============== WellNode members ============ @@ -761,6 +800,20 @@ namespace Opm return wells_->type[self_index_]; } + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + void WellNode::applyExplicitReinjectionControls(const std::vector&, + const std::vector&) + { + // Do nothing at well level. + } void WellNode::applyProdGroupControl(const ProductionSpecification::ControlMode control_mode, const double target, const bool forced) diff --git a/opm/core/WellsGroup.hpp b/opm/core/WellsGroup.hpp index d1b98aed..3c14950d 100644 --- a/opm/core/WellsGroup.hpp +++ b/opm/core/WellsGroup.hpp @@ -187,6 +187,18 @@ namespace Opm virtual double getTotalProductionFlow(const std::vector& phase_flows, const BlackoilPhases::PhaseIndex phase) = 0; + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + virtual void applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase) = 0; + protected: /// Calculates the correct rate for the given ProductionSpecification::ControlMode @@ -275,6 +287,18 @@ namespace Opm /// \param[in] phase The phase for which to sum up. virtual double getTotalProductionFlow(const std::vector& phase_flows, const BlackoilPhases::PhaseIndex phase); + + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + virtual void applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase); private: std::vector > children_; @@ -355,6 +379,18 @@ namespace Opm /// Returns the type of the well. WellType type() const; + + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + virtual void applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase); private: Wells* wells_; From 50cb06a592641676c98f8a92ee229c1ac429882b Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 12:58:29 +0200 Subject: [PATCH 08/22] Added helpful comment. --- opm/core/WellsGroup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index a23fbcbe..13caadaa 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -553,7 +553,7 @@ namespace Opm /// A vector containing surface rates by phase for each well. /// Is assumed to be ordered the same way as the related Wells-struct, /// with all phase rates of a single well adjacent in the array. - void WellsGroup::applyExplicitReinjectionControls(const std::vector&, + void WellsGroup::applyExplicitReinjectionControls(const std::vector& /*well_reservoirrates_phase*/, const std::vector& well_surfacerates_phase) { if (injSpec().control_mode_ == InjectionSpecification::REIN) { From 9c7097381d8fd9463bb942f8d9ce016cf679c61c Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 13:03:37 +0200 Subject: [PATCH 09/22] Made relevant methods in WellCollection and WellsManager to call the reinjection control functions --- opm/core/WellCollection.cpp | 18 ++++++++++++++++++ opm/core/WellCollection.hpp | 12 ++++++++++++ opm/core/WellsManager.cpp | 16 ++++++++++++++++ opm/core/WellsManager.hpp | 12 ++++++++++++ 4 files changed, 58 insertions(+) diff --git a/opm/core/WellCollection.cpp b/opm/core/WellCollection.cpp index f873c259..585c5eaa 100644 --- a/opm/core/WellCollection.cpp +++ b/opm/core/WellCollection.cpp @@ -122,4 +122,22 @@ namespace Opm roots_[i]->applyInjGroupControls(); } } + + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + + void WellCollection::applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase) + { + for (size_t i = 0; i < roots_.size(); ++i) { + roots_[i]->applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase); + } + } } diff --git a/opm/core/WellCollection.hpp b/opm/core/WellCollection.hpp index a5db76d7..3e71ab02 100644 --- a/opm/core/WellCollection.hpp +++ b/opm/core/WellCollection.hpp @@ -88,6 +88,18 @@ namespace Opm /// Applies all group controls (injection and production) void applyGroupControls(); + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + void applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase); + private: // To account for the possibility of a forest std::vector > roots_; diff --git a/opm/core/WellsManager.cpp b/opm/core/WellsManager.cpp index 6d3b782e..d0d79b9b 100644 --- a/opm/core/WellsManager.cpp +++ b/opm/core/WellsManager.cpp @@ -696,4 +696,20 @@ namespace Opm well_surfacerates_phase); } + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + + void WellsManager::applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase) + { + well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase); + } + } // namespace Opm diff --git a/opm/core/WellsManager.hpp b/opm/core/WellsManager.hpp index cdae43c1..61fb32c3 100644 --- a/opm/core/WellsManager.hpp +++ b/opm/core/WellsManager.hpp @@ -86,6 +86,18 @@ namespace Opm bool conditionsMet(const std::vector& well_bhp, const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase); + + /// Applies explicit reinjection controls. This must be called at each timestep to be correct. + /// \param[in] well_reservoirrates_phase + /// A vector containing reservoir rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + /// \param[in] well_surfacerates_phase + /// A vector containing surface rates by phase for each well. + /// Is assumed to be ordered the same way as the related Wells-struct, + /// with all phase rates of a single well adjacent in the array. + void applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, + const std::vector& well_surfacerates_phase); private: // Disable copying and assignment. From f2b847bb2c5ab12f77235c81283abf0218eb3aa9 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 13:15:02 +0200 Subject: [PATCH 10/22] Added explicit reinjection to spu_2p. --- examples/spu_2p.cpp | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index b7831ae7..c0cef127 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -517,11 +517,12 @@ main(int argc, char** argv) std::vector well_bhp; std::vector well_perfrates; std::vector fractional_flows; - std::vector well_resflows; + std::vector well_resflows_phase; if (wells->c_wells()) { const int nw = wells->c_wells()->number_of_wells; well_bhp.resize(nw, 0.0); well_perfrates.resize(wells->c_wells()->well_connpos[nw], 0.0); + well_resflows_phase.resize((wells->c_wells()->number_of_phases)*(wells->c_wells()->number_of_wells), 0.0); wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates); } for (; !simtimer.done(); ++simtimer) { @@ -544,6 +545,9 @@ main(int argc, char** argv) if (check_well_controls) { computeFractionalFlow(*props, allcells, state.saturation(), fractional_flows); } + if (check_well_controls) { + wells->applyExplicitReinjectionControls(well_resflows_phase, well_resflows_phase); + } bool well_control_passed = !check_well_controls; int well_control_iteration = 0; do { @@ -583,10 +587,10 @@ main(int argc, char** argv) Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), fractional_flows, well_perfrates, - well_resflows); + well_resflows_phase); std::cout << "Checking well conditions." << std::endl; // For testing we set surface := reservoir - well_control_passed = wells->conditionsMet(well_bhp, well_resflows, well_resflows); + well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase); ++well_control_iteration; if (!well_control_passed && well_control_iteration > max_well_control_iterations) { THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); @@ -598,6 +602,7 @@ main(int argc, char** argv) } } } while (!well_control_passed); + // Process transport sources (to include bdy terms and well flows). Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, wells->c_wells(), well_perfrates, reorder_src); From 82f0128a819b128a5379557dccc33ec0bb9f0d1a Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 13:15:53 +0200 Subject: [PATCH 11/22] Multiplied with reinjection_fraction_target, as is proper. --- opm/core/WellsGroup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 13caadaa..5dcac228 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -577,7 +577,7 @@ namespace Opm // as that would check if we're under group control, something we're not. const double children_guide_rate = children_[i]->injectionGuideRate(true); children_[i]->applyInjGroupControl(InjectionSpecification::RATE, - (children_guide_rate / my_guide_rate) * total_produced, + (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, false); } } From b7f16fa73279907efa35f23482894cd4a483bc08 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 13:58:25 +0200 Subject: [PATCH 12/22] Introduced a hack to make it possible to test reinjection without supporting surface rate controls. --- opm/core/WellsGroup.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 5dcac228..7b54e62f 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -576,9 +576,15 @@ namespace Opm // Note, we do _not_ want to call the applyProdGroupControl in this object, // as that would check if we're under group control, something we're not. const double children_guide_rate = children_[i]->injectionGuideRate(true); +#ifdef DIRTY_WELLCTRL_HACK + children_[i]->applyInjGroupControl(InjectionSpecification::RESV, + (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, + false); +#else children_[i]->applyInjGroupControl(InjectionSpecification::RATE, (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, false); +#endif } } } From 322801dfb53da4cdb9a98e6cd8de7b179a599150 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 14:26:30 +0200 Subject: [PATCH 13/22] Added explicit number of ControlMode. --- opm/core/ProductionSpecification.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/ProductionSpecification.hpp b/opm/core/ProductionSpecification.hpp index ff5ace8c..dd5bc7c0 100644 --- a/opm/core/ProductionSpecification.hpp +++ b/opm/core/ProductionSpecification.hpp @@ -11,7 +11,7 @@ namespace Opm enum ControlMode { - NONE, ORAT, WRAT, GRAT, LRAT, CRAT, RESV, PRBL, BHP, THP, GRUP, FLD + NONE = 0, ORAT = 1, WRAT=2, GRAT=3, LRAT=4, CRAT=5, RESV=6, PRBL=7, BHP=8, THP=9, GRUP=10, FLD=11 }; enum Procedure From 933d4adf441b265ed2d29d5d96b37eede3757d9f Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 9 May 2012 14:49:23 +0200 Subject: [PATCH 14/22] Modified spu_2p to agree with newton compressible pressure solver. --- examples/spu_2p.cpp | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 5c69dfb7..8f11e235 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -539,25 +539,39 @@ main(int argc, char** argv) if (rock_comp->isActive()) { rc.resize(num_cells); std::vector initial_pressure = state.pressure(); + std::vector initial_porevolume(num_cells); + computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume); + std::vector pressure_increment(num_cells); std::vector prev_pressure; - for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { - prev_pressure = state.pressure(); - for (int cell = 0; cell < num_cells; ++cell) { - rc[cell] = rock_comp->rockComp(state.pressure()[cell]); - } - state.pressure() = initial_pressure; - psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(), - state.pressure(), state.faceflux(), well_bhp, well_perfrates); + + + for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { + + for (int cell = 0; cell < num_cells; ++cell) { + rc[cell] = rock_comp->rockComp(state.pressure()[cell]); + } + computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + prev_pressure = state.pressure(); + + // compute pressure increment + psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, + prev_pressure, initial_porevolume, simtimer.currentStepLength(), + pressure_increment); + double max_change = 0.0; for (int cell = 0; cell < num_cells; ++cell) { - max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell])); + state.pressure()[cell] += pressure_increment[cell]; + max_change = std::max(max_change, std::fabs(pressure_increment[cell])); } + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; if (max_change < nl_pressure_tolerance) { break; } } - computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), + well_bhp, well_perfrates); + } else { psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), well_bhp, well_perfrates); From dbb34c73a60ec98c13621045b84f4046862e1174 Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 9 May 2012 15:06:13 +0200 Subject: [PATCH 15/22] Added Newton compressible fluid solver. --- opm/core/pressure/IncompTpfa.cpp | 100 ++++++++++++++++++++++++++++++ opm/core/pressure/IncompTpfa.hpp | 22 +++++++ opm/core/pressure/tpfa/ifs_tpfa.c | 53 ++++++++++++++++ opm/core/pressure/tpfa/ifs_tpfa.h | 17 +++++ 4 files changed, 192 insertions(+) diff --git a/opm/core/pressure/IncompTpfa.cpp b/opm/core/pressure/IncompTpfa.cpp index cd0a07b6..36bb2fcb 100644 --- a/opm/core/pressure/IncompTpfa.cpp +++ b/opm/core/pressure/IncompTpfa.cpp @@ -26,6 +26,7 @@ #include #include #include +#include namespace Opm { @@ -233,9 +234,108 @@ namespace Opm soln.well_flux = &well_rate[0]; soln.well_press = &well_bhp[0]; } + ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln); } + + void IncompTpfa::solveIncrement(const std::vector& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + const std::vector& porevol, + const std::vector& rock_comp, + const std::vector& prev_pressure, + const std::vector& initial_porevol, + const double dt, + std::vector& pressure_increment) + { + UnstructuredGrid* gg = const_cast(&grid_); + tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]); + + if (!omega.empty()) { + if (gpress_.empty()) { + THROW("Nozero omega argument given, but gravity was null in constructor."); + } + mim_ip_density_update(gg->number_of_cells, gg->cell_facepos, + &omega[0], + &gpress_[0], &gpress_omegaweighted_[0]); + } else { + if (!gpress_.empty()) { + THROW("Empty omega argument given, but gravity was non-null in constructor."); + } + } + + ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL }; + if (! src.empty()) { F.src = &src[0]; } + F.bc = bcs; + F.totmob = &totmob[0]; + F.wdp = &wdp[0]; + + if (rock_comp.empty()) { + ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_); + } else { + ifs_tpfa_assemble_comprock_increment(gg, &F, &trans_[0], &gpress_omegaweighted_[0], + &porevol[0], &rock_comp[0], dt, &prev_pressure[0], + &initial_porevol[0], h_); + } + + linsolver_.solve(h_->A, h_->b, &pressure_increment[0]); + + } + + void IncompTpfa::computeFaceFlux(const std::vector& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + std::vector& pressure, + std::vector& faceflux, + std::vector& well_bhp, + std::vector& well_rate) { + + UnstructuredGrid* gg = const_cast(&grid_); + tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]); + + if (!omega.empty()) { + if (gpress_.empty()) { + THROW("Nozero omega argument given, but gravity was null in constructor."); + } + mim_ip_density_update(gg->number_of_cells, gg->cell_facepos, + &omega[0], + &gpress_[0], &gpress_omegaweighted_[0]); + } else { + if (!gpress_.empty()) { + THROW("Empty omega argument given, but gravity was non-null in constructor."); + } + } + + ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL }; + if (! src.empty()) { F.src = &src[0]; } + F.bc = bcs; + F.totmob = &totmob[0]; + F.wdp = &wdp[0]; + + ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_); + + faceflux.resize(grid_.number_of_faces); + + ifs_tpfa_solution soln = { NULL, NULL, NULL, NULL }; + soln.cell_press = &pressure[0]; + soln.face_flux = &faceflux[0]; + + if(wells_ != NULL) { + well_bhp.resize(wells_->number_of_wells); + well_rate.resize(wells_->well_connpos[ wells_->number_of_wells ]); + soln.well_flux = &well_rate[0]; + soln.well_press = &well_bhp[0]; + } + + memcpy(h_->x, &pressure[0], grid_.number_of_cells * sizeof *(h_->x)); + + ifs_tpfa_press_flux(gg, &F, &trans_[0], h_, &soln); + } } // namespace Opm diff --git a/opm/core/pressure/IncompTpfa.hpp b/opm/core/pressure/IncompTpfa.hpp index f40d3e23..f14d9108 100644 --- a/opm/core/pressure/IncompTpfa.hpp +++ b/opm/core/pressure/IncompTpfa.hpp @@ -121,6 +121,28 @@ namespace Opm std::vector& well_bhp, std::vector& well_rate); + void solveIncrement(const std::vector& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + const std::vector& porevol, + const std::vector& rock_comp, + const std::vector& prev_pressure, + const std::vector& initial_porevol, + const double dt, + std::vector& pressure_increment); + + void computeFaceFlux(const std::vector& totmob, + const std::vector& omega, + const std::vector& src, + const std::vector& wdp, + const FlowBoundaryConditions* bcs, + std::vector& pressure, + std::vector& faceflux, + std::vector& well_bhp, + std::vector& well_rate); + /// Expose read-only reference to internal half-transmissibility. const ::std::vector& getHalfTrans() const { return htrans_; } diff --git a/opm/core/pressure/tpfa/ifs_tpfa.c b/opm/core/pressure/tpfa/ifs_tpfa.c index c3566c93..22a6e482 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.c +++ b/opm/core/pressure/tpfa/ifs_tpfa.c @@ -11,6 +11,16 @@ #include +void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v) { + size_t i,j; + for (j = 0; j < A->m; ++j) { + v[j] = 0; + for (i = (size_t) (A->ia[j]); i < (size_t) (A->ia[j+1]); ++i) { + v[j] += A->sa[i]*u[A->ja[i]]; + } + } +} + struct ifs_tpfa_impl { double *fgrav; /* Accumulated grav contrib/face */ @@ -725,6 +735,49 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid *G , return ok; } +/* ---------------------------------------------------------------------- */ +int +ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G , + const struct ifs_tpfa_forces *F , + const double *trans , + const double *gpress , + const double *porevol , + const double *rock_comp, + const double dt , + const double *prev_pressure, + const double *initial_porevolume, + struct ifs_tpfa_data *h ) +/* ---------------------------------------------------------------------- */ +{ + int c; + size_t j; + int system_singular, ok; + double *v; + + assemble_incompressible(G, F, trans, gpress, h, &system_singular, &ok); + + v = malloc(h->A->m * sizeof *v); + + mult_csr_matrix(h->A, prev_pressure, v); + + /* We want to solve a Newton step for the residual + * (porevol(pressure)-porevol(initial_pressure))/dt + residual_for_imcompressible + * + */ + + if (ok) { + for (c = 0; c < G->number_of_cells; ++c) { + j = csrmatrix_elm_index(c, c, h->A); + h->A->sa[j] += porevol[c] * rock_comp[c] / dt; + h->b[c] += -(porevol[c] - initial_porevolume[c])/dt - v[c]; + } + } + + free(v); + + return ok; +} + /* ---------------------------------------------------------------------- */ void diff --git a/opm/core/pressure/tpfa/ifs_tpfa.h b/opm/core/pressure/tpfa/ifs_tpfa.h index 8f8aed75..c06c2f66 100644 --- a/opm/core/pressure/tpfa/ifs_tpfa.h +++ b/opm/core/pressure/tpfa/ifs_tpfa.h @@ -61,6 +61,10 @@ struct ifs_tpfa_data * ifs_tpfa_construct(struct UnstructuredGrid *G, struct Wells *W); + +void mult_csr_matrix(const struct CSRMatrix* A, const double* u, double* v); + + int ifs_tpfa_assemble(struct UnstructuredGrid *G , const struct ifs_tpfa_forces *F , @@ -78,6 +82,19 @@ ifs_tpfa_assemble_comprock(struct UnstructuredGrid *G , const double dt , const double *pressure , struct ifs_tpfa_data *h ); +int +ifs_tpfa_assemble_comprock_increment(struct UnstructuredGrid *G , + const struct ifs_tpfa_forces *F , + const double *trans , + const double *gpress , + const double *porevol , + const double *rock_comp, + const double dt , + const double *prev_pressure , + const double *initial_porevolume, + struct ifs_tpfa_data *h ); + + void ifs_tpfa_press_flux(struct UnstructuredGrid *G , const struct ifs_tpfa_forces *F , From 9dabbd349f3c0688f22e4c4cec8ff0810087bd1a Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Wed, 9 May 2012 15:19:37 +0200 Subject: [PATCH 16/22] Added Newton compressible pressure solver in spu_2p. --- examples/spu_2p.cpp | 109 ++++++++++++++++++++++++-------------------- 1 file changed, 60 insertions(+), 49 deletions(-) diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index c0cef127..daa68b65 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -555,58 +555,69 @@ main(int argc, char** argv) if (rock_comp->isActive()) { rc.resize(num_cells); std::vector initial_pressure = state.pressure(); + std::vector initial_porevolume(num_cells); + computePorevolume(*grid->c_grid(), *props, *rock_comp, initial_pressure, initial_porevolume); + std::vector pressure_increment(num_cells); std::vector prev_pressure; - for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { - prev_pressure = state.pressure(); - for (int cell = 0; cell < num_cells; ++cell) { - rc[cell] = rock_comp->rockComp(state.pressure()[cell]); - } - state.pressure() = initial_pressure; - psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, simtimer.currentStepLength(), - state.pressure(), state.faceflux(), well_bhp, well_perfrates); - double max_change = 0.0; - for (int cell = 0; cell < num_cells; ++cell) { - max_change = std::max(max_change, std::fabs(state.pressure()[cell] - prev_pressure[cell])); - } - std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; - if (max_change < nl_pressure_tolerance) { - break; - } - } - computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); - } else { - psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), - well_bhp, well_perfrates); - } - pressure_timer.stop(); - double pt = pressure_timer.secsSinceStart(); - std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; - ptime += pt; + for (int iter = 0; iter < nl_pressure_maxiter; ++iter) { + + for (int cell = 0; cell < num_cells; ++cell) { + rc[cell] = rock_comp->rockComp(state.pressure()[cell]); + } + computePorevolume(*grid->c_grid(), *props, *rock_comp, state.pressure(), porevol); + prev_pressure = state.pressure(); + + // compute pressure increment + psolver.solveIncrement(totmob, omega, src, wdp, bcs.c_bcs(), porevol, rc, + prev_pressure, initial_porevolume, simtimer.currentStepLength(), + pressure_increment); - if (check_well_controls) { - Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), - fractional_flows, - well_perfrates, - well_resflows_phase); - std::cout << "Checking well conditions." << std::endl; - // For testing we set surface := reservoir - well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase); - ++well_control_iteration; - if (!well_control_passed && well_control_iteration > max_well_control_iterations) { - THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); - } - if (!well_control_passed) { - std::cout << "Well controls not passed, solving again." << std::endl; - } else { - std::cout << "Well conditions met." << std::endl; - } - } - } while (!well_control_passed); + double max_change = 0.0; + for (int cell = 0; cell < num_cells; ++cell) { + state.pressure()[cell] += pressure_increment[cell]; + max_change = std::max(max_change, std::fabs(pressure_increment[cell])); + } + + std::cout << "Pressure iter " << iter << " max change = " << max_change << std::endl; + if (max_change < nl_pressure_tolerance) { + break; + } + } + psolver.computeFaceFlux(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), + well_bhp, well_perfrates); + } else { + psolver.solve(totmob, omega, src, wdp, bcs.c_bcs(), state.pressure(), state.faceflux(), + well_bhp, well_perfrates); + } + pressure_timer.stop(); + double pt = pressure_timer.secsSinceStart(); + std::cout << "Pressure solver took: " << pt << " seconds." << std::endl; + ptime += pt; + + if (check_well_controls) { + Opm::computePhaseFlowRatesPerWell(*wells->c_wells(), + fractional_flows, + well_perfrates, + well_resflows_phase); + std::cout << "Checking well conditions." << std::endl; + // For testing we set surface := reservoir + well_control_passed = wells->conditionsMet(well_bhp, well_resflows_phase, well_resflows_phase); + ++well_control_iteration; + if (!well_control_passed && well_control_iteration > max_well_control_iterations) { + THROW("Could not satisfy well conditions in " << max_well_control_iterations << " tries."); + } + if (!well_control_passed) { + std::cout << "Well controls not passed, solving again." << std::endl; + } else { + std::cout << "Well conditions met." << std::endl; + } + } + } while (!well_control_passed); - // Process transport sources (to include bdy terms and well flows). - Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, - wells->c_wells(), well_perfrates, reorder_src); - if (!use_reorder) { + // Process transport sources (to include bdy terms and well flows). + Opm::computeTransportSource(*grid->c_grid(), src, state.faceflux(), 1.0, + wells->c_wells(), well_perfrates, reorder_src); + if (!use_reorder) { clear_transport_source(tsrc); for (int cell = 0; cell < num_cells; ++cell) { if (reorder_src[cell] > 0.0) { From 322831d0b176ed9bc21a137da4d9fa028d172dc0 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 15:43:43 +0200 Subject: [PATCH 17/22] Added initial support for VREP. --- opm/core/WellsGroup.cpp | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 7b54e62f..bd469473 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -553,7 +553,7 @@ namespace Opm /// A vector containing surface rates by phase for each well. /// Is assumed to be ordered the same way as the related Wells-struct, /// with all phase rates of a single well adjacent in the array. - void WellsGroup::applyExplicitReinjectionControls(const std::vector& /*well_reservoirrates_phase*/, + void WellsGroup::applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase) { if (injSpec().control_mode_ == InjectionSpecification::REIN) { @@ -587,6 +587,30 @@ namespace Opm #endif } } + else if (injSpec().control_mode_ == InjectionSpecification::VREP) { + double total_produced = 0.0; + if (phaseUsage().phase_used[BlackoilPhases::Aqua]) { + total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Aqua); + } + if (phaseUsage().phase_used[BlackoilPhases::Liquid]) { + total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Liquid); + } + if (phaseUsage().phase_used[BlackoilPhases::Vapour]) { + total_produced += getTotalProductionFlow(well_reservoirrates_phase, BlackoilPhases::Vapour); + } + + const double my_guide_rate = injectionGuideRate(true); + for (size_t i = 0; i < children_.size(); ++i) { + // Apply for all children. + // Note, we do _not_ want to call the applyProdGroupControl in this object, + // as that would check if we're under group control, something we're not. + const double children_guide_rate = children_[i]->injectionGuideRate(true); + children_[i]->applyInjGroupControl(InjectionSpecification::RESV, + (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, + false); + } + + } } // ============== WellNode members ============ From fba33e03d80b599e3eb8e0f14ed73e318e4d2452 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 15:47:01 +0200 Subject: [PATCH 18/22] Reading the voidage replacement fraction from GCONINJE. --- opm/core/eclipse/SpecialEclipseFields.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index 5c374e7f..62e4c9cd 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -1001,6 +1001,7 @@ struct GconinjeLine double surface_flow_max_rate_; // Surface flow rate target or upper limit double resv_flow_max_rate_; // Reservoir flow rate target or upper limit double reinjection_fraction_target_; + double voidage_replacement_fraction_; // Default values GconinjeLine() : @@ -1048,6 +1049,7 @@ struct GCONINJE : public SpecialBase gconinje_line.surface_flow_max_rate_ = double_data[0]; gconinje_line.resv_flow_max_rate_ = double_data[1]; gconinje_line.reinjection_fraction_target_ = double_data[2]; + gconinje_line.voidage_replacement_fraction_ = double_data[3]; // HACK! Ignore any further items if (num_read == num_to_read) { ignoreSlashLine(is); From 3bb9fcb58dd4b175fe43996812b0e74bd9f531e2 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 15:48:27 +0200 Subject: [PATCH 19/22] Corrected some default values in GCONINJE. --- opm/core/eclipse/SpecialEclipseFields.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index 62e4c9cd..269715ab 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -1007,7 +1007,8 @@ struct GconinjeLine GconinjeLine() : surface_flow_max_rate_(-1.0E20), resv_flow_max_rate_(-1E20), - reinjection_fraction_target_(-1E20) + reinjection_fraction_target_(1), + voidage_replacement_fraction_(1) { } }; From a0ad1f0acfcb4e3caa1f6ffb222e97ba7732e98d Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 15:54:25 +0200 Subject: [PATCH 20/22] Included voidage_replacment in injection specification. --- opm/core/InjectionSpecification.cpp | 3 ++- opm/core/InjectionSpecification.hpp | 1 + opm/core/WellsGroup.cpp | 2 ++ 3 files changed, 5 insertions(+), 1 deletion(-) diff --git a/opm/core/InjectionSpecification.cpp b/opm/core/InjectionSpecification.cpp index 691830c6..dca7fc76 100644 --- a/opm/core/InjectionSpecification.cpp +++ b/opm/core/InjectionSpecification.cpp @@ -8,7 +8,8 @@ namespace Opm surface_flow_max_rate_(-1e100), reservoir_flow_max_rate_(-1e100), BHP_limit_(-1e100), - reinjection_fraction_target_(-1e100), + reinjection_fraction_target_(1), + voidage_replacment_fraction_(1), guide_rate_(1.0), guide_rate_type_(NONE_GRT) { diff --git a/opm/core/InjectionSpecification.hpp b/opm/core/InjectionSpecification.hpp index 62c9e667..107296c6 100644 --- a/opm/core/InjectionSpecification.hpp +++ b/opm/core/InjectionSpecification.hpp @@ -32,6 +32,7 @@ namespace Opm double reservoir_flow_max_rate_; double BHP_limit_; double reinjection_fraction_target_; + double voidage_replacment_fraction_; double guide_rate_; GuideRateType guide_rate_type_; }; diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index bd469473..bd14e69d 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -1099,6 +1099,8 @@ namespace Opm injection_specification.control_mode_ = toInjectionControlMode(line.control_mode_); injection_specification.surface_flow_max_rate_ = line.surface_flow_max_rate_; injection_specification.reservoir_flow_max_rate_ = line.resv_flow_max_rate_; + injection_specification.reinjection_fraction_target_ = line.reinjection_fraction_target_; + injection_specification.voidage_replacment_fraction_ = line.voidage_replacement_fraction_; } } } From 26abaaf484c91495f8abc27ddaf56ec8512184b5 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 15:55:59 +0200 Subject: [PATCH 21/22] Used voidage_replacment correctly in applyExplicitReinjectionControls --- opm/core/WellsGroup.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index bd14e69d..22874090 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -606,7 +606,7 @@ namespace Opm // as that would check if we're under group control, something we're not. const double children_guide_rate = children_[i]->injectionGuideRate(true); children_[i]->applyInjGroupControl(InjectionSpecification::RESV, - (children_guide_rate / my_guide_rate) * total_produced * injSpec().reinjection_fraction_target_, + (children_guide_rate / my_guide_rate) * total_produced * injSpec().voidage_replacment_fraction_, false); } From 9f651f255c68aa6b43f8fd99a1d30c5078fbb5c1 Mon Sep 17 00:00:00 2001 From: Kjetil Olsen Lye Date: Wed, 9 May 2012 16:03:21 +0200 Subject: [PATCH 22/22] Added correct handling of VREP and REIN keyword. No longer generates a warning, just a friendly hint. --- opm/core/WellsGroup.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 22874090..006b7610 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -489,8 +489,9 @@ namespace Opm } return; } + case InjectionSpecification::VREP: case InjectionSpecification::REIN: - std::cout << "WARNING: Ignoring control type REIN" << std::endl; + std::cout << "Replacement keywords found, remember to call applyExplicitReinjectionControls." << std::endl; return; case InjectionSpecification::FLD: case InjectionSpecification::NONE: