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

AutoDiffMatrix: use operator += to add matrices

Browse Source
This commit is contained in:
Robert Kloefkorn 2016-02-13 16:02:48 +01:00
parent fadf5528e3
commit 98c49fd52f
2 changed files with 32 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
opm/autodiff/AutoDiffBlock.hpp
View File

@ -342,7 +342,9 @@ namespace Opm
jac[block] = D2*jac_[block];
}
else {
jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
//jac[block] = D2*jac_[block] + D1*rhs.jac_[block];
jac[block] = D2*jac_[block];
jac[block] += D1*rhs.jac_[block];
}
}
return function(val_ * rhs.val_, std::move(jac));
@ -366,6 +368,8 @@ namespace Opm
M D1(val_.matrix().asDiagonal());
M D2(rhs.val_.matrix().asDiagonal());
M D3((1.0/(rhs.val_*rhs.val_)).matrix().asDiagonal());
M D4;
M D5;
#pragma omp parallel for schedule(dynamic)
for (int block = 0; block < num_blocks; ++block) {
assert(jac_[block].rows() == rhs.jac_[block].rows());
@ -374,14 +378,31 @@ namespace Opm
jac[block] = M( D3.rows(), jac_[block].cols() );
}
else if( jac_[block].nonZeros() == 0 ) {
jac[block] = D3 * ( D1*rhs.jac_[block]) * (-1.0);
#pragma omp critical
if( D4.rows() == 0 ) {
D4 = D3 * D1 * (-1.0);
}
jac[block] = D4*rhs.jac_[block];
}
else if ( rhs.jac_[block].nonZeros() == 0 )
{
jac[block] = D3 * (D2*jac_[block]);
#pragma omp critical
if( D5.rows() == 0 ) {
D5 = D3 * D2 ;
}
jac[block] = D5*jac_[block];
}
else {
jac[block] = D3 * (D2*jac_[block] + (D1*rhs.jac_[block]*(-1.0)));
#pragma omp critical
if( D4.rows() == 0 ) {
D4 = D3 * D1 * (-1.0);
}
#pragma omp critical
if( D5.rows() == 0 ) {
D5 = D3 * D2 ;
}
jac[block] = D5*jac_[block];
jac[block] += D4*rhs.jac_[block];
}
}
return function(val_ / rhs.val_, std::move(jac));

14
opm/autodiff/BlackoilPropsAdFromDeck.cpp
View File

@ -98,7 +98,7 @@ namespace Opm
}
}
materialLawManager->initFromDeck(deck, eclState, compressedToCartesianIdx);
init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims,
init(deck, eclState, materialLawManager, grid.number_of_cells, grid.global_cell, grid.cartdims,
init_rock);
}
@ -399,7 +399,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dmudp_diag, po.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dmudr_diag, rs.derivative()[block], temp);
jacs[block] = jacs[block] + temp;
jacs[block] += temp;
}
return ADB::function(std::move(mu), std::move(jacs));
}
@ -463,7 +463,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dmudp_diag, pg.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dmudr_diag, rv.derivative()[block], temp);
jacs[block] = jacs[block] + temp;
jacs[block] += temp;
}
return ADB::function(std::move(mu), std::move(jacs));
}
@ -578,7 +578,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dbdp_diag, po.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dbdr_diag, rs.derivative()[block], temp);
jacs[block] = jacs[block] + temp;
jacs[block] += temp;
}
return ADB::function(std::move(b), std::move(jacs));
}
@ -643,7 +643,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
fastSparseProduct(dbdp_diag, pg.derivative()[block], jacs[block]);
ADB::M temp;
fastSparseProduct(dbdr_diag, rv.derivative()[block], temp);
jacs[block] = jacs[block] + temp;
jacs[block] += temp;
}
return ADB::function(std::move(b), std::move(jacs));
}
@ -824,7 +824,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < num_blocks; ++block) {
ADB::M temp;
fastSparseProduct(dkr1_ds2_diag, s[phase2]->derivative()[block], temp);
jacs[block] = jacs[block] + temp;
jacs[block] += temp;
}
}
ADB::V val = kr.col(phase1_pos);
@ -885,7 +885,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
for (int block = 0; block < nBlocks; ++block) {
ADB::M temp;
fastSparseProduct(dpc1_ds2_diag, s[phase2]->derivative()[block], temp);
jacs[block] = jacs[block] + temp;
jacs[block] += temp;
}
}
ADB::V val = pc.col(phase1_pos);