diff --git a/opm/core/InjectionSpecification.cpp b/opm/core/InjectionSpecification.cpp index 7ed0d225..fe6a1253 100644 --- a/opm/core/InjectionSpecification.cpp +++ b/opm/core/InjectionSpecification.cpp @@ -3,11 +3,15 @@ namespace Opm { InjectionSpecification::InjectionSpecification() - : injector_type_(WATER), control_mode_(NONE), surface_flow_max_rate_(1e100), - reinjection_fraction_target_(0.0), BHP_limit_(1e100), fluid_volume_max_rate_(1e100) + : injector_type_(WATER), + control_mode_(NONE), + surface_flow_max_rate_(1e100), + reinjection_fraction_target_(0.0), + fluid_volume_max_rate_(1e100), + BHP_limit_(1e100) { } -} // namespace Opm \ No newline at end of file +} // namespace Opm diff --git a/opm/core/WellsGroup.cpp b/opm/core/WellsGroup.cpp index 1ab1e034..6bf82b39 100644 --- a/opm/core/WellsGroup.cpp +++ b/opm/core/WellsGroup.cpp @@ -7,6 +7,7 @@ #include #include +#include namespace Opm { @@ -130,6 +131,7 @@ namespace Opm int number_of_leaf_nodes = numberOfLeafNodes(); + bool shut_down_on_exceed = false; double bhp_target = 1e100; double rate_target = 1e100; switch(wells->type[index_of_well]) { @@ -145,6 +147,7 @@ namespace Opm const ProductionSpecification& prod_spec = prodSpec(); bhp_target = prod_spec.BHP_limit_ / number_of_leaf_nodes; rate_target = prod_spec.fluid_volume_max_rate_ / number_of_leaf_nodes; + shut_down_on_exceed = prodSpec().procedure_ == ProductionSpecification::WELL; break; } } @@ -153,12 +156,26 @@ namespace Opm std::cout << "BHP not met" << std::endl; std::cout << "BHP limit was " << bhp_target << std::endl; std::cout << "Actual bhp was " << well_bhp[index_of_well] << std::endl; + + if(shut_down_on_exceed) { + // Shut down well + // Dirty hack for now + struct Wells* non_const_wells = const_cast(wells); + non_const_wells->ctrls[index_of_well]->target[0] = 0.0; + } return false; } if(well_rate[index_of_well] - rate_target > epsilon) { std::cout << "well_rate not met" << std::endl; std::cout << "target = " << rate_target << ", well_rate[index_of_well] = " << well_rate[index_of_well] << std::endl; std::cout << "Group name = " << name() << std::endl; + + if(shut_down_on_exceed) { + // Shut down well + // Dirty hack for now + struct Wells* non_const_wells = const_cast(wells); + non_const_wells->ctrls[index_of_well]->target[0] = 0.0; + } return false; } return true; diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index 15afd936..360acc60 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -1004,8 +1004,9 @@ struct GconinjeLine // Default values GconinjeLine() : - surface_flow_max_rate_(1.0E20), reinjection_fraction_target_(1E20), - resv_flow_max_rate_(1E20) + surface_flow_max_rate_(1.0E20), + resv_flow_max_rate_(1E20), + reinjection_fraction_target_(1E20) { } }; diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index 5dfbbdad..2ac7f8c5 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/opm/core/utility/miscUtilities.cpp @@ -406,12 +406,12 @@ namespace Opm void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector& saturations, - const std::vector& densities, std::vector& wdp) + const std::vector& densities, std::vector& wdp, bool per_grid_cell) { const size_t np = densities.size(); - + const int nw = wells.number_of_wells; // Simple for now: - for(int i = 0; i < wells.number_of_wells; i++) { + for(int i = 0; i < nw; i++) { double depth_ref = wells.depth_ref[i]; for(int j = wells.well_connpos[i]; j < wells.well_connpos[i+1]; j++) { int cell = wells.well_cells[j]; @@ -421,15 +421,25 @@ namespace Opm double saturation_sum = 0.0; for(size_t p = 0; p < np; p++) { - saturation_sum += saturations[np*cell + p]; + if(per_grid_cell) { + saturation_sum += saturations[i*nw*np + j*np + p]; + } + else { + saturation_sum += saturations[np*cell + p]; + } } if(saturation_sum == 0) { saturation_sum = 1.0; } double density = 0.0; for(size_t p = 0; p < np; p++) { - // Is this a smart way of doing it? - density += saturations[np*cell + p] * densities[p] / saturation_sum; + if(per_grid_cell) { + density += saturations[i*nw*np + j*np + p] * densities[p] / saturation_sum; + } + else { + // Is this a smart way of doing it? + density += saturations[np*cell + p] * densities[p] / saturation_sum; + } } // Is the sign correct? diff --git a/opm/core/utility/miscUtilities.hpp b/opm/core/utility/miscUtilities.hpp index 4334aba9..40e0c6be 100644 --- a/opm/core/utility/miscUtilities.hpp +++ b/opm/core/utility/miscUtilities.hpp @@ -180,13 +180,15 @@ namespace Opm /// \param[in] wells Wells that need their wdp calculated. /// \param[in] grid The associated grid to make cell lookups. /// \param[in] saturations A vector of weights for each cell for each phase - /// in the grid. So for cell i, + /// in the grid (or well, see per_grid_cell parameter). So for cell i, /// saturations[i*densities.size() + p] should give the weight /// of phase p in cell i. /// \param[in] densities Density for each phase. /// \param[out] wdp Will contain, for each well, the wdp of the well. + /// \param[in] per_grid_cell Whether or not the saturations are per grid cell or per + /// well cell. void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector& saturations, - const std::vector& densities, std::vector& wdp); + const std::vector& densities, std::vector& wdp, bool per_grid_cell = true); /// Computes (sums) the flow rate for each well. /// \param[in] wells The wells for which the flow rate should be computed.