diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 89c5da73..4677762e 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -95,9 +95,11 @@ list (APPEND MAIN_SOURCE_FILES opm/core/props/satfunc/SatFuncStone2.cpp opm/core/props/satfunc/SaturationPropsBasic.cpp opm/core/props/satfunc/SaturationPropsFromDeck.cpp + opm/core/simulator/BlackoilState.cpp opm/core/simulator/SimulatorCompressibleTwophase.cpp opm/core/simulator/SimulatorIncompTwophase.cpp opm/core/simulator/SimulatorReport.cpp + opm/core/simulator/SimulatorState.cpp opm/core/simulator/SimulatorTimer.cpp opm/core/tof/DGBasis.cpp opm/core/tof/TofReorder.cpp @@ -304,8 +306,10 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/simulator/SimulatorCompressibleTwophase.hpp opm/core/simulator/SimulatorIncompTwophase.hpp opm/core/simulator/SimulatorReport.hpp + opm/core/simulator/SimulatorState.hpp opm/core/simulator/SimulatorTimer.hpp opm/core/simulator/TwophaseState.hpp + opm/core/simulator/TwophaseState_impl.hpp opm/core/simulator/WellState.hpp opm/core/simulator/initState.hpp opm/core/simulator/initState_impl.hpp diff --git a/opm/core/simulator/BlackoilState.cpp b/opm/core/simulator/BlackoilState.cpp new file mode 100644 index 00000000..c8c3afa8 --- /dev/null +++ b/opm/core/simulator/BlackoilState.cpp @@ -0,0 +1,37 @@ +#include "BlackoilState.hpp" +#include + +using namespace Opm; + +void +BlackoilState::init(const UnstructuredGrid& g, int num_phases) { + SimulatorState::init(g, num_phases); + gor_.resize(g.number_of_cells, 0.) ; + // surfvol_ intentionally empty, left to initBlackoilSurfvol +} + +/// 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, + double epsilon) const { + const BlackoilState* that = dynamic_cast (&other); + bool equal = that != 0; + equal = equal && SimulatorState::equals (other, epsilon); + equal = equal && SimulatorState::vectorApproxEqual(this->surfacevol(), + that->surfacevol(), + epsilon); + equal = equal && SimulatorState::vectorApproxEqual(this->gasoilratio(), + that->gasoilratio(), + epsilon); + return equal; +} diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp index a13184d4..5cf8156e 100644 --- a/opm/core/simulator/BlackoilState.hpp +++ b/opm/core/simulator/BlackoilState.hpp @@ -22,37 +22,17 @@ #include #include -#include +#include #include -#include namespace Opm { /// Simulator state for a blackoil simulator. - class BlackoilState + class BlackoilState : public SimulatorState { public: - - void init(const UnstructuredGrid& g, const int num_phases) - { - num_phases_ = num_phases; - press_.resize(g.number_of_cells, 0.0); - fpress_.resize(g.number_of_faces, 0.0); - flux_.resize(g.number_of_faces, 0.0); - sat_.resize(num_phases * g.number_of_cells, 0.0); - for (int cell = 0; cell < g.number_of_cells; ++cell) { - // Defaulting the second saturation to 1.0. - // This will usually be oil in a water-oil case, - // gas in an oil-gas case. - // For proper initialization, one should not rely on this, - // but use available phase information instead. - sat_[num_phases*cell + 1] = 1.0; - } - gor_.resize(g.number_of_cells, 0.0); - } - - enum ExtremalSat { MinSat, MaxSat }; + virtual void init(const UnstructuredGrid& g, 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 @@ -60,80 +40,20 @@ namespace Opm /// are unchanged. void setFirstSat(const std::vector& cells, const Opm::BlackoilPropertiesInterface& props, - ExtremalSat es) - { - if (cells.empty()) { - return; - } - const int n = cells.size(); - assert(n > 0); - std::vector smin(num_phases_*n); - std::vector smax(num_phases_*n); - props.satRange(n, &cells[0], &smin[0], &smax[0]); - 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]; - } - } + ExtremalSat es); - int numPhases() const - { - return num_phases_; - } + virtual bool equals(const SimulatorState& other, + double epsilon = 1e-8) const; - - - bool equals(const BlackoilState& other, double epsilon = 1e-8) const { - bool equal = (num_phases_ == other.num_phases_); - - equal = equal && (vectorApproxEqual( pressure() , other.pressure() , epsilon)); - equal = equal && (vectorApproxEqual( facepressure() , other.facepressure() , epsilon)); - equal = equal && (vectorApproxEqual( faceflux() , other.faceflux() , epsilon)); - equal = equal && (vectorApproxEqual( surfacevol() , other.surfacevol() , epsilon)); - equal = equal && (vectorApproxEqual( saturation() , other.saturation() , epsilon)); - equal = equal && (vectorApproxEqual( gasoilratio() , other.gasoilratio() , epsilon)); - - return equal; - } - - - std::vector& pressure () { return press_ ; } - std::vector& facepressure() { return fpress_; } - std::vector& faceflux () { return flux_ ; } std::vector& surfacevol () { return surfvol_; } - std::vector& saturation () { return sat_ ; } std::vector& gasoilratio () { return gor_ ; } - const std::vector& pressure () const { return press_ ; } - const std::vector& facepressure() const { return fpress_; } - const std::vector& faceflux () const { return flux_ ; } const std::vector& surfacevol () const { return surfvol_; } - const std::vector& saturation () const { return sat_ ; } const std::vector& gasoilratio () const { return gor_ ; } private: - int num_phases_; - std::vector press_ ; - std::vector fpress_; - std::vector flux_ ; std::vector surfvol_; - std::vector sat_ ; std::vector gor_ ; - - - static bool vectorApproxEqual(const std::vector& v1, const std::vector& v2 , double epsilon) { - if (v1.size() != v2.size()) - return false; - - for (size_t i = 0; i < v1.size(); i++) - if (std::abs(v1[i] - v2[i]) > epsilon * (std::abs(v1[i]) + std::abs(v2[i]))) - return false; - - return true; - } - }; } // namespace Opm diff --git a/opm/core/simulator/SimulatorState.cpp b/opm/core/simulator/SimulatorState.cpp new file mode 100644 index 00000000..3681c9d6 --- /dev/null +++ b/opm/core/simulator/SimulatorState.cpp @@ -0,0 +1,96 @@ +#include +#include + +#include + +using namespace Opm; + +bool +SimulatorState::equals (const SimulatorState& other, + double epsilon) const { + bool equal = (num_phases_ == other.num_phases_); + + // if we use &=, then all the tests will be run regardless + equal = equal && vectorApproxEqual( pressure() , other.pressure() , epsilon); + equal = equal && vectorApproxEqual( facepressure() , other.facepressure() , epsilon); + equal = equal && vectorApproxEqual( faceflux() , other.faceflux() , epsilon); + equal = equal && vectorApproxEqual( saturation() , other.saturation() , epsilon); + + return equal; +} + +bool +SimulatorState::vectorApproxEqual(const std::vector& v1, + const std::vector& v2, + double epsilon) { + if (v1.size() != v2.size()) { + return false; + } + + for (size_t i = 0; i < v1.size(); i++) { + const double diff = std::abs(v1[i] - v2[i]); + const double scale = std::abs(v1[i]) + std::abs(v2[i]); + if (diff > epsilon * scale) { + return false; + } + } + + return true; +} + +void +SimulatorState::init(const UnstructuredGrid& g, int num_phases) +{ + num_phases_ = num_phases; + press_.resize(g.number_of_cells, 0.0); + fpress_.resize(g.number_of_faces, 0.0); + flux_.resize(g.number_of_faces, 0.0); + sat_.resize(num_phases_ * g.number_of_cells, 0.0); + for (int cell = 0; cell < g.number_of_cells; ++cell) { + // Defaulting the second saturation to 1.0. + // This will usually be oil in a water-oil case, + // gas in an oil-gas case. + // For proper initialization, one should not rely on this, + // but use available phase information instead. + sat_[num_phases_*cell + 1] = 1.0; + } +} + +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]); + 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]; + } +} + +// 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 new file mode 100644 index 00000000..0baeeb1d --- /dev/null +++ b/opm/core/simulator/SimulatorState.hpp @@ -0,0 +1,78 @@ +// Copyright (C) 2013 Uni Research AS +// This file is licensed under the GNU General Public License v3.0 + +#ifndef OPM_SIMULATORSTATE_HEADER_INCLUDED +#define OPM_SIMULATORSTATE_HEADER_INCLUDED + +#include + +// forward declaration +struct UnstructuredGrid; + +namespace Opm +{ + class SimulatorState + { + public: + + virtual void init(const UnstructuredGrid& g, int num_phases); + + enum ExtremalSat { MinSat, MaxSat }; + + protected: + /** + * 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: + int numPhases() const { return num_phases_; } + + std::vector& pressure () { return press_ ; } + std::vector& facepressure() { return fpress_; } + std::vector& faceflux () { return flux_ ; } + std::vector& saturation () { return sat_ ; } + + const std::vector& pressure () const { return press_ ; } + const std::vector& facepressure() const { return fpress_; } + const std::vector& faceflux () const { return flux_ ; } + const std::vector& saturation () const { return sat_ ; } + + /** + * Compare this state with another, to see if they are different + * only within a small margin. + */ + virtual bool equals(const SimulatorState& other, + double epsilon = 1e-8) const; + private: + int num_phases_; + std::vector press_ ; + std::vector fpress_; + std::vector flux_ ; + std::vector sat_ ; + + protected: + /** + * Check if two vectors are equal within a margin. + * + * @param epsilon Relative difference that is tolerated for the + * vectors to still be considered equal. + * + * @return True if every element is within the margin, false if + * there is at least one that is not. + */ + static bool vectorApproxEqual(const std::vector& v1, + const std::vector& v2, + double epsilon); + }; + +} // namespace Opm + +#endif // OPM_SIMULATORSTATE_HEADER_INCLUDED diff --git a/opm/core/simulator/TwophaseState.hpp b/opm/core/simulator/TwophaseState.hpp index e1eaa174..21e51d5c 100644 --- a/opm/core/simulator/TwophaseState.hpp +++ b/opm/core/simulator/TwophaseState.hpp @@ -20,77 +20,26 @@ #ifndef OPM_TWOPHASESTATE_HEADER_INCLUDED #define OPM_TWOPHASESTATE_HEADER_INCLUDED -#include #include -#include +#include namespace Opm { /// Simulator state for a two-phase simulator. - class TwophaseState + class TwophaseState : public SimulatorState { public: - - void init(const UnstructuredGrid& g, int num_phases) - { - num_phases_ = num_phases; - press_.resize(g.number_of_cells, 0.0); - fpress_.resize(g.number_of_faces, 0.0); - flux_.resize(g.number_of_faces, 0.0); - sat_.resize(num_phases_ * g.number_of_cells, 0.0); - for (int cell = 0; cell < g.number_of_cells; ++cell) { - // Defaulting the second saturation to 1.0. - // This will usually be oil in a water-oil case, - // gas in an oil-gas case. - // For proper initialization, one should not rely on this, - // but use available phase information instead. - sat_[num_phases_*cell + 1] = 1.0; - } - } - - enum ExtremalSat { MinSat, MaxSat }; - void setFirstSat(const std::vector& cells, const Opm::IncompPropertiesInterface& props, - ExtremalSat es) - { - const 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]); - 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]; - } - } + ExtremalSat es); - int numPhases() const - { - return num_phases_; - } - - std::vector& pressure () { return press_ ; } - std::vector& facepressure() { return fpress_; } - std::vector& faceflux () { return flux_ ; } - std::vector& saturation () { return sat_ ; } - - const std::vector& pressure () const { return press_ ; } - const std::vector& facepressure() const { return fpress_; } - const std::vector& faceflux () const { return flux_ ; } - const std::vector& saturation () const { return sat_ ; } - - private: - int num_phases_; - std::vector press_ ; - std::vector fpress_; - std::vector flux_ ; - std::vector sat_ ; + virtual bool equals (const SimulatorState& other, + double epsilon = 1e-8) const; }; } // namespace Opm +#include "TwophaseState_impl.hpp" #endif // OPM_TWOPHASESTATE_HEADER_INCLUDED diff --git a/opm/core/simulator/TwophaseState_impl.hpp b/opm/core/simulator/TwophaseState_impl.hpp new file mode 100644 index 00000000..d54ecae5 --- /dev/null +++ b/opm/core/simulator/TwophaseState_impl.hpp @@ -0,0 +1,22 @@ +#ifndef OPM_TWOPHASESTATE_HEADER_INCLUDED +#error Do not include this file directly! +#endif + +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 { + return dynamic_cast (&other) + ? SimulatorState::equals (other, epsilon) + : false; +} + +} // namespace Opm