From 09e234c68ccb281d381886ebde8bdc98d7f04ed0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 8 Nov 2010 14:12:10 +0100 Subject: [PATCH 2/5] Created a new utility class, UniformTableLinear. --- dune/porsol/common/UniformTableLinear.hpp | 222 ++++++++++++++++++++++ 1 file changed, 222 insertions(+) create mode 100644 dune/porsol/common/UniformTableLinear.hpp diff --git a/dune/porsol/common/UniformTableLinear.hpp b/dune/porsol/common/UniformTableLinear.hpp new file mode 100644 index 000000000..5b416cbdc --- /dev/null +++ b/dune/porsol/common/UniformTableLinear.hpp @@ -0,0 +1,222 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + 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_UNIFORMTABLELINEAR_HEADER_INCLUDED +#define OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED + +#include +#include +#include +#include + +#include +#include + +namespace Dune { + namespace utils { + + + /// @brief Exception used for domain errors. + struct OutsideDomainException : public std::exception {}; + + + /// @brief This class uses linear interpolation to compute the value + /// (and its derivative) of a function f sampled at uniform points. + /// @tparam T the range type of the function (should be an algebraic ring type) + template + class UniformTableLinear + { + public: + /// @brief Default constructor. + UniformTableLinear(); + + /// @brief Useful constructor. + /// @param xmin the x value corresponding to the first y value. + /// @param xmax the x value corresponding to the last y value. + /// @param y_values vector of range values. + UniformTableLinear(double xmin, + double xmax, + const std::vector& y_values); + + /// @brief Get the domain. + /// @return the domain as a pair of doubles. + std::pair domain(); + + /// @brief Rescale the domain. + /// @param new_domain the new domain as a pair of doubles. + void rescaleDomain(std::pair new_domain); + + /// @brief Evaluate the value at x. + /// @param x a domain value + /// @return f(x) + double operator()(const double x) const; + + /// @brief Evaluate the derivative at x. + /// @param x a domain value + /// @return f'(x) + double derivative(const double x) const; + + /// @brief Equality operator. + /// @param other another UniformTableLinear. + /// @return true if they are represented exactly alike. + bool operator==(const UniformTableLinear& other) const; + + /// @brief Policies for how to behave when trying to evaluate outside the domain. + enum RangePolicy {Throw = 0, ClosestValue = 1, Extrapolate = 2}; + + /// @brief Sets the behavioural policy for evaluation to the left of the domain. + /// @param rp the policy + void setLeftPolicy(RangePolicy rp); + + /// @brief Sets the behavioural policy for evaluation to the right of the domain. + /// @param rp the policy + void setRightPolicy(RangePolicy rp); + + protected: + double xmin_; + double xmax_; + double xdelta_; + std::vector y_values_; + RangePolicy left_; + RangePolicy right_; + }; + + + // Member implementations. + + template + inline + UniformTableLinear + ::UniformTableLinear() + : left_(ClosestValue), right_(ClosestValue) + { + } + + template + inline + UniformTableLinear + ::UniformTableLinear(double xmin, + double xmax, + const std::vector& y_values) + : xmin_(xmin), xmax_(xmax), y_values_(y_values), + left_(ClosestValue), right_(ClosestValue) + { + ASSERT(xmax > xmin); + ASSERT(y_values.size() > 1); + xdelta_ = (xmax - xmin)/(y_values.size() - 1); + } + + template + inline std::pair + UniformTableLinear + ::domain() + { + return std::make_pair(xmin_, xmax_); + } + + template + inline void + UniformTableLinear + ::rescaleDomain(std::pair new_domain) + { + xmin_ = new_domain.first; + xmax_ = new_domain.second; + xdelta_ = (xmax_ - xmin_)/(y_values_.size() - 1); + } + + template + inline double + UniformTableLinear + ::operator()(const double xparam) const + { + // Implements ClosestValue policy. + double x = std::min(xparam, xmax_); + x = std::max(x, xmin_); + + // Lookup is easy since we are uniform in x. + double pos = (x - xmin_)/xdelta_; + double posi = std::floor(pos); + int left = int(posi); + if (left == int(y_values_.size()) - 1) { + // We are at xmax_ + return y_values_.back(); + } + double w = pos - posi; + return (1.0 - w)*y_values_[left] + w*y_values_[left + 1]; + } + + template + inline double + UniformTableLinear + ::derivative(const double xparam) const + { + // Implements ClosestValue policy. + double x = std::min(xparam, xmax_); + x = std::max(x, xmin_); + + // Lookup is easy since we are uniform in x. + double pos = (x - xmin_)/xdelta_; + double posi = std::floor(pos); + int left = int(posi); + if (left == int(y_values_.size()) - 1) { + // We are at xmax_ + --left; + } + return (y_values_[left + 1] - y_values_[left])/xdelta_; + } + + + template + inline bool + UniformTableLinear + ::operator==(const UniformTableLinear& other) const + { + return xmin_ == other.xmin_ + && xdelta_ == other.xdelta_ + && y_values_ == other.y_values_ + && left_ == other.left_ + && right_ == other.right_; + } + + template + inline void + UniformTableLinear + ::setLeftPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + left_ = rp; + } + + template + inline void + UniformTableLinear + ::setRightPolicy(RangePolicy rp) + { + if (rp != ClosestValue) { + THROW("Only ClosestValue RangePolicy implemented."); + } + right_ = rp; + } + + } // namespace utils +} // namespace Dune + +#endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED From b8b9581ae44a6aa8966784501d81c1202b08b956 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 10 Nov 2010 13:31:32 +0100 Subject: [PATCH 3/5] Implemented FluidMatrixInteractionBlackoil init(), kr() and a test prog. --- dune/porsol/common/UniformTableLinear.hpp | 28 +++++++++- .../common/buildUniformMonotoneTable.hpp | 52 +++++++++++++++++++ 2 files changed, 79 insertions(+), 1 deletion(-) create mode 100644 dune/porsol/common/buildUniformMonotoneTable.hpp diff --git a/dune/porsol/common/UniformTableLinear.hpp b/dune/porsol/common/UniformTableLinear.hpp index 5b416cbdc..8aa42d6bc 100644 --- a/dune/porsol/common/UniformTableLinear.hpp +++ b/dune/porsol/common/UniformTableLinear.hpp @@ -46,7 +46,7 @@ namespace Dune { /// @brief Default constructor. UniformTableLinear(); - /// @brief Useful constructor. + /// @brief Construct from vector of y-values. /// @param xmin the x value corresponding to the first y value. /// @param xmax the x value corresponding to the last y value. /// @param y_values vector of range values. @@ -54,6 +54,16 @@ namespace Dune { double xmax, const std::vector& y_values); + /// @brief Construct from array of y-values. + /// @param xmin the x value corresponding to the first y value. + /// @param xmax the x value corresponding to the last y value. + /// @param y_values array of range values. + /// @param num_y_values the number of values in y_values. + UniformTableLinear(double xmin, + double xmax, + const T* y_values, + int num_y_values); + /// @brief Get the domain. /// @return the domain as a pair of doubles. std::pair domain(); @@ -122,6 +132,22 @@ namespace Dune { xdelta_ = (xmax - xmin)/(y_values.size() - 1); } + template + inline + UniformTableLinear + ::UniformTableLinear(double xmin, + double xmax, + const T* y_values, + int num_y_values) + : xmin_(xmin), xmax_(xmax), + y_values_(y_values, y_values + num_y_values), + left_(ClosestValue), right_(ClosestValue) + { + ASSERT(xmax > xmin); + ASSERT(y_values_.size() > 1); + xdelta_ = (xmax - xmin)/(y_values_.size() - 1); + } + template inline std::pair UniformTableLinear diff --git a/dune/porsol/common/buildUniformMonotoneTable.hpp b/dune/porsol/common/buildUniformMonotoneTable.hpp new file mode 100644 index 000000000..03c060bec --- /dev/null +++ b/dune/porsol/common/buildUniformMonotoneTable.hpp @@ -0,0 +1,52 @@ +/* + Copyright 2010 SINTEF ICT, Applied Mathematics. + + 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_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED +#define OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED + +#include +#include + +namespace Dune { + namespace utils { + + template + void buildUniformMonotoneTable(const std::vector& xv, + const std::vector& yv, + const int samples, + UniformTableLinear& table) + { + MonotCubicInterpolator interp(xv, yv); + std::vector uniform_yv(samples); + double xmin = xv[0]; + double xmax = xv.back(); + for (int i = 0; i < samples; ++i) { + double w = double(i)/double(samples - 1); + double x = (1.0 - w)*xmin + w*xmax; + uniform_yv[i] = interp(x); + } + table = UniformTableLinear(xmin, xmax, uniform_yv); + } + + } // namespace utils +} // namespace Dune + + + +#endif // OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED From 5b7052eb114201b15548298af8230ff959759557 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 22 Nov 2010 15:00:26 +0100 Subject: [PATCH 4/5] A large number of additions to start testing compressible tpfa-solver. --- dune/porsol/common/UniformTableLinear.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/dune/porsol/common/UniformTableLinear.hpp b/dune/porsol/common/UniformTableLinear.hpp index 8aa42d6bc..644ad7629 100644 --- a/dune/porsol/common/UniformTableLinear.hpp +++ b/dune/porsol/common/UniformTableLinear.hpp @@ -32,10 +32,6 @@ namespace Dune { namespace utils { - /// @brief Exception used for domain errors. - struct OutsideDomainException : public std::exception {}; - - /// @brief This class uses linear interpolation to compute the value /// (and its derivative) of a function f sampled at uniform points. /// @tparam T the range type of the function (should be an algebraic ring type) From afc31d3b085b0feebc18a3ad60bf78f1e2b6352d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 1 Feb 2011 12:40:05 +0100 Subject: [PATCH 5/5] Added output operator for easy dumping of tables. --- dune/porsol/common/UniformTableLinear.hpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/dune/porsol/common/UniformTableLinear.hpp b/dune/porsol/common/UniformTableLinear.hpp index 644ad7629..cf3bef256 100644 --- a/dune/porsol/common/UniformTableLinear.hpp +++ b/dune/porsol/common/UniformTableLinear.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -101,6 +102,8 @@ namespace Dune { std::vector y_values_; RangePolicy left_; RangePolicy right_; + template + friend std::ostream& operator<<(std::ostream& os, const UniformTableLinear& t); }; @@ -238,6 +241,19 @@ namespace Dune { right_ = rp; } + + template + inline std::ostream& operator<<(std::ostream& os, const UniformTableLinear& t) + { + int n = t.y_values_.size(); + for (int i = 0; i < n; ++i) { + double f = double(i)/double(n - 1); + os << (1.0 - f)*t.xmin_ + f*t.xmax_ + << " " << t.y_values_[i] << '\n'; + } + return os; + } + } // namespace utils } // namespace Dune