From 98c49fd52fc407a9dad595525ff27128616a76f5 Mon Sep 17 00:00:00 2001 From: Robert Kloefkorn Date: Sat, 13 Feb 2016 16:02:48 +0100 Subject: [PATCH] AutoDiffMatrix: use operator += to add matrices --- opm/autodiff/AutoDiffBlock.hpp | 29 ++++++++++++++++++++---- opm/autodiff/BlackoilPropsAdFromDeck.cpp | 14 ++++++------ 2 files changed, 32 insertions(+), 11 deletions(-) diff --git a/opm/autodiff/AutoDiffBlock.hpp b/opm/autodiff/AutoDiffBlock.hpp index d8c33322a..3fee8b74e 100644 --- a/opm/autodiff/AutoDiffBlock.hpp +++ b/opm/autodiff/AutoDiffBlock.hpp @@ -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)); diff --git a/opm/autodiff/BlackoilPropsAdFromDeck.cpp b/opm/autodiff/BlackoilPropsAdFromDeck.cpp index ecc5f2f09..9dfb96fbf 100644 --- a/opm/autodiff/BlackoilPropsAdFromDeck.cpp +++ b/opm/autodiff/BlackoilPropsAdFromDeck.cpp @@ -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);