From 0e306b45ef2e97e9c26f89cf33044acc93288c55 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Mon, 28 Jun 2010 22:47:55 +0000 Subject: [PATCH 01/38] For kicks and giggles, add an (untested) C+BLAS/LAPACK implementation of the 'ip_simple' mimetic inner product. Suggested by: jrn. --- blas_lapack.h | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) create mode 100644 blas_lapack.h diff --git a/blas_lapack.h b/blas_lapack.h new file mode 100644 index 000000000..10603f7a9 --- /dev/null +++ b/blas_lapack.h @@ -0,0 +1,31 @@ +// C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} +void dgemm_(char *transA , char *transB , + const MAT_SIZE_T* m, const MAT_SIZE_T* n , const MAT_SIZE_T* k , + const double* a1, const double* A , const MAT_SIZE_T* ldA, + const double* B, const MAT_SIZE_T* ldB, + const double* a2, double* C , const MAT_SIZE_T* ldC); + +// C <- a1*A*A' + a2*C *or* C <- a1*A'*A + a2*C +void dsyrk_(char *uplo, char *trans, + const MAT_SIZE_T *n , const MAT_SIZE_T *k , + const double *a1 , const double *A , const MAT_SIZE_T *ldA, + const double *a2 , double *C , const MAT_SIZE_T *ldC); + + +// B <- a*op(A)*B *or* B <- a*B*op(A) where op(X) \in {X, X.', X'} +void dtrmm_(char *, char *, char *, char *, + const MAT_SIZE_T *m, const MAT_SIZE_T* n , + const double *a, + const double *A, const MAT_SIZE_T* ldA, + double *B, const MAT_SIZE_T* ldB); + + +void dgeqrf_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , + double *A , const MAT_SIZE_T *ld , + double *tau , double *work, + const MAT_SIZE_T *lwork, MAT_SIZE_T *info); + + +void dorgqr_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , const MAT_SIZE_T *k , + double *A , const MAT_SIZE_T *ld , const double *tau, + double *work, const MAT_SIZE_T *lwork, MAT_SIZE_T *info); From a586f5040f629b51bf14b7964f7e94810d9284e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 29 Jun 2010 15:32:03 +0000 Subject: [PATCH 02/38] Appease the -pedantic -ansi gods. --- blas_lapack.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/blas_lapack.h b/blas_lapack.h index 10603f7a9..3c73d1918 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -1,18 +1,18 @@ -// C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} +/* C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} */ void dgemm_(char *transA , char *transB , const MAT_SIZE_T* m, const MAT_SIZE_T* n , const MAT_SIZE_T* k , const double* a1, const double* A , const MAT_SIZE_T* ldA, const double* B, const MAT_SIZE_T* ldB, const double* a2, double* C , const MAT_SIZE_T* ldC); -// C <- a1*A*A' + a2*C *or* C <- a1*A'*A + a2*C +/* C <- a1*A*A' + a2*C *or* C <- a1*A'*A + a2*C */ void dsyrk_(char *uplo, char *trans, const MAT_SIZE_T *n , const MAT_SIZE_T *k , const double *a1 , const double *A , const MAT_SIZE_T *ldA, const double *a2 , double *C , const MAT_SIZE_T *ldC); -// B <- a*op(A)*B *or* B <- a*B*op(A) where op(X) \in {X, X.', X'} +/* B <- a*op(A)*B *or* B <- a*B*op(A) where op(X) \in {X, X.', X'} */ void dtrmm_(char *, char *, char *, char *, const MAT_SIZE_T *m, const MAT_SIZE_T* n , const double *a, From 4cb92b1e1694ad5f86c23437aad1a91be2a5a7c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Sat, 3 Jul 2010 11:07:56 +0000 Subject: [PATCH 03/38] The LAPACK and BLAS operators do not modify their (char*) parameters. Declare these operators as taking (const char*)s. --- blas_lapack.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/blas_lapack.h b/blas_lapack.h index 3c73d1918..4231e5736 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -1,19 +1,19 @@ /* C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} */ -void dgemm_(char *transA , char *transB , +void dgemm_(const char *transA , const char *transB , const MAT_SIZE_T* m, const MAT_SIZE_T* n , const MAT_SIZE_T* k , const double* a1, const double* A , const MAT_SIZE_T* ldA, const double* B, const MAT_SIZE_T* ldB, const double* a2, double* C , const MAT_SIZE_T* ldC); /* C <- a1*A*A' + a2*C *or* C <- a1*A'*A + a2*C */ -void dsyrk_(char *uplo, char *trans, +void dsyrk_(const char *uplo, const char *trans, const MAT_SIZE_T *n , const MAT_SIZE_T *k , const double *a1 , const double *A , const MAT_SIZE_T *ldA, const double *a2 , double *C , const MAT_SIZE_T *ldC); /* B <- a*op(A)*B *or* B <- a*B*op(A) where op(X) \in {X, X.', X'} */ -void dtrmm_(char *, char *, char *, char *, +void dtrmm_(const char *, const char *, const char *, const char *, const MAT_SIZE_T *m, const MAT_SIZE_T* n , const double *a, const double *A, const MAT_SIZE_T* ldA, From f4e9a3ed2f78ac5f4d66d0ff1c8fd8ca77ee310e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 3 Aug 2010 17:01:33 +0000 Subject: [PATCH 04/38] Initial implementation of hybrid system infrastructure. Actual assembly and system solve not currently implemented. This is untested. --- blas_lapack.h | 33 +++++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/blas_lapack.h b/blas_lapack.h index 4231e5736..8066447b3 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -1,3 +1,6 @@ +#ifndef BLAS_LAPACK_H_INCLUDED +#define BLAS_LAPACK_H_INCLUDED + /* C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} */ void dgemm_(const char *transA , const char *transB , const MAT_SIZE_T* m, const MAT_SIZE_T* n , const MAT_SIZE_T* k , @@ -5,6 +8,7 @@ void dgemm_(const char *transA , const char *transB , const double* B, const MAT_SIZE_T* ldB, const double* a2, double* C , const MAT_SIZE_T* ldC); + /* C <- a1*A*A' + a2*C *or* C <- a1*A'*A + a2*C */ void dsyrk_(const char *uplo, const char *trans, const MAT_SIZE_T *n , const MAT_SIZE_T *k , @@ -12,14 +16,6 @@ void dsyrk_(const char *uplo, const char *trans, const double *a2 , double *C , const MAT_SIZE_T *ldC); -/* B <- a*op(A)*B *or* B <- a*B*op(A) where op(X) \in {X, X.', X'} */ -void dtrmm_(const char *, const char *, const char *, const char *, - const MAT_SIZE_T *m, const MAT_SIZE_T* n , - const double *a, - const double *A, const MAT_SIZE_T* ldA, - double *B, const MAT_SIZE_T* ldB); - - void dgeqrf_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , double *A , const MAT_SIZE_T *ld , double *tau , double *work, @@ -29,3 +25,24 @@ void dgeqrf_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , void dorgqr_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , const MAT_SIZE_T *k , double *A , const MAT_SIZE_T *ld , const double *tau, double *work, const MAT_SIZE_T *lwork, MAT_SIZE_T *info); + + +/* y <- a1*op(A)*x + a2*y */ +void dgemv_(const char *trans, + const MAT_SIZE_T *m , const MAT_SIZE_T *n, + const double *a1 , const double *A, const MAT_SIZE_T *ldA , + const double *x, const MAT_SIZE_T *incX, + const double *a2 , double *y, const MAT_SIZE_T *incY); + + +/* y <- a*x + y */ +void daxpy_(const MAT_SIZE_T *n, const double *a, + const double *x, const MAT_SIZE_T *incx, + double *y, const MAT_SIZE_T *incy); + +/* s <- x' * y */ +double ddot_(const MAT_SIZE_T *n, const double *x, const MAT_SIZE_T *incx, + const double *y, const MAT_SIZE_T *incy); + + +#endif /* BLAS_LAPACK_H_INCLUDED */ From 0bf40b266dc051b7090862a346139dd03af664dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 2 Sep 2010 16:25:29 +0000 Subject: [PATCH 06/38] Add local definition of MAT_SIZE_T, contingent upon preprocessor symbol 'MATLAB_MEX_FILE' that is automatically defined by MATLAB's MEX function. Add declarations for factorisation, lin-sys solution, and matrix inversion for (symmetric) positive definite full matrices in full and packed formats. Will be used in the coarse-system assembly process. --- blas_lapack.h | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/blas_lapack.h b/blas_lapack.h index 8066447b3..d57f0e6a2 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -1,6 +1,16 @@ #ifndef BLAS_LAPACK_H_INCLUDED #define BLAS_LAPACK_H_INCLUDED +#if defined(MATLAB_MEX_FILE) && MATLAB_MEX_FILE +#include +#undef MAT_SIZE_T +#define MAT_SIZE_T mwSignedIndex +#endif + +#ifndef MAT_SIZE_T +#define MAT_SIZE_T int +#endif + /* C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} */ void dgemm_(const char *transA , const char *transB , const MAT_SIZE_T* m, const MAT_SIZE_T* n , const MAT_SIZE_T* k , @@ -27,6 +37,24 @@ void dorgqr_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , const MAT_SIZE_T * double *work, const MAT_SIZE_T *lwork, MAT_SIZE_T *info); +/* A <- chol(A) */ +void dpotrf_(const char *uplo, const MAT_SIZE_T *n, + double *A , const MAT_SIZE_T *lda, + MAT_SIZE_T *info); + +/* B <- (A \ (A' \ B)), when A=DPOTRF(A_orig) */ +void dpotrs_(const char *uplo, const MAT_SIZE_T *n , const MAT_SIZE_T *nrhs, + double *A , const MAT_SIZE_T *lda, + double *B , const MAT_SIZE_T *ldb, MAT_SIZE_T *info); + +/* A <- chol(A), packed format. */ +void dpptrf_(const char *uplo, const MAT_SIZE_T *n, + double *Ap , MAT_SIZE_T *info); + +/* A <- (A \ (A' \ eye(n)) when A=DPPTRF(A_orig) (packed format). */ +void dpptri_(const char *uplo, const MAT_SIZE_T *n, + double *Ap , MAT_SIZE_T *info); + /* y <- a1*op(A)*x + a2*y */ void dgemv_(const char *trans, const MAT_SIZE_T *m , const MAT_SIZE_T *n, From 877ede21eb55edc4ff6b955d6c0c67c9f38b3c9a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 12 Oct 2010 07:25:46 +0000 Subject: [PATCH 07/38] Added copyright block to all source code files. --- blas_lapack.h | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/blas_lapack.h b/blas_lapack.h index d57f0e6a2..5dca57223 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -1,3 +1,22 @@ +/* + 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 BLAS_LAPACK_H_INCLUDED #define BLAS_LAPACK_H_INCLUDED From e0b68dded2ef233f79b6787f381ba82518f67c58 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Tue, 12 Oct 2010 07:44:02 +0000 Subject: [PATCH 08/38] Made all C headers includeable from C++. --- blas_lapack.h | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/blas_lapack.h b/blas_lapack.h index 5dca57223..d12621ec9 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -20,6 +20,10 @@ #ifndef BLAS_LAPACK_H_INCLUDED #define BLAS_LAPACK_H_INCLUDED +#ifdef __cplusplus +extern "C" { +#endif + #if defined(MATLAB_MEX_FILE) && MATLAB_MEX_FILE #include #undef MAT_SIZE_T @@ -92,4 +96,9 @@ double ddot_(const MAT_SIZE_T *n, const double *x, const MAT_SIZE_T *incx, const double *y, const MAT_SIZE_T *incy); + +#ifdef __cplusplus +} +#endif + #endif /* BLAS_LAPACK_H_INCLUDED */ From 933fb4229f454a475fbd49dd4eba96f8e33eaaef Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 13 Oct 2010 18:35:15 +0200 Subject: [PATCH 09/38] Use canonical include guards. Suggested by atgeirr. Template: OPM__HEADER_INCLUDED --- blas_lapack.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/blas_lapack.h b/blas_lapack.h index d12621ec9..ac66f625f 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -17,8 +17,8 @@ along with OPM. If not, see . */ -#ifndef BLAS_LAPACK_H_INCLUDED -#define BLAS_LAPACK_H_INCLUDED +#ifndef OPM_BLAS_LAPACK_HEADER_INCLUDED +#define OPM_BLAS_LAPACK_HEADER_INCLUDED #ifdef __cplusplus extern "C" { @@ -101,4 +101,4 @@ double ddot_(const MAT_SIZE_T *n, const double *x, const MAT_SIZE_T *incx, } #endif -#endif /* BLAS_LAPACK_H_INCLUDED */ +#endif /* OPM_BLAS_LAPACK_HEADER_INCLUDED */ From dc58bb8ffc5c61a5904fca0a75ff4b97ddb12f01 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 28 Oct 2010 10:51:59 +0200 Subject: [PATCH 10/38] Declare DGETRF and DGETRS for compressible support. --- blas_lapack.h | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/blas_lapack.h b/blas_lapack.h index ac66f625f..4c0230b61 100644 --- a/blas_lapack.h +++ b/blas_lapack.h @@ -59,7 +59,18 @@ void dorgqr_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , const MAT_SIZE_T * double *A , const MAT_SIZE_T *ld , const double *tau, double *work, const MAT_SIZE_T *lwork, MAT_SIZE_T *info); +/* A <- LU(A) */ +void dgetrf_(const MAT_SIZE_T *m , const MAT_SIZE_T *n , + double *A , const MAT_SIZE_T *ld, + MAT_SIZE_T *ipiv, MAT_SIZE_T *info); +/* B <- A \ B, when A is LU(A) from dgetrf() */ +void dgetrs_(const char *trans, const MAT_SIZE_T *n, + const MAT_SIZE_T *nrhs , + const double *A , const MAT_SIZE_T *lda, + const MAT_SIZE_T *ipiv , double *B, + const MAT_SIZE_T *ldb , MAT_SIZE_T *info); + /* A <- chol(A) */ void dpotrf_(const char *uplo, const MAT_SIZE_T *n, double *A , const MAT_SIZE_T *lda, From 146402119b77bc4daadf78a8013775fa7278e850 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 29 Oct 2010 15:08:09 +0200 Subject: [PATCH 11/38] Move source files to sub-dir 'src'. --- blas_lapack.h => src/blas_lapack.h | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename blas_lapack.h => src/blas_lapack.h (100%) diff --git a/blas_lapack.h b/src/blas_lapack.h similarity index 100% rename from blas_lapack.h rename to src/blas_lapack.h 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 12/38] 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 13/38] 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 14/38] 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 15/38] 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 From 63656013cc825390fda0943d3f28673f52f5a414 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 18 Oct 2011 20:41:13 +0200 Subject: [PATCH 16/38] Delete trailing whitespace. --- src/blas_lapack.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/blas_lapack.h b/src/blas_lapack.h index 4c0230b61..7051783a8 100644 --- a/src/blas_lapack.h +++ b/src/blas_lapack.h @@ -70,7 +70,7 @@ void dgetrs_(const char *trans, const MAT_SIZE_T *n, const double *A , const MAT_SIZE_T *lda, const MAT_SIZE_T *ipiv , double *B, const MAT_SIZE_T *ldb , MAT_SIZE_T *info); - + /* A <- chol(A) */ void dpotrf_(const char *uplo, const MAT_SIZE_T *n, double *A , const MAT_SIZE_T *lda, From 6a21bc7ecf97035c9dfd2fc84119c1148ceb1dc8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 8 Dec 2011 12:52:57 +0100 Subject: [PATCH 17/38] Move opmpressure/src into core library directory structure. --- {src => opmcore/linalg}/blas_lapack.h | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {src => opmcore/linalg}/blas_lapack.h (100%) diff --git a/src/blas_lapack.h b/opmcore/linalg/blas_lapack.h similarity index 100% rename from src/blas_lapack.h rename to opmcore/linalg/blas_lapack.h From 9682593496ce9259bb55a1f91b4ed5d9c2638ba2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 12 Dec 2011 11:13:54 +0100 Subject: [PATCH 18/38] Moved code from opmcore/ to opm/core/ --- {opmcore => opm/core}/linalg/blas_lapack.h | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename {opmcore => opm/core}/linalg/blas_lapack.h (100%) diff --git a/opmcore/linalg/blas_lapack.h b/opm/core/linalg/blas_lapack.h similarity index 100% rename from opmcore/linalg/blas_lapack.h rename to opm/core/linalg/blas_lapack.h From f475a14f9d8bede973a1338c8bc956333a2626f6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Wed, 21 Dec 2011 13:29:15 +0100 Subject: [PATCH 19/38] Now fluid cpp files compile successfully. --- .../porsol/common => opm/core/utility}/UniformTableLinear.hpp | 4 ++-- .../common => opm/core/utility}/buildUniformMonotoneTable.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) rename {dune/porsol/common => opm/core/utility}/UniformTableLinear.hpp (98%) rename {dune/porsol/common => opm/core/utility}/buildUniformMonotoneTable.hpp (94%) diff --git a/dune/porsol/common/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp similarity index 98% rename from dune/porsol/common/UniformTableLinear.hpp rename to opm/core/utility/UniformTableLinear.hpp index cf3bef256..c88407704 100644 --- a/dune/porsol/common/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -26,8 +26,8 @@ #include #include -#include -#include +#include +#include namespace Dune { namespace utils { diff --git a/dune/porsol/common/buildUniformMonotoneTable.hpp b/opm/core/utility/buildUniformMonotoneTable.hpp similarity index 94% rename from dune/porsol/common/buildUniformMonotoneTable.hpp rename to opm/core/utility/buildUniformMonotoneTable.hpp index 03c060bec..09501f9a3 100644 --- a/dune/porsol/common/buildUniformMonotoneTable.hpp +++ b/opm/core/utility/buildUniformMonotoneTable.hpp @@ -20,8 +20,8 @@ #ifndef OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED #define OPM_BUILDUNIFORMMONOTONETABLE_HEADER_INCLUDED -#include -#include +#include +#include namespace Dune { namespace utils { From 43a542ed028d06112339626fab4d832cc287fbf2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bj=C3=B8rn=20Spjelkavik?= Date: Thu, 19 Jan 2012 13:50:57 +0100 Subject: [PATCH 20/38] Changed namespace Dune -> namespace Opm. --- opm/core/utility/UniformTableLinear.hpp | 4 ++-- opm/core/utility/buildUniformMonotoneTable.hpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index c88407704..6c4870711 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -29,7 +29,7 @@ #include #include -namespace Dune { +namespace Opm { namespace utils { @@ -255,6 +255,6 @@ namespace Dune { } } // namespace utils -} // namespace Dune +} // namespace Opm #endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED diff --git a/opm/core/utility/buildUniformMonotoneTable.hpp b/opm/core/utility/buildUniformMonotoneTable.hpp index 09501f9a3..bf4ab69ed 100644 --- a/opm/core/utility/buildUniformMonotoneTable.hpp +++ b/opm/core/utility/buildUniformMonotoneTable.hpp @@ -23,7 +23,7 @@ #include #include -namespace Dune { +namespace Opm { namespace utils { template @@ -45,7 +45,7 @@ namespace Dune { } } // namespace utils -} // namespace Dune +} // namespace Opm From e54e818639046ae6d801c6cf72166ba0891e95ac Mon Sep 17 00:00:00 2001 From: "Jostein R. Natvig" Date: Wed, 25 Jan 2012 10:44:37 +0100 Subject: [PATCH 21/38] Remove unnecessary include statement. --- opm/core/utility/UniformTableLinear.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 6c4870711..b7e8d751f 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -27,7 +27,6 @@ #include #include -#include namespace Opm { namespace utils { From d2404762fe9937ea0bce8df211badb28670386ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 1 Mar 2012 14:36:10 +0100 Subject: [PATCH 22/38] Added interface and test for lapack tridiagonal solver. --- opm/core/linalg/blas_lapack.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/opm/core/linalg/blas_lapack.h b/opm/core/linalg/blas_lapack.h index 7051783a8..c82e67c8d 100644 --- a/opm/core/linalg/blas_lapack.h +++ b/opm/core/linalg/blas_lapack.h @@ -71,6 +71,16 @@ void dgetrs_(const char *trans, const MAT_SIZE_T *n, const MAT_SIZE_T *ipiv , double *B, const MAT_SIZE_T *ldb , MAT_SIZE_T *info); +/* B <- A \ B, tridiagonal A with bands DL, D, DU */ +void dgtsv_(const MAT_SIZE_T *n , + const MAT_SIZE_T *nrhs , + double *DL , + double *D , + double *DU , + double *B , + const MAT_SIZE_T *ldb , + MAT_SIZE_T *info); + /* A <- chol(A) */ void dpotrf_(const char *uplo, const MAT_SIZE_T *n, double *A , const MAT_SIZE_T *lda, From f9b6a024eeadf55f1074d335a6dd74294816b829 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 1 Mar 2012 15:22:26 +0100 Subject: [PATCH 23/38] Untabify. --- opm/core/linalg/blas_lapack.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/linalg/blas_lapack.h b/opm/core/linalg/blas_lapack.h index c82e67c8d..7da165076 100644 --- a/opm/core/linalg/blas_lapack.h +++ b/opm/core/linalg/blas_lapack.h @@ -73,13 +73,13 @@ void dgetrs_(const char *trans, const MAT_SIZE_T *n, /* B <- A \ B, tridiagonal A with bands DL, D, DU */ void dgtsv_(const MAT_SIZE_T *n , - const MAT_SIZE_T *nrhs , + const MAT_SIZE_T *nrhs , double *DL , double *D , double *DU , double *B , const MAT_SIZE_T *ldb , - MAT_SIZE_T *info); + MAT_SIZE_T *info); /* A <- chol(A) */ void dpotrf_(const char *uplo, const MAT_SIZE_T *n, From b7bd56ab6ab2b350eb70541e51c1256adb7bf40f Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Thu, 15 Mar 2012 16:17:16 +0100 Subject: [PATCH 24/38] Added solver for band matrix and a test example. --- opm/core/linalg/blas_lapack.h | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/opm/core/linalg/blas_lapack.h b/opm/core/linalg/blas_lapack.h index 7da165076..c33ef4e77 100644 --- a/opm/core/linalg/blas_lapack.h +++ b/opm/core/linalg/blas_lapack.h @@ -34,6 +34,7 @@ extern "C" { #define MAT_SIZE_T int #endif + /* C <- a1*op(A)*op(B) + a2*C where op(X) in {X, X.'} */ void dgemm_(const char *transA , const char *transB , const MAT_SIZE_T* m, const MAT_SIZE_T* n , const MAT_SIZE_T* k , @@ -81,6 +82,18 @@ void dgtsv_(const MAT_SIZE_T *n , const MAT_SIZE_T *ldb , MAT_SIZE_T *info); +/* B <- A \ B, band matrix A stored in AB with kl subdiagonals, ku superdiagonals */ +void dgbsv_(const MAT_SIZE_T *n , + const MAT_SIZE_T *kl , + const MAT_SIZE_T *ku , + const MAT_SIZE_T *nrhs , + double *AB , + const MAT_SIZE_T *ldab , + MAT_SIZE_T *ipiv , + double *B , + const MAT_SIZE_T *ldb , + MAT_SIZE_T *info); + /* A <- chol(A) */ void dpotrf_(const char *uplo, const MAT_SIZE_T *n, double *A , const MAT_SIZE_T *lda, From 4ffda4a95e7776bc099abe899381e56dfe49045c Mon Sep 17 00:00:00 2001 From: Xavier Raynaud Date: Fri, 23 Mar 2012 15:44:32 +0100 Subject: [PATCH 25/38] Added general linear lapack solver. Updated test for band matrices. --- opm/core/linalg/blas_lapack.h | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/opm/core/linalg/blas_lapack.h b/opm/core/linalg/blas_lapack.h index c33ef4e77..be2b956c4 100644 --- a/opm/core/linalg/blas_lapack.h +++ b/opm/core/linalg/blas_lapack.h @@ -94,6 +94,16 @@ void dgbsv_(const MAT_SIZE_T *n , const MAT_SIZE_T *ldb , MAT_SIZE_T *info); +/* B <- A \ B, general solver */ +void dgesv_(const MAT_SIZE_T *n, + const MAT_SIZE_T *nrhs , + double *A , + const MAT_SIZE_T *lda , + MAT_SIZE_T *piv , + double *B , + const MAT_SIZE_T *ldb , + MAT_SIZE_T *info); + /* A <- chol(A) */ void dpotrf_(const char *uplo, const MAT_SIZE_T *n, double *A , const MAT_SIZE_T *lda, From ca046f982797a04d4e0e9a5156f564815abed8aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 11:45:52 +0200 Subject: [PATCH 26/38] Moved UniformTableLinear and related func out of subnamespace utils. --- opm/core/utility/UniformTableLinear.hpp | 3 -- .../utility/buildUniformMonotoneTable.hpp | 34 +++++++++---------- 2 files changed, 16 insertions(+), 21 deletions(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index b7e8d751f..058c5f685 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -29,8 +29,6 @@ #include namespace Opm { - namespace utils { - /// @brief This class uses linear interpolation to compute the value /// (and its derivative) of a function f sampled at uniform points. @@ -253,7 +251,6 @@ namespace Opm { return os; } - } // namespace utils } // namespace Opm #endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED diff --git a/opm/core/utility/buildUniformMonotoneTable.hpp b/opm/core/utility/buildUniformMonotoneTable.hpp index bf4ab69ed..3ebad7c99 100644 --- a/opm/core/utility/buildUniformMonotoneTable.hpp +++ b/opm/core/utility/buildUniformMonotoneTable.hpp @@ -24,27 +24,25 @@ #include namespace Opm { - 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); + 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 Opm From 292772b140d181d78c5c070a70cda0212b8deb1a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 19 Apr 2012 14:11:00 +0200 Subject: [PATCH 27/38] Reinject UniformTableLinear into utils subnamespace for backwards compatibility. --- opm/core/utility/UniformTableLinear.hpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 058c5f685..11c2d4ca9 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -251,6 +251,12 @@ namespace Opm { return os; } + namespace utils + { + using Opm::UniformTableLinear; + } + + } // namespace Opm #endif // OPM_UNIFORMTABLELINEAR_HEADER_INCLUDED From 2a21bf4a7595e460788cf988bd812980e7abb827 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= Date: Fri, 31 Aug 2012 17:01:07 +0200 Subject: [PATCH 28/38] Introduced posibility to change number of sample points for pvt. Did change the PVTW calculation so derivatives are exact. Extended the test functions for pvt and relperm --- opm/core/utility/UniformTableLinear.hpp | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 11c2d4ca9..60b3433dc 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -188,9 +188,14 @@ namespace Opm { UniformTableLinear ::derivative(const double xparam) const { - // Implements ClosestValue policy. + // Implements derivative consistent + // with ClosestValue policy for function + double value; + if(xparam > xmax_ || xparam < xmin_){ + value=0.0; + }else{ double x = std::min(xparam, xmax_); - x = std::max(x, xmin_); + x = std::max(x, xmin_); // Lookup is easy since we are uniform in x. double pos = (x - xmin_)/xdelta_; @@ -200,7 +205,9 @@ namespace Opm { // We are at xmax_ --left; } - return (y_values_[left + 1] - y_values_[left])/xdelta_; + value = (y_values_[left + 1] - y_values_[left])/xdelta_; + } + return value; } From 731f984c87da1fe8d27f59eceb26e67f086a84e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Mon, 3 Sep 2012 15:07:03 +0200 Subject: [PATCH 29/38] Formatting fixes. --- opm/core/utility/UniformTableLinear.hpp | 37 ++++++++++++------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 60b3433dc..2a0fc55f7 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -188,26 +188,25 @@ namespace Opm { UniformTableLinear ::derivative(const double xparam) const { - // Implements derivative consistent - // with ClosestValue policy for function - double value; - if(xparam > xmax_ || xparam < xmin_){ - value=0.0; - }else{ - 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; + // Implements derivative consistent + // with ClosestValue policy for function + double value; + if (xparam > xmax_ || xparam < xmin_) { + value = 0.0; + } else { + 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; + } + value = (y_values_[left + 1] - y_values_[left])/xdelta_; } - value = (y_values_[left + 1] - y_values_[left])/xdelta_; - } - return value; + return value; } From a23267a62f4e9b6b8733cff2f60bf66cd7c35cb7 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Thu, 24 Jan 2013 12:39:18 +0100 Subject: [PATCH 30/38] Restructure tests directory to unit tests Every program that relies on manual inspection has been moved to a new (hopefully short-lived) directory called not-unit/; every remaining file has been given the prefix test_ to indicate that this is the executable test to be run. --- tests/test_sparsevector.cpp | 109 ++++++++++++++++++++++++++++++++++++ 1 file changed, 109 insertions(+) create mode 100644 tests/test_sparsevector.cpp diff --git a/tests/test_sparsevector.cpp b/tests/test_sparsevector.cpp new file mode 100644 index 000000000..a6f89be09 --- /dev/null +++ b/tests/test_sparsevector.cpp @@ -0,0 +1,109 @@ +//=========================================================================== +// +// File: sparsevector_test.cpp +// +// Created: Mon Jun 29 21:00:53 2009 +// +// Author(s): Atgeirr F Rasmussen +// BÃ¥rd Skaflestad +// +// $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 . +*/ + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + +#define BOOST_TEST_MODULE SparseVectorTest +#include + +#include + +using namespace Opm; + +BOOST_AUTO_TEST_CASE(construction_and_queries) +{ + const SparseVector sv1; + BOOST_CHECK(sv1.empty()); + BOOST_CHECK_EQUAL(sv1.size(), 0); + BOOST_CHECK_EQUAL(sv1.nonzeroSize(), 0); + + const int size = 100; + const int num_elem = 9; + const int elem[num_elem] = { 9, 8, 7, 6, 5, 4, 3, 2, 1 }; + const int indices[num_elem] = { 1, 2, 3, 5, 8, 13, 21, 34, 55 }; + const SparseVector sv2(size, elem, elem + num_elem, indices, indices + num_elem); + BOOST_CHECK(!sv2.empty()); + BOOST_CHECK_EQUAL(sv2.size(), size); + BOOST_CHECK_EQUAL(sv2.element(0), 0); + BOOST_CHECK_EQUAL(sv2.element(1), 9); + BOOST_CHECK_EQUAL(sv2.element(2), 8); + BOOST_CHECK_EQUAL(sv2.element(3), 7); + BOOST_CHECK_EQUAL(sv2.element(4), 0); + BOOST_CHECK_EQUAL(sv2.element(5), 6); + BOOST_CHECK_EQUAL(sv2.element(55), 1); + BOOST_CHECK_EQUAL(sv2.element(99), 0); + BOOST_CHECK_EQUAL(sv2.nonzeroSize(), num_elem); + for (int i = 0; i < num_elem; ++i) { + BOOST_CHECK_EQUAL(sv2.nonzeroElement(i), elem[i]); + } + const SparseVector sv2_again(size, elem, elem + num_elem, indices, indices + num_elem); + BOOST_CHECK(sv2 == sv2_again); + SparseVector sv2_append(size, elem, elem + num_elem - 1, indices, indices + num_elem - 1); + BOOST_CHECK_EQUAL(sv2_append.nonzeroSize(), num_elem - 1); + sv2_append.addElement(elem[num_elem - 1], indices[num_elem - 1]); + BOOST_CHECK(sv2 == sv2_append); + SparseVector sv2_append2(size); + for (int i = 0; i < num_elem; ++i) { + sv2_append2.addElement(elem[i], indices[i]); + } + BOOST_CHECK(sv2 == sv2_append2); + sv2_append2.clear(); + SparseVector sv_empty; + BOOST_CHECK(sv2_append2 == sv_empty); + + // Tests that only run in debug mode. +#ifndef NDEBUG + // One element too few. + BOOST_CHECK_THROW(const SparseVector sv3(size, elem, elem + num_elem - 1, indices, indices + num_elem), std::exception); + // One element too many. + BOOST_CHECK_THROW(const SparseVector sv4(size, elem, elem + num_elem, indices, indices + num_elem - 1), std::exception); + // Indices out of range. + BOOST_CHECK_THROW(const SparseVector sv5(4, elem, elem + num_elem, indices, indices + num_elem), std::exception); + // Indices not strictly increasing. Cheating by using the elements as indices. + BOOST_CHECK_THROW(const SparseVector sv5(size, elem, elem + num_elem, elem, elem + num_elem), std::exception); + + // Do not ask for out of range indices. + BOOST_CHECK_THROW(sv1.element(0), std::exception); + BOOST_CHECK_THROW(sv2.element(-1), std::exception); + BOOST_CHECK_THROW(sv2.element(sv2.size()), std::exception); + BOOST_CHECK_THROW(sv2.nonzeroElement(sv2.nonzeroSize()), std::exception); +#endif +} + From 4e9e2f0438e2de325a64188b06caddefbfdc9c70 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 22 Mar 2013 21:39:32 +0100 Subject: [PATCH 31/38] Added test for class UniformTableLinear. Moved from opm-porsol. --- tests/test_uniformtablelinear.cpp | 75 +++++++++++++++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 tests/test_uniformtablelinear.cpp diff --git a/tests/test_uniformtablelinear.cpp b/tests/test_uniformtablelinear.cpp new file mode 100644 index 000000000..544c3004e --- /dev/null +++ b/tests/test_uniformtablelinear.cpp @@ -0,0 +1,75 @@ +/* + 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 . +*/ +#include + +#include + +#if defined(HAVE_DYNAMIC_BOOST_TEST) +#define BOOST_TEST_DYN_LINK +#endif +#define NVERBOSE // to suppress our messages when throwing + + +#define BOOST_TEST_MODULE UniformTableLinearTests +#include +#include + + + +BOOST_AUTO_TEST_CASE(table_operations) +{ + // Make a simple table. + double yva[] = { 1.0, -1.0, 3.0, 4.0, 2.0 }; + const int numvals = sizeof(yva)/sizeof(yva[0]); + std::vector yv(yva, yva + numvals); + const double xmin = 1.0; + const double xdelta = 2.5; + const double xmax = xmin + (numvals - 1)*xdelta; + Opm::utils::UniformTableLinear t1(xmin, xmax, yv); + Opm::utils::UniformTableLinear t1_copy1(1.0, 11.0, yv); + Opm::utils::UniformTableLinear t1_copy2(t1); + + // Check equality. + BOOST_CHECK(t1 == t1_copy1); + BOOST_CHECK(t1 == t1_copy2); + + // Check some evaluations. + for (int i = 0; i < numvals; ++i) { + BOOST_CHECK_EQUAL(t1(xmin + i*xdelta), yv[i]); + } + BOOST_CHECK_EQUAL(t1(2.25), 0.0); + BOOST_CHECK_EQUAL(t1(9.75), 3.0); + BOOST_CHECK_EQUAL(t1.derivative(9.75), -2.0/xdelta); + // Until we implement anything but the ClosestValue end policy, we only test that. + BOOST_CHECK_EQUAL(t1(xmin - 1.0), yv[0]); + BOOST_CHECK_EQUAL(t1(xmax + 1.0), yv.back()); + + // Domains. + BOOST_CHECK_EQUAL(t1.domain().first, xmin); + BOOST_CHECK_EQUAL(t1.domain().second, xmin + (numvals-1)*xdelta); + std::pair new_domain(-100.0, 20.0); + t1.rescaleDomain(new_domain); + BOOST_CHECK_EQUAL(t1.domain().first, new_domain.first); + BOOST_CHECK_EQUAL(t1.domain().second, new_domain.second); + for (int i = 0; i < numvals; ++i) { + BOOST_CHECK_EQUAL(t1(-100.0 + i*120.0/(double(numvals - 1))), yv[i]); + } + BOOST_CHECK_EQUAL(t1(-85.0), 0.0); + BOOST_CHECK(std::fabs(t1.derivative(0.0) + 2.0/30.0) < 1e-14); +} From b6db8534d8693b0e71ca8ca8fa3dd703473a2357 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Sat, 23 Mar 2013 23:46:05 +0100 Subject: [PATCH 32/38] Whitespace cleanup. --- tests/test_uniformtablelinear.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_uniformtablelinear.cpp b/tests/test_uniformtablelinear.cpp index 544c3004e..5d40c42a9 100644 --- a/tests/test_uniformtablelinear.cpp +++ b/tests/test_uniformtablelinear.cpp @@ -51,7 +51,7 @@ BOOST_AUTO_TEST_CASE(table_operations) // Check some evaluations. for (int i = 0; i < numvals; ++i) { - BOOST_CHECK_EQUAL(t1(xmin + i*xdelta), yv[i]); + BOOST_CHECK_EQUAL(t1(xmin + i*xdelta), yv[i]); } BOOST_CHECK_EQUAL(t1(2.25), 0.0); BOOST_CHECK_EQUAL(t1(9.75), 3.0); @@ -68,7 +68,7 @@ BOOST_AUTO_TEST_CASE(table_operations) BOOST_CHECK_EQUAL(t1.domain().first, new_domain.first); BOOST_CHECK_EQUAL(t1.domain().second, new_domain.second); for (int i = 0; i < numvals; ++i) { - BOOST_CHECK_EQUAL(t1(-100.0 + i*120.0/(double(numvals - 1))), yv[i]); + BOOST_CHECK_EQUAL(t1(-100.0 + i*120.0/(double(numvals - 1))), yv[i]); } BOOST_CHECK_EQUAL(t1(-85.0), 0.0); BOOST_CHECK(std::fabs(t1.derivative(0.0) + 2.0/30.0) < 1e-14); From ba2df2d25ba79759eaca263de80891c22ac3b935 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 20 Jun 2013 23:14:39 +0200 Subject: [PATCH 33/38] Disable build kluge that is no longer pertinent The header was introduced (commit 82369f9) as a work-around for a particular interaction in the Autotools-based setup of OPM-Core and the Dune core modules. Notably, Dune's "Enable" trick for Boost failed on some older Autoconf systems. Now that we're using CMake, however, that kluge is no longer needed because we (OPM-Core) always #define HAVE_BOOST 1 i.e., as an explict true/false value. Therefore, we need no longer include . The header will be removed at a later time. --- tests/test_uniformtablelinear.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_uniformtablelinear.cpp b/tests/test_uniformtablelinear.cpp index 5d40c42a9..ed7d0ea40 100644 --- a/tests/test_uniformtablelinear.cpp +++ b/tests/test_uniformtablelinear.cpp @@ -18,8 +18,6 @@ */ #include -#include - #if defined(HAVE_DYNAMIC_BOOST_TEST) #define BOOST_TEST_DYN_LINK #endif From e542fd6104f2db516542f24bdc75b21163c95e4a Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Thu, 1 Aug 2013 09:54:30 +0200 Subject: [PATCH 34/38] Graduate these tests to unit tests Although they don't use Boost::UnitTest, they can at least pass, so we can use them to detect simple compilation and runtime errors, although we miss the semantic check. (If you have time, please make them proper unit tests) --- tests/test_cubic.cpp | 55 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 tests/test_cubic.cpp diff --git a/tests/test_cubic.cpp b/tests/test_cubic.cpp new file mode 100644 index 000000000..278e586cc --- /dev/null +++ b/tests/test_cubic.cpp @@ -0,0 +1,55 @@ +//=========================================================================== +// +// File: monotcubicinterpolator_test.cpp +// +// Created: Tue Dec 8 12:25:30 2009 +// +// Author(s): Atgeirr F Rasmussen +// Bård Skaflestad +// +// $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 . +*/ + +#include "config.h" +#include +using namespace Opm; + +int main() +{ + const int num_v = 3; + double xv[num_v] = {0.0, 1.0, 2.0}; + double fv[num_v] = {10.0, 21.0, 2.0}; + std::vector x(xv, xv + num_v); + std::vector f(fv, fv + num_v); + MonotCubicInterpolator interp(x, f); + interp.evaluate(-1.0); + interp.evaluate(0.0); + interp.evaluate(0.0001); + interp.evaluate(0.5); + interp.evaluate(1.0); + interp.evaluate(2.0); + interp.evaluate(4.0); +} From e38548ebf60e489c25d995a891d4bb51dc654388 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Thu, 1 Aug 2013 10:19:44 +0200 Subject: [PATCH 35/38] Convert cubic interpolator to use Boost::UnitTest --- tests/test_cubic.cpp | 34 ++++++++++++++++++++++++++-------- 1 file changed, 26 insertions(+), 8 deletions(-) diff --git a/tests/test_cubic.cpp b/tests/test_cubic.cpp index 278e586cc..ecbb3f686 100644 --- a/tests/test_cubic.cpp +++ b/tests/test_cubic.cpp @@ -16,6 +16,7 @@ /* Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. Copyright 2009, 2010 Statoil ASA. + Portions Copyright 2013 Uni Research AS. This file is part of The Open Reservoir Simulator Project (OpenRS). @@ -34,10 +35,25 @@ */ #include "config.h" + +/* --- Boost.Test boilerplate --- */ +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE CubicTest +#include +#include + +/* --- our own headers --- */ #include using namespace Opm; -int main() +BOOST_AUTO_TEST_SUITE () + +BOOST_AUTO_TEST_CASE (cubic) { const int num_v = 3; double xv[num_v] = {0.0, 1.0, 2.0}; @@ -45,11 +61,13 @@ int main() std::vector x(xv, xv + num_v); std::vector f(fv, fv + num_v); MonotCubicInterpolator interp(x, f); - interp.evaluate(-1.0); - interp.evaluate(0.0); - interp.evaluate(0.0001); - interp.evaluate(0.5); - interp.evaluate(1.0); - interp.evaluate(2.0); - interp.evaluate(4.0); + BOOST_REQUIRE_CLOSE (interp.evaluate(-1.0), 10., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(0.0), 10., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(0.0001), 10.0011, 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(0.5), 17.375, 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(1.0), 21., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(2.0), 2., 0.00001); + BOOST_REQUIRE_CLOSE (interp.evaluate(4.0), 2., 0.00001); } + +BOOST_AUTO_TEST_SUITE_END() From b48c9df650e1ce216742ccd9574347df39b5f914 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 28 Aug 2013 13:59:03 +0200 Subject: [PATCH 36/38] convert THROW to OPM_THROW --- opm/core/utility/UniformTableLinear.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 2a0fc55f7..16dce38bf 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -228,7 +228,7 @@ namespace Opm { ::setLeftPolicy(RangePolicy rp) { if (rp != ClosestValue) { - THROW("Only ClosestValue RangePolicy implemented."); + OPM_THROW(std::runtime_error, "Only ClosestValue RangePolicy implemented."); } left_ = rp; } @@ -239,7 +239,7 @@ namespace Opm { ::setRightPolicy(RangePolicy rp) { if (rp != ClosestValue) { - THROW("Only ClosestValue RangePolicy implemented."); + OPM_THROW(std::runtime_error, "Only ClosestValue RangePolicy implemented."); } right_ = rp; } From 1f7458e13616d05b5928e65a98f86234e89c7098 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 28 Aug 2013 14:00:35 +0200 Subject: [PATCH 37/38] convert users of the ASSERT and the ASSERT2 macros to standard assert() --- opm/core/utility/UniformTableLinear.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 16dce38bf..2ede43b3d 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -123,8 +123,8 @@ namespace Opm { : xmin_(xmin), xmax_(xmax), y_values_(y_values), left_(ClosestValue), right_(ClosestValue) { - ASSERT(xmax > xmin); - ASSERT(y_values.size() > 1); + assert(xmax > xmin); + assert(y_values.size() > 1); xdelta_ = (xmax - xmin)/(y_values.size() - 1); } @@ -139,8 +139,8 @@ namespace Opm { y_values_(y_values, y_values + num_y_values), left_(ClosestValue), right_(ClosestValue) { - ASSERT(xmax > xmin); - ASSERT(y_values_.size() > 1); + assert(xmax > xmin); + assert(y_values_.size() > 1); xdelta_ = (xmax - xmin)/(y_values_.size() - 1); } From 8db816739546f3a7ff8912bbb9ff032969c17ff5 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Thu, 8 Oct 2015 11:42:15 +0200 Subject: [PATCH 38/38] use the error macros from opm-common --- opm/core/utility/UniformTableLinear.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/UniformTableLinear.hpp b/opm/core/utility/UniformTableLinear.hpp index 2ede43b3d..dc7f75e2a 100644 --- a/opm/core/utility/UniformTableLinear.hpp +++ b/opm/core/utility/UniformTableLinear.hpp @@ -26,7 +26,7 @@ #include #include -#include +#include namespace Opm {