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

NewtonBlackoilInterleaved: add solver option for single precision.

Browse Source
This commit is contained in:
Robert Kloefkorn 2016-02-08 13:28:28 +01:00
parent 37dc877074
commit 31812bd2f9
2 changed files with 47 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
opm/autodiff/NewtonIterationBlackoilInterleaved.cpp
View File

@ -49,6 +49,22 @@
#include <opm/common/utility/platform_dependent/reenable_warnings.h> #include <opm/common/utility/platform_dependent/reenable_warnings.h>
namespace Dune
{
namespace FMatrixHelp {
//! invert Matrix without changing the original matrix
//! return inverse matrix (specialization for n <= 3 in dune/common/fmatrix.hh)
template <typename K, int n>
static inline K invertMatrix (const FieldMatrix<K,n,n> &matrix, FieldMatrix<K,n,n> &inverse)
{
inverse.invert();
return K(0);
}
}
}
namespace Opm namespace Opm
{ {
@ -70,12 +86,28 @@ namespace Opm
/// solving the reduced system (after eliminating well variables) /// solving the reduced system (after eliminating well variables)
/// as a block-structured matrix (one block for all cell variables) for a fixed /// as a block-structured matrix (one block for all cell variables) for a fixed
/// number of cell variables np . /// number of cell variables np .
template <int np> template <int np, class ScalarT = double >
class NewtonIterationBlackoilInterleavedImpl : public NewtonIterationBlackoilInterface class NewtonIterationBlackoilInterleavedImpl : public NewtonIterationBlackoilInterface
{ {
typedef double Scalar; typedef ScalarT Scalar;
typedef Dune::FieldVector<Scalar, np > VectorBlockType; typedef Dune::FieldVector<Scalar, np > VectorBlockType;
typedef Dune::FieldMatrix<Scalar, np, np> MatrixBlockType;
class MatrixBlock : public Dune::FieldMatrix<Scalar, np, np>
{
typedef Dune::FieldMatrix<Scalar, np, np> BaseType;
public:
using BaseType :: operator= ;
using BaseType :: rows;
using BaseType :: cols;
MatrixBlock() : BaseType( Scalar( 0 ) ) {}
void invert()
{
MatrixBlock matrix( *this );
Dune::FMatrixHelp::invertMatrix(matrix, *this);
}
};
typedef MatrixBlock MatrixBlockType;
typedef Dune::BCRSMatrix <MatrixBlockType> Mat; typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
typedef Dune::BlockVector<VectorBlockType> Vector; typedef Dune::BlockVector<VectorBlockType> Vector;
@ -406,7 +438,7 @@ namespace Opm
namespace detail { namespace detail {
template< int NP > template< int NP, class Scalar >
struct NewtonIncrement struct NewtonIncrement
{ {
template <class NewtonIncVector> template <class NewtonIncVector>
@ -421,18 +453,18 @@ namespace Opm
assert( np < int(newtonIncrements.size()) ); assert( np < int(newtonIncrements.size()) );
// create NewtonIncrement with fixed np // create NewtonIncrement with fixed np
if( ! newtonIncrements[ NP ] ) if( ! newtonIncrements[ NP ] )
newtonIncrements[ NP ].reset( new NewtonIterationBlackoilInterleavedImpl< NP >( param, parallelInformation ) ); newtonIncrements[ NP ].reset( new NewtonIterationBlackoilInterleavedImpl< NP, Scalar >( param, parallelInformation ) );
return *(newtonIncrements[ NP ]); return *(newtonIncrements[ NP ]);
} }
else else
{ {
return NewtonIncrement< NP-1 >::get(newtonIncrements, param, parallelInformation, np ); return NewtonIncrement< NP-1, Scalar >::get(newtonIncrements, param, parallelInformation, np );
} }
} }
}; };
template<> template<class Scalar>
struct NewtonIncrement< 0 > struct NewtonIncrement< 0, Scalar >
{ {
template <class NewtonIncVector> template <class NewtonIncVector>
static const NewtonIterationBlackoilInterface& static const NewtonIterationBlackoilInterface&
@ -455,9 +487,12 @@ namespace Opm
const int np = residual.material_balance_eq.size(); const int np = residual.material_balance_eq.size();
// maxNumberEquations_ denotes the currently maximal number of equations // maxNumberEquations_ denotes the currently maximal number of equations
// covered, this is mostly to reduce compile time. Adjust accordingly to cover // covered, this is mostly to reduce compile time. Adjust accordingly to cover
const bool singlePrecision = false ; // residual.singlePrecision ;
// more cases // more cases
const NewtonIterationBlackoilInterface& newtonIncrement = const NewtonIterationBlackoilInterface& newtonIncrement = singlePrecision ?
detail::NewtonIncrement< maxNumberEquations_ > :: get( newtonIncrement_, parameters_, parallelInformation_, np ); detail::NewtonIncrement< maxNumberEquations_, float > :: get( newtonIncrementSinglePrecision_, parameters_, parallelInformation_, np ) :
detail::NewtonIncrement< maxNumberEquations_, double > :: get( newtonIncrement_, parameters_, parallelInformation_, np );
// compute newton increment // compute newton increment
SolutionVector dx = newtonIncrement.computeNewtonIncrement( residual ); SolutionVector dx = newtonIncrement.computeNewtonIncrement( residual );

1
opm/autodiff/NewtonIterationBlackoilInterleaved.hpp
View File

@ -101,6 +101,7 @@ namespace Opm
static const int maxNumberEquations_ = 6 ; static const int maxNumberEquations_ = 6 ;
mutable std::array< std::unique_ptr< NewtonIterationBlackoilInterface >, maxNumberEquations_+1 > newtonIncrement_; mutable std::array< std::unique_ptr< NewtonIterationBlackoilInterface >, maxNumberEquations_+1 > newtonIncrement_;
mutable std::array< std::unique_ptr< NewtonIterationBlackoilInterface >, maxNumberEquations_+1 > newtonIncrementSinglePrecision_;
NewtonIterationBlackoilInterleavedParameters parameters_; NewtonIterationBlackoilInterleavedParameters parameters_;
boost::any parallelInformation_; boost::any parallelInformation_;
mutable int iterations_; mutable int iterations_;