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 sim_simple.cpp

Browse Source
This commit is contained in:
babrodtk 2015-08-27 14:03:15 +02:00
parent 95e9ca6d2a
commit 171cbbe3bb
2 changed files with 16 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
CMakeLists_files.cmake
View File

@ -83,11 +83,11 @@ endif()
# originally generated with the command:
# find tutorials examples -name '*.c*' -printf '\t%p\n' | sort
list (APPEND EXAMPLE_SOURCE_FILES
# examples/find_zero.cpp
examples/find_zero.cpp
examples/flow.cpp
# examples/flow_solvent.cpp
# examples/sim_2p_incomp_ad.cpp
# examples/sim_simple.cpp
examples/flow_solvent.cpp
examples/sim_2p_incomp_ad.cpp
examples/sim_simple.cpp
# examples/opm_init_check.cpp
)

20
examples/sim_simple.cpp
View File

@ -66,6 +66,7 @@ phaseMobility(const Opm::IncompPropertiesInterface& props,
{
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();
@ -79,8 +80,8 @@ phaseMobility(const Opm::IncompPropertiesInterface& props,
V kro = kr.rightCols<1>();
V dkrw = dkr.leftCols<1>(); // Left column is top-left of dkr/ds 2x2 matrix.
V dkro = -dkr.rightCols<1>(); // Right column is bottom-right of dkr/ds 2x2 matrix.
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);
@ -89,8 +90,8 @@ phaseMobility(const Opm::IncompPropertiesInterface& props,
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)) };
@ -116,6 +117,7 @@ try
typedef Opm::AutoDiffBlock<double> ADB;
typedef ADB::V V;
typedef ADB::M M;
typedef Eigen::SparseMatrix<double> S;
Opm::time::StopWatch clock;
clock.start();
@ -211,13 +213,14 @@ try
// residual.derived()[0].
#if HAVE_SUITESPARSE_UMFPACK_H
typedef Eigen::UmfPackLU<M> LinSolver;
typedef Eigen::UmfPackLU<S> LinSolver;
#else
typedef Eigen::BiCGSTAB<M> LinSolver;
typedef Eigen::BiCGSTAB<S> LinSolver;
#endif // HAVE_SUITESPARSE_UMFPACK_H
LinSolver solver;
M pmatr = residual.derivative()[0];
S pmatr;
residual.derivative()[0].toSparse(pmatr);
pmatr.coeffRef(0,0) *= 2.0;
pmatr.makeCompressed();
solver.compute(pmatr);
@ -275,7 +278,8 @@ try
std::cout << "res_norm[" << it << "] = "
<< res_norm << std::endl;
M smatr = transport_residual.derivative()[0];
S smatr;
transport_residual.derivative()[0].toSparse(smatr);
smatr.makeCompressed();
solver.compute(smatr);
if (solver.info() != Eigen::Success) {