OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

Added the use of boost::any

Browse Source
This commit is contained in:
babrodtk 2015-08-28 11:34:46 +02:00
parent 624790e7e3
commit df1d0f795d
1 changed files with 67 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

96
opm/autodiff/AutoDiffMatrix.hpp
View File

@ -24,6 +24,7 @@
#include <Eigen/Eigen> #include <Eigen/Eigen>
#include <Eigen/Sparse> #include <Eigen/Sparse>
#include <boost/any.hpp>
#include <opm/core/utility/platform_dependent/reenable_warnings.h> #include <opm/core/utility/platform_dependent/reenable_warnings.h>
@ -38,6 +39,9 @@ namespace Opm
class AutoDiffMatrix class AutoDiffMatrix
{ {
public: public:
typedef std::vector<double> Diag;
typedef Eigen::SparseMatrix<double> Sparse;
AutoDiffMatrix() AutoDiffMatrix()
: type_(Z), : type_(Z),
rows_(0), rows_(0),
@ -72,7 +76,7 @@ namespace Opm
: type_(D), : type_(D),
rows_(d.rows()), rows_(d.rows()),
cols_(d.cols()), cols_(d.cols()),
d_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()) data_(Diag(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()))
{ {
} }
@ -82,7 +86,7 @@ namespace Opm
: type_(S), : type_(S),
rows_(s.rows()), rows_(s.rows()),
cols_(s.cols()), cols_(s.cols()),
s_(s) data_(s)
{ {
} }
@ -114,8 +118,7 @@ namespace Opm
std::swap(type_, other.type_); std::swap(type_, other.type_);
std::swap(rows_, other.rows_); std::swap(rows_, other.rows_);
std::swap(cols_, other.cols_); std::swap(cols_, other.cols_);
d_.swap(other.d_); data_.swap(other.data_);
s_.swap(other.s_);
} }
@ -258,13 +261,14 @@ namespace Opm
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.type_ = D; retval.type_ = D;
retval.d_.assign(rows_, rhs); retval.data_ = boost::any(Diag(rows_, rhs));
return retval; return retval;
} }
case D: case D:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
for (double& elem : retval.d_) { Diag& d = boost::any_cast<Diag&>(retval.data_);
for (double& elem : d) {
elem *= rhs; elem *= rhs;
} }
return retval; return retval;
@ -272,7 +276,8 @@ namespace Opm
case S: case S:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.s_ *= rhs; Sparse& s = boost::any_cast<Sparse&>(retval.data_);
s *= rhs;
return retval; return retval;
} }
default: default:
@ -294,13 +299,14 @@ namespace Opm
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.type_ = D; retval.type_ = D;
retval.d_.assign(rows_, 1.0/rhs); boost::any_cast<Diag&>(retval.data_).assign(rows_, 1.0/rhs);
return retval; return retval;
} }
case D: case D:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
for (double& elem : retval.d_) { Diag& d = boost::any_cast<Diag&>(retval.data_);
for (double& elem : d) {
elem /= rhs; elem /= rhs;
} }
return retval; return retval;
@ -308,7 +314,8 @@ namespace Opm
case S: case S:
{ {
AutoDiffMatrix retval(*this); AutoDiffMatrix retval(*this);
retval.s_ /= rhs; Sparse& s = boost::any_cast<Sparse&>(retval.data_);
s /= rhs;
return retval; return retval;
} }
default: default:
@ -330,9 +337,15 @@ namespace Opm
case I: case I:
return rhs; return rhs;
case D: case D:
return Eigen::Map<const Eigen::VectorXd>(d_.data(), rows_) * rhs; {
const Diag& d = boost::any_cast<const Diag&>(data_);
return Eigen::Map<const Eigen::VectorXd>(d.data(), rows_) * rhs;
}
case S: case S:
return s_ * rhs; {
const Sparse& s = boost::any_cast<const Sparse&>(data_);
return s * rhs;
}
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
} }
@ -351,7 +364,7 @@ namespace Opm
retval.type_ = D; retval.type_ = D;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.d_.assign(lhs.rows_, 2.0); retval.data_ = boost::any(Diag(lhs.rows_, 2.0));
return retval; return retval;
} }
@ -361,8 +374,9 @@ namespace Opm
assert(lhs.type_ == D); assert(lhs.type_ == D);
assert(rhs.type_ == I); assert(rhs.type_ == I);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
Diag& d = boost::any_cast<Diag&>(retval.data_);
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
retval.d_[r] += 1.0; d[r] += 1.0;
} }
return retval; return retval;
} }
@ -372,8 +386,10 @@ namespace Opm
assert(lhs.type_ == D); assert(lhs.type_ == D);
assert(rhs.type_ == D); assert(rhs.type_ == D);
AutoDiffMatrix retval = lhs; AutoDiffMatrix retval = lhs;
Diag& d_lhs = boost::any_cast<Diag&>(retval.data_);
const Diag& d_rhs = boost::any_cast<const Diag&>(rhs.data_);
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
retval.d_[r] += rhs.d_[r]; d_lhs[r] += d_rhs[r];
} }
return retval; return retval;
} }
@ -387,7 +403,8 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.s_ = lhs.s_ + ident; const Sparse& s = boost::any_cast<const Sparse&>(lhs.data_);
retval.data_ = static_cast<Sparse>(s + ident);
return retval; return retval;
} }
@ -396,11 +413,12 @@ namespace Opm
assert(lhs.type_ == S); assert(lhs.type_ == S);
assert(rhs.type_ == D); assert(rhs.type_ == D);
AutoDiffMatrix retval; AutoDiffMatrix retval;
Eigen::SparseMatrix<double> diag = spdiag(rhs.d_); Sparse diag = spdiag(boost::any_cast<const Diag&>(rhs.data_));
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.s_ = lhs.s_ + diag; const Sparse& s = boost::any_cast<const Sparse&>(lhs.data_);
retval.data_ = static_cast<Sparse>(s + diag);
return retval; return retval;
} }
@ -412,7 +430,9 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.s_ = lhs.s_ + rhs.s_; const Sparse& s_lhs = boost::any_cast<const Sparse&>(lhs.data_);
const Sparse& s_rhs = boost::any_cast<const Sparse&>(rhs.data_);
retval.data_ = static_cast<Sparse>(s_lhs + s_rhs);
return retval; return retval;
} }
@ -428,9 +448,12 @@ namespace Opm
retval.type_ = D; retval.type_ = D;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
retval.d_.resize(lhs.rows_); retval.data_ = boost::any(Diag(lhs.rows_));
Diag& d = boost::any_cast<Diag&>(retval.data_);
const Diag& d_lhs = boost::any_cast<const Diag&>(lhs.data_);
const Diag& d_rhs = boost::any_cast<const Diag&>(rhs.data_);
for (int r = 0; r < lhs.rows_; ++r) { for (int r = 0; r < lhs.rows_; ++r) {
retval.d_[r] = lhs.d_[r] * rhs.d_[r]; d[r] = d_lhs[r] * d_rhs[r];
} }
return retval; return retval;
} }
@ -443,7 +466,11 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastDiagSparseProduct(lhs.d_, rhs.s_, retval.s_); retval.data_ = boost::any(Sparse(retval.rows_, retval.cols_)); //FIXME: Superfluous?
const Diag& a = boost::any_cast<const Diag&>(lhs.data_);
const Sparse& b = boost::any_cast<const Sparse&>(rhs.data_);
Sparse& c = boost::any_cast<Sparse&>(retval.data_);
fastDiagSparseProduct(a, b, c);
return retval; return retval;
} }
@ -455,7 +482,11 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastSparseDiagProduct(lhs.s_, rhs.d_, retval.s_); retval.data_ = boost::any(Sparse(retval.rows_, retval.cols_)); //FIXME: Superfluous?
const Sparse& a = boost::any_cast<const Sparse&>(lhs.data_);
const Diag& b = boost::any_cast<const Diag&>(rhs.data_);
Sparse& c = boost::any_cast<Sparse&>(retval.data_);
fastSparseDiagProduct(a, b, c);
return retval; return retval;
} }
@ -467,7 +498,11 @@ namespace Opm
retval.type_ = S; retval.type_ = S;
retval.rows_ = lhs.rows_; retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_; retval.cols_ = rhs.cols_;
fastSparseProduct(lhs.s_, rhs.s_, retval.s_); retval.data_ = boost::any(Sparse(retval.rows_, retval.cols_)); //FIXME: Superfluous?
const Sparse& a = boost::any_cast<const Sparse&>(lhs.data_);
const Sparse& b = boost::any_cast<const Sparse&>(rhs.data_);
Sparse& c = boost::any_cast<Sparse&>(retval.data_);
fastSparseProduct(a, b, c);
return retval; return retval;
} }
@ -483,10 +518,10 @@ namespace Opm
s = spdiag(Eigen::VectorXd::Ones(rows_)); s = spdiag(Eigen::VectorXd::Ones(rows_));
return; return;
case D: case D:
s = spdiag(d_); s = spdiag(boost::any_cast<const Diag&>(data_));
return; return;
case S: case S:
s = s_; s = boost::any_cast<const Sparse&>(data_);
return; return;
} }
} }
@ -512,7 +547,7 @@ namespace Opm
case D: case D:
return rows_; return rows_;
case S: case S:
return s_.nonZeros(); return boost::any_cast<const Sparse&>(data_).nonZeros();
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
} }
@ -527,9 +562,9 @@ namespace Opm
case I: case I:
return (row == col) ? 1.0 : 0.0; return (row == col) ? 1.0 : 0.0;
case D: case D:
return (row == col) ? d_[row] : 0.0; return (row == col) ? boost::any_cast<const Diag&>(data_)[row] : 0.0;
case S: case S:
return s_.coeff(row, col); return boost::any_cast<const Sparse&>(data_).coeff(row, col);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << type_);
} }
@ -540,8 +575,11 @@ namespace Opm
MatrixType type_; MatrixType type_;
int rows_; int rows_;
int cols_; int cols_;
boost::any data_;
/*
std::vector<double> d_; std::vector<double> d_;
Eigen::SparseMatrix<double> s_; Eigen::SparseMatrix<double> s_;
*/
template <class V> template <class V>
static inline static inline