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

Merge pull request #582 from dr-robertk/PR/fast-sparse-add

Browse Source 
Fast sparse add.
This commit is contained in:
Atgeirr Flø Rasmussen 2016-02-17 19:55:15 +01:00
commit 8dfc8d82fc
7 changed files with 127 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists_files.cmake
View File

@ -138,7 +138,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/autodiff/BlackoilSolventState.hpp
opm/autodiff/BlackoilMultiSegmentModel.hpp
opm/autodiff/BlackoilMultiSegmentModel_impl.hpp
opm/autodiff/fastSparseProduct.hpp
opm/autodiff/fastSparseOperations.hpp
opm/autodiff/DuneMatrix.hpp
opm/autodiff/ExtractParallelGridInformationToISTL.hpp
opm/autodiff/FlowMain.hpp

5
opm/autodiff/AutoDiffBlock.hpp
View File

@ -24,7 +24,7 @@
#include <Eigen/Eigen>
#include <Eigen/Sparse>
#include <opm/autodiff/fastSparseProduct.hpp>
#include <opm/autodiff/fastSparseOperations.hpp>
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
@ -342,7 +342,8 @@ namespace Opm
jac[block] = D2*jac_[block];
}
else {
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));

18
opm/autodiff/AutoDiffMatrix.hpp
View File

@ -28,7 +28,7 @@
#include <opm/common/utility/platform_dependent/reenable_warnings.h>
#include <opm/common/ErrorMacros.hpp>
#include <opm/autodiff/fastSparseProduct.hpp>
#include <opm/autodiff/fastSparseOperations.hpp>
#include <vector>
@ -264,7 +264,13 @@ namespace Opm
AutoDiffMatrix& operator+=(const AutoDiffMatrix& rhs)
{
*this = *this + rhs;
if( type_ == Sparse && rhs.type_ == Sparse )
{
fastSparseAdd( sparse_, rhs.sparse_ );
}
else {
*this = *this + rhs;
}
return *this;
}
@ -275,7 +281,13 @@ namespace Opm
AutoDiffMatrix& operator-=(const AutoDiffMatrix& rhs)
{
*this = *this + (rhs * -1.0);
if( type_ == Sparse && rhs.type_ == Sparse )
{
fastSparseSubstract( sparse_, rhs.sparse_ );
}
else {
*this = *this + (rhs * -1.0);
}
return *this;
}

4
opm/autodiff/BlackoilModelBase_impl.hpp
View File

@ -1244,11 +1244,11 @@ namespace detail {
if (phase == Oil && active_[Gas]) {
const int gaspos = pu.phase_pos[Gas];
tmp = tmp - rv_perfcells * cmix_s[gaspos] / d;
tmp -= rv_perfcells * cmix_s[gaspos] / d;
}
if (phase == Gas && active_[Oil]) {
const int oilpos = pu.phase_pos[Oil];
tmp = tmp - rs_perfcells * cmix_s[oilpos] / d;
tmp -= rs_perfcells * cmix_s[oilpos] / d;
}
volumeRatio += tmp / b_perfcells[phase];
}

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);

99
opm/autodiff/fastSparseProduct.hpp → opm/autodiff/fastSparseOperations.hpp
View File

@ -177,7 +177,106 @@ inline void fastSparseDiagProduct(const Eigen::SparseMatrix<double>& lhs,
}
}
template<typename Lhs, typename Rhs>
inline bool
equalSparsityPattern(const Lhs& lhs, const Rhs& rhs)
{
// if both matrices have equal storage and non zeros match, we can check sparsity pattern
bool equal = (Lhs::IsRowMajor == Rhs::IsRowMajor) && (lhs.nonZeros() == rhs.nonZeros());
// check complete sparsity pattern
if( equal )
{
typedef typename Eigen::internal::remove_all<Lhs>::type::Index Index;
const Index outerSize = lhs.outerSize();
if( outerSize != rhs.outerSize() )
{
return false;
}
// outer indices
const Index* rhsOuter = rhs.outerIndexPtr();
const Index* lhsOuter = lhs.outerIndexPtr();
for(Index i=0; i<=outerSize; ++i )
{
if( lhsOuter[ i ] != rhsOuter[ i ] ) {
return false ;
}
}
// inner indices
const Index* rhsInner = rhs.innerIndexPtr();
const Index* lhsInner = lhs.innerIndexPtr();
const Index nnz = lhs.nonZeros();
for( Index i=0; i<nnz; ++i)
{
if( lhsInner[ i ] != rhsInner[ i ] ) {
return false;
}
}
}
return equal;
}
// this function substracts two sparse matrices
// if the sparsity pattern is the same a faster add/substract is performed
template<typename Lhs, typename Rhs>
inline void
fastSparseAdd(Lhs& lhs, const Rhs& rhs)
{
if( equalSparsityPattern( lhs, rhs ) )
{
typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
typedef typename Eigen::internal::remove_all<Lhs>::type::Index Index;
const Index nnz = lhs.nonZeros();
// fast add using only the data pointers
const Scalar* rhsV = rhs.valuePtr();
Scalar* lhsV = lhs.valuePtr();
for(Index i=0; i<nnz; ++i )
{
lhsV[ i ] += rhsV[ i ];
}
}
else
{
// default Eigen operator+=
lhs = lhs + rhs;
}
}
// this function substracts two sparse matrices
// if the sparsity pattern is the same a faster add/substract is performed
template<typename Lhs, typename Rhs>
inline void
fastSparseSubstract(Lhs& lhs, const Rhs& rhs)
{
if( equalSparsityPattern( lhs, rhs ) )
{
typedef typename Eigen::internal::remove_all<Lhs>::type::Scalar Scalar;
typedef typename Eigen::internal::remove_all<Lhs>::type::Index Index;
const Index nnz = lhs.nonZeros();
// fast add using only the data pointers
const Scalar* rhsV = rhs.valuePtr();
Scalar* lhsV = lhs.valuePtr();
for(Index i=0; i<nnz; ++i )
{
lhsV[ i ] -= rhsV[ i ];
}
}
else
{
// default Eigen operator-=
lhs = lhs - rhs;
}
}
} // end namespace Opm

1
opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.cpp
View File

@ -99,7 +99,6 @@ namespace {
}
typedef AutoDiffBlock<double>::V V;
typedef AutoDiffBlock<double>::M M;
const V& gpot = geo.gravityPotential();
const V& trans = geo.transmissibility();