From a214c10595c37b8a184df7aedc6037dbcc9c6215 Mon Sep 17 00:00:00 2001 From: Joakim Hove Date: Wed, 17 Feb 2016 10:04:13 +0100 Subject: [PATCH] Changes in SimulatorState: 1. Added method setCellDataComponent() 2. Removed setFirstSat() Implemented saturation initialisation using setCellDataComponent() instead of setFirstSat(). This way the template has been removed from the SimulatorState class. --- opm/core/simulator/BlackoilState.cpp | 10 ---- opm/core/simulator/BlackoilState.hpp | 8 --- opm/core/simulator/SimulatorState.cpp | 64 ++++++++++------------- opm/core/simulator/SimulatorState.hpp | 19 ++----- opm/core/simulator/TwophaseState.hpp | 4 -- opm/core/simulator/TwophaseState_impl.hpp | 7 --- opm/core/simulator/initState_impl.hpp | 46 +++++++++++++--- tutorials/tutorial3.cpp | 14 ++++- tutorials/tutorial4.cpp | 15 +++++- 9 files changed, 100 insertions(+), 87 deletions(-) diff --git a/opm/core/simulator/BlackoilState.cpp b/opm/core/simulator/BlackoilState.cpp index 46b9d3a5..5084a1fd 100644 --- a/opm/core/simulator/BlackoilState.cpp +++ b/opm/core/simulator/BlackoilState.cpp @@ -23,16 +23,6 @@ BlackoilState::init(const UnstructuredGrid& g, int num_phases) { init(g.number_of_cells, g.number_of_faces, num_phases); } -/// Set the first saturation to either its min or max value in -/// the indicated cells. The second saturation value s2 is set -/// to (1.0 - s1) for each cell. Any further saturation values -/// are unchanged. -void -BlackoilState::setFirstSat(const std::vector& cells, - const Opm::BlackoilPropertiesInterface& props, - ExtremalSat es) { - SimulatorState::setFirstSat(cells, props, es); -} bool BlackoilState::equals(const SimulatorState& other, diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp index b6622e19..ce7655a8 100644 --- a/opm/core/simulator/BlackoilState.hpp +++ b/opm/core/simulator/BlackoilState.hpp @@ -39,14 +39,6 @@ namespace Opm virtual void init(int number_of_cells, int number_of_faces, int num_phases); - /// Set the first saturation to either its min or max value in - /// the indicated cells. The second saturation value s2 is set - /// to (1.0 - s1) for each cell. Any further saturation values - /// are unchanged. - void setFirstSat(const std::vector& cells, - const Opm::BlackoilPropertiesInterface& props, - ExtremalSat es); - virtual bool equals(const SimulatorState& other, double epsilon = 1e-8) const; diff --git a/opm/core/simulator/SimulatorState.cpp b/opm/core/simulator/SimulatorState.cpp index ccd9864a..45d7d3c8 100644 --- a/opm/core/simulator/SimulatorState.cpp +++ b/opm/core/simulator/SimulatorState.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -80,42 +81,35 @@ SimulatorState::registerFaceData( const std::string& name, const int components, return pos ; } -template void -SimulatorState::setFirstSat(const std::vector& cells, - const Props& props, - ExtremalSat es) -{ - if (cells.empty()) { - return; - } - int n = cells.size(); - std::vector smin(num_phases_*n); - std::vector smax(num_phases_*n); - props.satRange(n, &cells[0], &smin[0], &smax[0]); - std::vector< double >& sat = saturation(); - const double* svals = (es == MinSat) ? &smin[0] : &smax[0]; - for (int ci = 0; ci < n; ++ci) { - const int cell = cells[ci]; - sat[num_phases_*cell] = svals[num_phases_*ci]; - sat[num_phases_*cell + 1] = 1.0 - sat[num_phases_*cell]; +void SimulatorState::setCellDataComponent( const std::string& name , size_t component , const std::vector& cells , const std::vector& values) { + const auto iter = std::find( cellDataNames_.begin() , cellDataNames_.end() , name); + int id = iter - cellDataNames_.begin(); + auto& data = cellData_[id]; + if (component >= num_phases_) + throw std::invalid_argument("Invalid component"); + + if (cells.size() != values.size()) + throw std::invalid_argument("size mismatch between cells and values"); + + /* This is currently quite broken; the setCellDataComponent + method assumes that the number of components in the field + we are currently focusing on has num_phases components in + total. This restriction should be lifted by allowing a per + field number of components. + */ + if (data.size() != num_phases_ * num_cells_) + throw std::invalid_argument("Can currently only be used on fields with num_components == num_phases (i.e. saturation...) "); + + for (size_t i = 0; i < cells.size(); i++) { + if (cells[i] < num_cells_) { + auto field_index = cells[i] * num_phases_ + component; + auto value = values[i]; + + data[field_index] = value; + } else { + throw std::invalid_argument("Invalid cell number"); } + } } -// template instantiations for all known (to this library) subclasses -// of SimulatorState that will call this method. notice that there are -// no empty angle brackets after "template" -- that would have been -// specialization instead -#include -#include -template void -SimulatorState::setFirstSat ( - const std::vector &cells, - const IncompPropertiesInterface &props, - ExtremalSat es); - -template void -SimulatorState::setFirstSat ( - const std::vector &cells, - const BlackoilPropertiesInterface &props, - ExtremalSat es); diff --git a/opm/core/simulator/SimulatorState.hpp b/opm/core/simulator/SimulatorState.hpp index 4ef06d95..e3dd9d3c 100644 --- a/opm/core/simulator/SimulatorState.hpp +++ b/opm/core/simulator/SimulatorState.hpp @@ -17,8 +17,6 @@ namespace Opm virtual void init(int number_of_cells, int number_of_faces, int num_phases); - enum ExtremalSat { MinSat, MaxSat }; - protected: /// \brief pressure per cell. static const int pressureId_ = 0; @@ -32,19 +30,12 @@ namespace Opm /// \brief The fluxes at the faces. static const int faceFluxId_ = 1; - /** - * Initialize the first saturation to maximum value. This method - * should be considered deprecated. Avoid to use it! - * - * \tparam Props Fluid and rock properties that pertain to this - * kind of simulation. Currently, only Blackoil- - * and IncompPropertiesInterface are supported. - */ - template - void setFirstSat(const std::vector& cells, - const Props& props, - ExtremalSat es); public: + /// Will set the values of component nr @component in the + /// field @key. All the cells in @cells will be set to the + /// values in @values. + void setCellDataComponent( const std::string& key , size_t component , const std::vector& cells , const std::vector& values); + int numPhases() const { return num_phases_; } int numCells () const { return num_cells_; } int numFaces () const { return num_faces_; } diff --git a/opm/core/simulator/TwophaseState.hpp b/opm/core/simulator/TwophaseState.hpp index 21e51d5c..c8bcf23b 100644 --- a/opm/core/simulator/TwophaseState.hpp +++ b/opm/core/simulator/TwophaseState.hpp @@ -30,10 +30,6 @@ namespace Opm class TwophaseState : public SimulatorState { public: - void setFirstSat(const std::vector& cells, - const Opm::IncompPropertiesInterface& props, - ExtremalSat es); - virtual bool equals (const SimulatorState& other, double epsilon = 1e-8) const; }; diff --git a/opm/core/simulator/TwophaseState_impl.hpp b/opm/core/simulator/TwophaseState_impl.hpp index d54ecae5..6dbf9b36 100644 --- a/opm/core/simulator/TwophaseState_impl.hpp +++ b/opm/core/simulator/TwophaseState_impl.hpp @@ -4,13 +4,6 @@ namespace Opm { -inline void -TwophaseState::setFirstSat(const std::vector& cells, - const Opm::IncompPropertiesInterface& props, - ExtremalSat es) { - SimulatorState::setFirstSat(cells, props, es); -} - inline bool TwophaseState::equals (const SimulatorState& other, double epsilon) const { diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index d411d613..8134b8b9 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -76,6 +76,34 @@ namespace Opm #endif /* __clang__ */ enum WaterInit { WaterBelow, WaterAbove }; + enum ExtremalSat { MinSat, MaxSat }; + + /// Will initialize the first and second component of the + /// SATURATION field in all the cells in the set @cells. The + /// @props object will be queried, and depending on the value + /// @satType either the minimum or the maximum saturation is + /// applied to thee first component in the SATURATION field. + /// For the second component (1 - first_sat) is used. + + template + static void initSaturation(const std::vector& cells , const Props& props , State& state , ExtremalSat satType) { + std::vector min_sat(cells.size()); + std::vector max_sat(cells.size()); + std::vector second_sat(cells.size()); + std::vector* first_sat; + + props.satRange(cells.size() ,cells.data() , min_sat.data() , max_sat.data()); + if (satType == MinSat) { + first_sat = &min_sat; + } else { + first_sat = &max_sat; + } + + std::transform( first_sat->begin() , first_sat->end() , second_sat.begin() , [](double s) { return 1 - s; }); + state.setCellDataComponent( "SATURATION" , 0 , cells , *first_sat ); + state.setCellDataComponent( "SATURATION" , 1 , cells , second_sat ); + } + // Initialize saturations so that there is water below woc, // and oil above. @@ -106,9 +134,9 @@ namespace Opm cellsBelowAbove(number_of_cells, begin_cell_centroids, dimensions, woc, oil, water); } - // Set saturations. - state.setFirstSat(oil, props, State::MinSat); - state.setFirstSat(water, props, State::MaxSat); + + initSaturation( oil , props , state , MinSat ); + initSaturation( water , props , state , MaxSat ); } @@ -380,7 +408,10 @@ namespace Opm for (int i = 0; i < num_cells; ++i) { all_cells[i] = i; } - state.setFirstSat(all_cells, props, State::MinSat); + + initSaturation( all_cells , props , state , MinSat ); + + const bool convection_testcase = param.getDefault("convection_testcase", false); const bool segregation_testcase = param.getDefault("segregation_testcase", false); if (convection_testcase) { @@ -394,7 +425,8 @@ namespace Opm left_cells.push_back(cell); } } - state.setFirstSat(left_cells, props, State::MaxSat); + + initSaturation( left_cells , props , state , MaxSat ); const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; std::fill(state.pressure().begin(), state.pressure().end(), init_p); } else if (segregation_testcase) { @@ -501,7 +533,7 @@ namespace Opm for (int i = 0; i < num_cells; ++i) { all_cells[i] = i; } - state.setFirstSat(all_cells, props, State::MinSat); + initSaturation(all_cells , props , state , MinSat); const bool convection_testcase = param.getDefault("convection_testcase", false); if (convection_testcase) { // Initialise water saturation to max in the 'left' part. @@ -514,7 +546,7 @@ namespace Opm left_cells.push_back(cell); } } - state.setFirstSat(left_cells, props, State::MaxSat); + initSaturation(left_cells , props , state , MaxSat ); const double init_p = param.getDefault("ref_pressure", 100.0)*unit::barsa; std::fill(state.pressure().begin(), state.pressure().end(), init_p); } else if (param.has("water_oil_contact")) { diff --git a/tutorials/tutorial3.cpp b/tutorials/tutorial3.cpp index 9103cb2c..7f39159f 100644 --- a/tutorials/tutorial3.cpp +++ b/tutorials/tutorial3.cpp @@ -267,7 +267,19 @@ try /// \internal [two-phase state] TwophaseState state; state.init(grid.number_of_cells , grid.number_of_faces, 2); - state.setFirstSat(allcells, props, TwophaseState::MinSat); + { + std::vector min_sat(allcells.size()); + std::vector max_sat(allcells.size()); + std::vector second_sat(allcells.size()); + + props.satRange(allcells.size() ,allcells.data() , min_sat.data() , max_sat.data()); + + std::transform( min_sat.begin() , min_sat.end() , second_sat.begin() , [](double s) { return 1 - s; }); + state.setCellDataComponent( "SATURATION" , 0 , allcells , min_sat ); + state.setCellDataComponent( "SATURATION" , 1 , allcells , second_sat ); + + } + /// \internal [two-phase state] /// \endinternal diff --git a/tutorials/tutorial4.cpp b/tutorials/tutorial4.cpp index 926cd8c4..1d8efd9f 100644 --- a/tutorials/tutorial4.cpp +++ b/tutorials/tutorial4.cpp @@ -214,7 +214,20 @@ try /// \internal[two-phase state] TwophaseState state; state.init(grid.number_of_cells , grid.number_of_faces, 2); - state.setFirstSat(allcells, props, TwophaseState::MinSat); + { + std::vector min_sat(allcells.size()); + std::vector max_sat(allcells.size()); + std::vector second_sat(allcells.size()); + + props.satRange(allcells.size() ,allcells.data() , min_sat.data() , max_sat.data()); + + std::transform( min_sat.begin() , min_sat.end() , second_sat.begin() , [](double s) { return 1 - s; }); + state.setCellDataComponent( "SATURATION" , 0 , allcells , min_sat ); + state.setCellDataComponent( "SATURATION" , 1 , allcells , second_sat ); + + } + + /// \internal[two-phase state] /// \endinternal