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

Fixed TransportSolverThophaseAd.cpp

Browse Source
This commit is contained in:
babrodtk 2015-08-27 13:43:02 +02:00
parent 6deb3e2c4a
commit 51b85276ec
3 changed files with 43 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
CMakeLists_files.cmake
View File

@ -35,8 +35,8 @@ list (APPEND MAIN_SOURCE_FILES
opm/autodiff/GridHelpers.cpp
opm/autodiff/ImpesTPFAAD.cpp
opm/autodiff/SimulatorFullyImplicitBlackoilOutput.cpp
# opm/autodiff/SimulatorIncompTwophaseAd.cpp
# opm/autodiff/TransportSolverTwophaseAd.cpp
opm/autodiff/SimulatorIncompTwophaseAd.cpp
opm/autodiff/TransportSolverTwophaseAd.cpp
opm/autodiff/BlackoilPropsAdFromDeck.cpp
# opm/autodiff/SolventPropsAdFromDeck.cpp
opm/autodiff/BlackoilModelParameters.cpp

34
opm/autodiff/AutoDiffMatrix.hpp
View File

@ -265,6 +265,40 @@ namespace Opm
AutoDiffMatrix operator/(const double rhs) const
{
switch (type_) {
case Z:
return *this;
case I:
{
AutoDiffMatrix retval(*this);
retval.type_ = D;
retval.d_.assign(rows_, 1.0/rhs);
return retval;
}
case D:
{
AutoDiffMatrix retval(*this);
for (double& elem : retval.d_) {
elem /= rhs;
}
return retval;
}
case S:
{
AutoDiffMatrix retval(*this);
retval.s_ /= rhs;
return retval;
}
}
}
Eigen::VectorXd operator*(const Eigen::VectorXd& rhs) const
{
assert(cols_ == rhs.size());

12
opm/autodiff/TransportSolverTwophaseAd.cpp
View File

@ -96,6 +96,7 @@ namespace Opm
{
typedef Eigen::Array<double, Eigen::Dynamic, 2, Eigen::RowMajor> TwoCol;
typedef Eigen::Array<double, Eigen::Dynamic, 4, Eigen::RowMajor> FourCol;
typedef Eigen::SparseMatrix<double> S;
typedef typename ADB::V V;
typedef typename ADB::M M;
const int nc = props.numCells();
@ -117,8 +118,8 @@ namespace Opm
// dkro/dsw = col(1) - col(3)
V dkrw = dkr.leftCols<1>() - dkr.rightCols<2>().leftCols<1>();
V dkro = dkr.leftCols<2>().rightCols<1>() - dkr.rightCols<1>();
M krwjac(nc,nc);
M krojac(nc,nc);
S krwjac(nc,nc);
S krojac(nc,nc);
auto sizes = Eigen::ArrayXi::Ones(nc);
krwjac.reserve(sizes);
krojac.reserve(sizes);
@ -127,8 +128,8 @@ namespace Opm
krojac.insert(c,c) = dkro(c);
}
const double* mu = props.viscosity();
std::vector<M> dmw = { krwjac/mu[0] };
std::vector<M> dmo = { krojac/mu[1] };
std::vector<M> dmw = { M(krwjac)/mu[0] };
std::vector<M> dmo = { M(krojac)/mu[1] };
std::vector<ADB> pmobc = { ADB::function(krw / mu[0], std::move(dmw)) ,
ADB::function(kro / mu[1], std::move(dmo)) };
@ -229,7 +230,8 @@ namespace Opm
std::cout << "Residual l2-norm = " << res_norm << std::endl;
// Solve linear system.
Eigen::SparseMatrix<double, Eigen::RowMajor> smatr = transport_residual.derivative()[0];
Eigen::SparseMatrix<double, Eigen::RowMajor> smatr;
transport_residual.derivative()[0].toSparse(smatr);
assert(smatr.isCompressed());
V ds(nc);
LinearSolverInterface::LinearSolverReport rep