OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Separate handling of NR vectors out to new class NewtonVectorCollection.

Browse Source 
The class is parametrised on policies for setting the size of
individual vectors, adding two vectors, and assembling local
contributions into the residual vector.

Provide default implementations of these policies, suitable for
base-vectors that implement (some of) the std::vector<T> interface.
This commit is contained in:
Bård Skaflestad 2011-10-04 21:27:57 +02:00
parent 5207b20058
commit 321f84b4f6
1 changed files with 98 additions and 106 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

204
src/JacobianSystem.hpp
View File

@ -36,149 +36,136 @@
#ifndef OPM_JACOBIANSYSTEM_HPP_HEADER #ifndef OPM_JACOBIANSYSTEM_HPP_HEADER
#define OPM_JACOBIANSYSTEM_HPP_HEADER #define OPM_JACOBIANSYSTEM_HPP_HEADER
#include <algorithm> #include <cassert>
#include <array>
#include <cmath> #include <cmath>
#include <cstddef> #include <cstddef>
#include <algorithm>
#include <array>
#include <functional> #include <functional>
#include <numeric> #include <numeric>
namespace Opm { namespace Opm {
namespace ImplicitTransportDefault { namespace ImplicitTransportDefault {
template <typename T> template <class BaseVec>
class MaxAbs : public std::binary_function<T, T, T> { class VectorAdder {
public: public:
T operator()(const T& x, const T& y) { // y += x
return std::max(std::abs(x), std::abs(y)); static void
add(const BaseVec& x, BaseVec& y) {
typedef typename BaseVec::value_type VT;
::std::transform(x.begin(), x.end(),
y.begin(),
y.begin(),
::std::plus<VT>());
} }
}; };
template <class Vector> template <class BaseVec>
class MaxNorm { class VectorNegater {
public: public:
typename Vector::value_type // x *= -1
norm(const Vector& v) { static void
typedef typename Vector::value_type VT; negate(BaseVec& x) {
typedef typename BaseVec::value_type VT;
return std::accumulate(v.begin(), v.end(), ::std::transform(x.begin(), x.end(),
VT(0), MaxAbs<VT>()); x.begin(),
::std::negate<VT>());
} }
}; };
template <typename T> template <class BaseVec>
class SumAbs : public std::binary_function<T, T, T> { class VectorZero {
public: static void
T operator()(const T& x, const T& y) { zero(BaseVec& x) {
return std::abs(x) + std::abs(y); typedef typename BaseVec::value_type VT;
::std::fill(x.begin(), x.end(), VT(0.0));
} }
}; };
template <class Vector> template <class BaseVec>
class TaxiCabNorm { class VectorBlockAssembler {
public: public:
typename Vector::value_type template <class Block>
norm(const Vector& v) { static void
typedef typename Vector::value_type VT; assemble(::std::size_t ndof,
::std::size_t i ,
const Block& b ,
BaseVec& vec ) {
return std::accumulate(v.begin(), v.end(), for (::std::size_t d = 0; d < ndof; ++d) {
VT(0), SumAbs<VT>()); vec[i*ndof + d] += b[d];
}
};
template <class Vector>
class EuclidianNorm {
public:
typename Vector::value_type
norm(const Vector& v) {
typedef typename Vector::value_type VT;
typedef typename Vector::iterator VI;
VT ret2 = 0;
for (VI i = v.begin(), e = v.end(); i != e; ++i) {
VT x = std::abs(*i);
ret2 += x * x;
} }
return std::sqrt(ret2);
} }
}; };
template <class Vector, template <class BaseVec>
template<class> class Norm = MaxNorm> class VectorSizeSetter {
class DefaultNewtonVector {
typedef typename Vector::value_type VT;
public: public:
void VectorSizeSetter(BaseVec& v) : v_(v) {}
resize(size_t m) { v_.resize(m); }
void void
negate() { setSize(::std::size_t ndof, ::std::size_t m) {
std::transform(v_.begin(), v_.resize(ndof * m);
v_.end (),
v_.begin(),
std::negate<VT>());
} }
typename Vector::value_type
norm() const { return Norm<Vector>::norm(v_); }
typename Vector::value_type&
operator[](size_t i) { return v_[i]; }
const typename Vector::value_type&
operator[](size_t i) const { return v_[i]; }
Vector& wrappedVector() { return v_; }
const Vector& wrappedVector() const { return v_; }
private: private:
Vector v_; BaseVec& v_;
}; };
template <class Vector> template <class BaseVec ,
class NewtonVectors { template <class> class VSzSetter = VectorSizeSetter ,
typedef typename std::array<Vector, 3> VC ; template <class> class VAdd = VectorAdder ,
typedef typename VC ::iterator VCI; template <class> class VBlkAsm = VectorBlockAssembler>
class NewtonVectorCollection {
enum { Residual = 0, Increment = 1, Solution = 2 }; enum { Residual = 0, Increment = 1, Solution = 2 };
public: public:
void void
setSize(size_t m) { setSize(::std::size_t ndof, ::std::size_t m) {
for (VCI i = v_.begin(), e = v_.end(); i != e; ++i) { typedef typename ::std::array<BaseVec, 3>::iterator VAI;
i->resize(m);
}
}
template <class Block> for (VAI i = vcoll_.begin(), e = vcoll_.end(); i != e; ++i) {
void VSzSetter<BaseVec>(*i).setSize(ndof, m);
assembleBlock(size_t n, size_t i, const Block& b) {
for (size_t k = 0; k < n; ++k) {
residual()[i*n + k] += b[k];
} }
ndof_ = ndof;
} }
void void
addIncrement() { addIncrement() {
std::transform(v_[ Solution ].begin(), VAdd<BaseVec>::add(vcoll_[ Increment ], vcoll_[ Solution ]);
v_[ Solution ].end (),
v_[ Increment ].begin(),
v_[ Solution ].begin(),
std::plus<typename Vector::value_type>());
} }
Vector& residual () { return v_[ Residual ]; } template <class Block>
const Vector& residual () const { return v_[ Residual ]; } void
assembleBlock(::std::size_t ndof,
::std::size_t i ,
const Block& b ) {
Vector& increment() { return v_[ Increment ]; } assert (ndof_ > 0 );
const Vector& increment() const { return v_[ Increment ]; } assert (ndof == ndof_);
Vector& solution () { return v_[ Solution ]; } VBlkAsm<BaseVec>::assemble(ndof, i, b, vcoll_[ Residual ]);
const Vector& solution () const { return v_[ Solution ]; } }
typedef BaseVec vector_type;
const vector_type& increment() const { return vcoll_[ Increment ]; }
const vector_type& residual () const { return vcoll_[ Residual ]; }
const vector_type& solution () const { return vcoll_[ Solution ]; }
// Write access for Newton solver purposes
vector_type& writableIncrement() { return vcoll_[ Increment ]; }
vector_type& writableResidual () { return vcoll_[ Residual ]; }
vector_type& writableSolution () { return vcoll_[ Solution ]; }
private: private:
VC v_; ::std::size_t ndof_ ;
::std::array<BaseVec, 3> vcoll_;
}; };
template <class Matrix> template <class Matrix>
@ -203,30 +190,35 @@ namespace Opm {
matrix(); matrix();
} */; } */;
template <class Matrix, class Vector, class LinSolve> template <class Matrix ,
class NVecCollection>
class JacobianSystem { class JacobianSystem {
public: public:
JacobianSystem(LinSolve& solver) : solver_(solver) {} JacobianSystem() {}
typedef Matrix matrix_type; typedef Matrix matrix_type;
typedef MatrixBlockAssembler<Matrix> assembler_type;
typedef typename NVecCollection::vector_type vector_type;
MatrixBlockAssembler<Matrix>& matrix() { return mba_ ; } assembler_type& matasm() { return mba_ ; }
NewtonVectors <Vector>& vector() { return vecs_; } NVecCollection& vector() { return sysvec_ ; }
const matrix_type& matrix() const { return mba_.matrix(); }
void void
solve() { setSize(::std::size_t ndof,
solver_.solve(mba_ .matrix (), ::std::size_t m ,
vecs_.residual (), ::std::size_t nnz = 0) {
vecs_.increment());
mba_ .setSize(ndof, m, m, nnz);
sysvec_.setSize(ndof, m );
} }
private: private:
JacobianSystem (const JacobianSystem&); JacobianSystem (const JacobianSystem&);
JacobianSystem& operator=(const JacobianSystem&); JacobianSystem& operator=(const JacobianSystem&);
MatrixBlockAssembler<Matrix> mba_ ; MatrixBlockAssembler<Matrix> mba_ ; // Coefficient matrix
NewtonVectors <Vector> vecs_ ; NVecCollection sysvec_; // Residual
LinSolve& solver_;
}; };
} }
} }