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

Fixed bug in diagonal*vector product

Browse Source
This commit is contained in:
babrodtk 2015-09-01 14:10:18 +02:00
parent 9855d7340f
commit 4c82e9abc7
1 changed files with 14 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
opm/autodiff/AutoDiffMatrix.hpp
View File

@ -84,6 +84,7 @@ namespace Opm
diag_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()), diag_(d.diagonal().array().data(), d.diagonal().array().data() + d.rows()),
sparse_() sparse_()
{ {
assert(rows_ == cols_);
} }
@ -145,9 +146,9 @@ namespace Opm
case I: case I:
return addII(*this, rhs); return addII(*this, rhs);
case D: case D:
return rhs + (*this); return addDI(rhs, *this);
case S: case S:
return rhs + (*this); return addSI(rhs, *this);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
} }
@ -160,7 +161,7 @@ namespace Opm
case D: case D:
return addDD(*this, rhs); return addDD(*this, rhs);
case S: case S:
return rhs + (*this); return addSD(rhs, *this);
default: default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_); OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
} }
@ -189,18 +190,7 @@ namespace Opm
case Z: case Z:
return AutoDiffMatrix(rows_, rhs.cols_); return AutoDiffMatrix(rows_, rhs.cols_);
case I: case I:
switch (rhs.type_) { return rhs;
case Z:
return rhs;
case I:
return rhs;
case D:
return rhs;
case S:
return rhs;
default:
OPM_THROW(std::logic_error, "Invalid AutoDiffMatrix type encountered: " << rhs.type_);
}
case D: case D:
switch (rhs.type_) { switch (rhs.type_) {
case Z: case Z:
@ -251,7 +241,7 @@ namespace Opm
AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs) AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs)
{ {
*this = *this + rhs * -1.0; *this = *this + (rhs * -1.0);
return *this; return *this;
} }
@ -341,7 +331,10 @@ namespace Opm
case I: case I:
return rhs; return rhs;
case D: case D:
return Eigen::Map<const Eigen::VectorXd>(diag_.data(), rows_) * rhs; {
const Eigen::VectorXd d = Eigen::Map<const Eigen::VectorXd>(diag_.data(), rows_);
return d.cwiseProduct(rhs);
}
case S: case S:
return sparse_ * rhs; return sparse_ * rhs;
default: default:
@ -393,11 +386,8 @@ namespace Opm
{ {
assert(lhs.type_ == S); assert(lhs.type_ == S);
assert(rhs.type_ == I); assert(rhs.type_ == I);
AutoDiffMatrix retval; AutoDiffMatrix retval = lhs;
retval.type_ = S; retval.sparse_ += spdiag(Eigen::VectorXd::Ones(lhs.rows_));;
retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_;
retval.sparse_ = lhs.sparse_ + spdiag(Eigen::VectorXd::Ones(lhs.rows_));;
return retval; return retval;
} }
@ -405,11 +395,8 @@ namespace Opm
{ {
assert(lhs.type_ == S); assert(lhs.type_ == S);
assert(rhs.type_ == D); assert(rhs.type_ == D);
AutoDiffMatrix retval; AutoDiffMatrix retval = lhs;
retval.type_ = S; retval.sparse_ += spdiag(rhs.diag_);
retval.rows_ = lhs.rows_;
retval.cols_ = rhs.cols_;
retval.sparse_ = lhs.sparse_ + spdiag(rhs.diag_);
return retval; return retval;
} }