From 8dc4adbe51d8fa83552524c44b85b5256c42ea3c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 23 Aug 2010 09:23:09 +0000 Subject: [PATCH 1/6] Move solvers/common solvers/euler and solvers/mimetic to dune-porsol --- dune/porsol/common/linearInterpolation.hpp | 90 ++++++++++++++++++++++ 1 file changed, 90 insertions(+) create mode 100644 dune/porsol/common/linearInterpolation.hpp diff --git a/dune/porsol/common/linearInterpolation.hpp b/dune/porsol/common/linearInterpolation.hpp new file mode 100644 index 00000000..830545f8 --- /dev/null +++ b/dune/porsol/common/linearInterpolation.hpp @@ -0,0 +1,90 @@ +//=========================================================================== +// +// File: linearInterpolation.hpp +// +// Created: Tue Sep 9 12:49:39 2008 +// +// Author(s): Atgeirr F Rasmussen +// +// $Date$ +// +// $Revision$ +// +//=========================================================================== + +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. + + This file is part of The Open Reservoir Simulator Project (OpenRS). + + OpenRS 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. + + OpenRS 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 OpenRS. If not, see . +*/ + +#ifndef OPENRS_LINEARINTERPOLATION_HEADER +#define OPENRS_LINEARINTERPOLATION_HEADER + + +#include +#include + +namespace Dune +{ + + /** Linear interpolation. + * Given an increasing vector xv of parameter values and + * a vector yv of point values of the same size, + * the function returns ... + */ + template + T linearInterpolation(const std::vector& xv, + const std::vector& yv, + double x) + { + std::vector::const_iterator lb = std::lower_bound(xv.begin(), xv.end(), x); + int lb_ix = lb - xv.begin(); + if (lb_ix == 0) { + return yv[0]; + } else if (lb_ix == int(xv.size())) { + return yv.back(); + } else { + double w = (x - xv[lb_ix - 1])/(xv[lb_ix] - xv[lb_ix - 1]); + return (1.0 - w)*yv[lb_ix - 1] + w*yv[lb_ix]; + } + } + + /// @brief + /// @todo Doc me! + /// @tparam + /// @param + /// @return + template + T linearInterpolationDerivative(const std::vector& xv, + const std::vector& yv, + double x) + { + double epsilon = 1e-4; // @@ Ad hoc, should choose based on xv. + double x_low = std::max(xv[0], x - epsilon); + double x_high = std::min(xv.back(), x + epsilon); + T low = linearInterpolation(xv, yv, x_low); + T high = linearInterpolation(xv, yv, x_high); + return (high - low)/(x_high - x_low); + } + + +} // namespace Dune + + + +#endif // OPENRS_LINEARINTERPOLATION_HEADER From df226cc2a4b0fc5505fc9dd1f898deda990abaa5 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/6] 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 00000000..5b416cbd --- /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 4f4e837a87b23dcbc5f1ab55da076a340e4f16b2 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/6] 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 5b416cbd..8aa42d6b 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 00000000..03c060be --- /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 ebb352dc042e33ac8ee6a85c1896b636879fcb95 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/6] 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 8aa42d6b..644ad762 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 6f9ef5a32cd5f46be5ac65dab5cc39b4257ae544 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/6] 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 644ad762..cf3bef25 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 From 409c2751316cc16a543fc8f3ecb3c4b2d2e407f0 Mon Sep 17 00:00:00 2001 From: convert-repo <> Date: Wed, 21 Dec 2011 11:21:09 +0000 Subject: [PATCH 6/6] update tags --- .hgtags | 1 + 1 file changed, 1 insertion(+) create mode 100644 .hgtags diff --git a/.hgtags b/.hgtags new file mode 100644 index 00000000..df0776bb --- /dev/null +++ b/.hgtags @@ -0,0 +1 @@ +b3cc85c2b5ab3f5b283244624d06a1afab9f8d67 Version-InitialImport-svn935