From 321f84b4f625b8200a83a70b784c72652702a8ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 4 Oct 2011 21:27:57 +0200 Subject: [PATCH] Separate handling of NR vectors out to new class NewtonVectorCollection. 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 interface. --- src/JacobianSystem.hpp | 204 ++++++++++++++++++++--------------------- 1 file changed, 98 insertions(+), 106 deletions(-) diff --git a/src/JacobianSystem.hpp b/src/JacobianSystem.hpp index a55f42c2..ab4bac22 100644 --- a/src/JacobianSystem.hpp +++ b/src/JacobianSystem.hpp @@ -36,149 +36,136 @@ #ifndef OPM_JACOBIANSYSTEM_HPP_HEADER #define OPM_JACOBIANSYSTEM_HPP_HEADER -#include -#include +#include #include #include + +#include +#include #include #include namespace Opm { namespace ImplicitTransportDefault { - template - class MaxAbs : public std::binary_function { + template + class VectorAdder { public: - T operator()(const T& x, const T& y) { - return std::max(std::abs(x), std::abs(y)); + // y += x + 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()); } }; - template - class MaxNorm { + template + class VectorNegater { public: - typename Vector::value_type - norm(const Vector& v) { - typedef typename Vector::value_type VT; + // x *= -1 + static void + negate(BaseVec& x) { + typedef typename BaseVec::value_type VT; - return std::accumulate(v.begin(), v.end(), - VT(0), MaxAbs()); + ::std::transform(x.begin(), x.end(), + x.begin(), + ::std::negate()); } }; - template - class SumAbs : public std::binary_function { - public: - T operator()(const T& x, const T& y) { - return std::abs(x) + std::abs(y); + template + class VectorZero { + static void + zero(BaseVec& x) { + typedef typename BaseVec::value_type VT; + + ::std::fill(x.begin(), x.end(), VT(0.0)); } }; - template - class TaxiCabNorm { + template + class VectorBlockAssembler { public: - typename Vector::value_type - norm(const Vector& v) { - typedef typename Vector::value_type VT; + template + static void + assemble(::std::size_t ndof, + ::std::size_t i , + const Block& b , + BaseVec& vec ) { - return std::accumulate(v.begin(), v.end(), - VT(0), SumAbs()); - } - }; - - template - 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; + for (::std::size_t d = 0; d < ndof; ++d) { + vec[i*ndof + d] += b[d]; } - - return std::sqrt(ret2); } }; - template class Norm = MaxNorm> - class DefaultNewtonVector { - typedef typename Vector::value_type VT; + template + class VectorSizeSetter { public: - void - resize(size_t m) { v_.resize(m); } + VectorSizeSetter(BaseVec& v) : v_(v) {} void - negate() { - std::transform(v_.begin(), - v_.end (), - v_.begin(), - std::negate()); + setSize(::std::size_t ndof, ::std::size_t m) { + v_.resize(ndof * m); } - typename Vector::value_type - norm() const { return Norm::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: - Vector v_; + BaseVec& v_; }; - template - class NewtonVectors { - typedef typename std::array VC ; - typedef typename VC ::iterator VCI; - + template class VSzSetter = VectorSizeSetter , + template class VAdd = VectorAdder , + template class VBlkAsm = VectorBlockAssembler> + class NewtonVectorCollection { enum { Residual = 0, Increment = 1, Solution = 2 }; public: void - setSize(size_t m) { - for (VCI i = v_.begin(), e = v_.end(); i != e; ++i) { - i->resize(m); - } - } + setSize(::std::size_t ndof, ::std::size_t m) { + typedef typename ::std::array::iterator VAI; - template - void - assembleBlock(size_t n, size_t i, const Block& b) { - for (size_t k = 0; k < n; ++k) { - residual()[i*n + k] += b[k]; + for (VAI i = vcoll_.begin(), e = vcoll_.end(); i != e; ++i) { + VSzSetter(*i).setSize(ndof, m); } + + ndof_ = ndof; } void addIncrement() { - std::transform(v_[ Solution ].begin(), - v_[ Solution ].end (), - v_[ Increment ].begin(), - v_[ Solution ].begin(), - std::plus()); + VAdd::add(vcoll_[ Increment ], vcoll_[ Solution ]); } - Vector& residual () { return v_[ Residual ]; } - const Vector& residual () const { return v_[ Residual ]; } + template + void + assembleBlock(::std::size_t ndof, + ::std::size_t i , + const Block& b ) { - Vector& increment() { return v_[ Increment ]; } - const Vector& increment() const { return v_[ Increment ]; } + assert (ndof_ > 0 ); + assert (ndof == ndof_); - Vector& solution () { return v_[ Solution ]; } - const Vector& solution () const { return v_[ Solution ]; } + VBlkAsm::assemble(ndof, i, b, vcoll_[ Residual ]); + } + + 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: - VC v_; + ::std::size_t ndof_ ; + ::std::array vcoll_; }; template @@ -203,30 +190,35 @@ namespace Opm { matrix(); } */; - template + template class JacobianSystem { public: - JacobianSystem(LinSolve& solver) : solver_(solver) {} + JacobianSystem() {} - typedef Matrix matrix_type; + typedef Matrix matrix_type; + typedef MatrixBlockAssembler assembler_type; + typedef typename NVecCollection::vector_type vector_type; - MatrixBlockAssembler& matrix() { return mba_ ; } - NewtonVectors & vector() { return vecs_; } + assembler_type& matasm() { return mba_ ; } + NVecCollection& vector() { return sysvec_ ; } + const matrix_type& matrix() const { return mba_.matrix(); } void - solve() { - solver_.solve(mba_ .matrix (), - vecs_.residual (), - vecs_.increment()); + setSize(::std::size_t ndof, + ::std::size_t m , + ::std::size_t nnz = 0) { + + mba_ .setSize(ndof, m, m, nnz); + sysvec_.setSize(ndof, m ); } private: JacobianSystem (const JacobianSystem&); JacobianSystem& operator=(const JacobianSystem&); - MatrixBlockAssembler mba_ ; - NewtonVectors vecs_ ; - LinSolve& solver_; + MatrixBlockAssembler mba_ ; // Coefficient matrix + NVecCollection sysvec_; // Residual }; } }