My Project
JacobianSystem.hpp
1 /*===========================================================================
2 //
3 // File: JacobianSystem.hpp
4 //
5 // Created: 2011-09-30 19:23:31+0200
6 //
7 // Authors: Ingeborg S. Ligaarden <Ingeborg.Ligaarden@sintef.no>
8 // Jostein R. Natvig <Jostein.R.Natvig@sintef.no>
9 // Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
10 // Atgeirr F. Rasmussen <atgeirr@sintef.no>
11 // Bård Skaflestad <Bard.Skaflestad@sintef.no>
12 //
13 //==========================================================================*/
14 
15 
16 /*
17  Copyright 2011 SINTEF ICT, Applied Mathematics.
18  Copyright 2011 Statoil ASA.
19 
20  This file is part of the Open Porous Media Project (OPM).
21 
22  OPM is free software: you can redistribute it and/or modify
23  it under the terms of the GNU General Public License as published by
24  the Free Software Foundation, either version 3 of the License, or
25  (at your option) any later version.
26 
27  OPM is distributed in the hope that it will be useful,
28  but WITHOUT ANY WARRANTY; without even the implied warranty of
29  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30  GNU General Public License for more details.
31 
32  You should have received a copy of the GNU General Public License
33  along with OPM. If not, see <http://www.gnu.org/licenses/>.
34 */
35 
36 #ifndef OPM_JACOBIANSYSTEM_HPP_HEADER
37 #define OPM_JACOBIANSYSTEM_HPP_HEADER
38 
39 #include <cassert>
40 #include <cmath>
41 #include <cstddef>
42 
43 #include <algorithm>
44 //#include <array>
45 #include <functional>
46 #include <numeric>
47 
48 namespace Opm {
49  namespace ImplicitTransportDefault {
50  template <class BaseVec>
51  class VectorAdder {
52  public:
53  // y += x
54  static void
55  add(const BaseVec& x, BaseVec& y) {
56  typedef typename BaseVec::value_type VT;
57  ::std::transform(x.begin(), x.end(),
58  y.begin(),
59  y.begin(),
60  ::std::plus<VT>());
61  }
62  };
63 
64  template <class BaseVec>
65  class VectorNegater {
66  public:
67  // x *= -1
68  static void
69  negate(BaseVec& x) {
70  typedef typename BaseVec::value_type VT;
71  ::std::transform(x.begin(), x.end(),
72  x.begin(),
73  ::std::negate<VT>());
74  }
75  };
76 
77  template <class BaseVec>
78  class VectorZero {
79  public:
80  static void
81  zero(BaseVec& x) {
82  typedef typename BaseVec::value_type VT;
83 
84  ::std::fill(x.begin(), x.end(), VT(0.0));
85  }
86  };
87 
88  template <class BaseVec>
89  class VectorAssign {
90  public:
91  // y <- x
92  static void
93  assign(const BaseVec& x, BaseVec& y) {
94  ::std::copy(x.begin(), x.end(), y.begin());
95  }
96 
97  // y <- a*x
98  template <class Scalar>
99  static void
100  assign(const Scalar& a, const BaseVec& x, BaseVec& y) {
101  ::std::transform(x.begin(), x.end(),
102  y.begin(),
103  [&a](auto z) { return a * z; });
104  }
105  };
106 
107  template <class Matrix>
108  class MatrixZero;
109 
110  template <class BaseVec>
112  public:
113  template <class Block>
114  static void
115  assemble(::std::size_t ndof,
116  ::std::size_t i ,
117  const Block& b ,
118  BaseVec& vec ) {
119 
120  for (::std::size_t d = 0; d < ndof; ++d) {
121  vec[i*ndof + d] += b[d];
122  }
123  }
124  };
125 
126  template <class BaseVec>
128  public:
129  VectorSizeSetter(BaseVec& v) : v_(v) {}
130 
131  void
132  setSize(::std::size_t ndof, ::std::size_t m) {
133  v_.resize(ndof * m);
134  }
135 
136  private:
137  BaseVec& v_;
138  };
139 
140  template <class BaseVec ,
141  template <class> class VSzSetter = VectorSizeSetter ,
142  template <class> class VAdd = VectorAdder ,
143  template <class> class VBlkAsm = VectorBlockAssembler>
145  enum { Residual = 0, Increment = 1, Solution = 2 };
146 
147  public:
148  void
149  setSize(::std::size_t ndof, ::std::size_t m) {
150 #if 0
151  typedef typename ::std::array<BaseVec, 3>::iterator VAI;
152 
153  for (VAI i = vcoll_.begin(), e = vcoll_.end(); i != e; ++i) {
154  VSzSetter<BaseVec>(*i).setSize(ndof, m);
155  }
156 #endif
157  for (::std::size_t i = 0; i < sizeof (vcoll_) / sizeof (vcoll_[0]); ++i) {
158  VSzSetter<BaseVec>(vcoll_[i]).setSize(ndof, m);
159  }
160 
161  ndof_ = ndof;
162  }
163 
164  void
165  addIncrement() {
166  VAdd<BaseVec>::add(vcoll_[ Increment ], vcoll_[ Solution ]);
167  }
168 
169  template <class Block>
170  void
171  assembleBlock(::std::size_t ndof,
172  ::std::size_t i ,
173  const Block& b ) {
174 
175  assert (ndof_ > 0 );
176  assert (ndof == ndof_);
177 
178  VBlkAsm<BaseVec>::assemble(ndof, i, b, vcoll_[ Residual ]);
179  }
180 
181  typedef BaseVec vector_type;
182 
183  const vector_type& increment() const { return vcoll_[ Increment ]; }
184  const vector_type& residual () const { return vcoll_[ Residual ]; }
185  const vector_type& solution () const { return vcoll_[ Solution ]; }
186 
187  // Write access for Newton solver purposes
188  vector_type& writableIncrement() { return vcoll_[ Increment ]; }
189  vector_type& writableResidual () { return vcoll_[ Residual ]; }
190  vector_type& writableSolution () { return vcoll_[ Solution ]; }
191 
192  private:
193  ::std::size_t ndof_ ;
194  BaseVec vcoll_[3];
195  //::std::array<BaseVec, 3> vcoll_;
196  };
197 
198  template <class Matrix>
199  class MatrixBlockAssembler
200  /* {
201  public:
202  template <class Block>
203  void
204  assembleBlock(size_t n, size_t i, size j, const Block& b);
205 
206  template <class Connections>
207  void
208  createBlockRow(size_t i, const Connections& conn, size_t n);
209 
210  void
211  finalizeStructure();
212 
213  void
214  setSize(size_t ndof, size_t m, size_t n, size_t nnz = 0);
215 
216  Matrix& matrix() ;
217  const Matrix& matrix() const;
218  } */;
219 
220  template <class Matrix ,
221  class NVecCollection>
223  public:
224  JacobianSystem() {}
225 
226  typedef Matrix matrix_type;
227  typedef MatrixBlockAssembler<Matrix> assembler_type;
228  typedef NVecCollection vector_collection_type;
229  typedef typename NVecCollection::vector_type vector_type;
230 
231  assembler_type& matasm() { return mba_ ; }
232  NVecCollection& vector() { return sysvec_ ; }
233  const matrix_type& matrix() const { return mba_.matrix(); }
234  matrix_type& writableMatrix() { return mba_.matrix(); }
235 
236  void
237  setSize(::std::size_t ndof,
238  ::std::size_t m ,
239  ::std::size_t nnz = 0) {
240 
241  mba_ .setSize(ndof, m, m, nnz);
242  sysvec_.setSize(ndof, m );
243  }
244 
245  private:
247  JacobianSystem& operator=(const JacobianSystem&);
248 
249  MatrixBlockAssembler<Matrix> mba_ ; // Coefficient matrix
250  NVecCollection sysvec_; // Residual
251  };
252  }
253 }
254 
255 #endif /* OPM_JACOBIANSYSTEM_HPP_HEADER */
Definition: JacobianSystem.hpp:222
Definition: JacobianSystem.hpp:108
Definition: JacobianSystem.hpp:51
Definition: JacobianSystem.hpp:89
Definition: JacobianSystem.hpp:111
Definition: JacobianSystem.hpp:65
Definition: JacobianSystem.hpp:127
Definition: JacobianSystem.hpp:78
Inverting small matrices.
Definition: ImplicitAssembly.hpp:43