ImplicitAssembly.hpp
Go to the documentation of this file.
1 /*===========================================================================
2 //
3 // File: ImplicitAssembly.hpp
4 //
5 // Created: 2011-09-28 10:00:46+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_IMPLICITASSEMBLY_HPP_HEADER
37 #define OPM_IMPLICITASSEMBLY_HPP_HEADER
38 
39 #include <cstddef>
40 #include <algorithm>
41 #include <vector>
42 
43 namespace Opm {
44  template <class Model>
45  class ImplicitAssembly {
46  enum { DofPerCell = Model::DofPerCell };
47 
48  public:
49  ImplicitAssembly(Model& model)
50  : model_(model),
51  nconn_(-1) ,
52  asm_buffer_()
53  {}
54 
55  template <class Grid ,
56  class JacobianSystem>
57  void
58  createSystem(const Grid& g ,//const SourceTerms* tsrc,
59  JacobianSystem& sys) const {
60 
61  ::std::size_t m = g.number_of_cells;
62  ::std::size_t nnz = g.number_of_cells
63  + countConnections(g);
64 // + tsrc->pf;
65 
66  sys.setSize(DofPerCell, m, nnz);
67 
68  for (int c = 0; c < g.number_of_cells; ++c) {
69  this->createRowStructure(g, c, sys);
70  }
71  /* is done by modifiying grid
72  // add elements for periodic boundary
73  for (int pf = 0; pf<tsrc->pf; ++pf){
74  int nconn = 1;
75  ::std::vector<int> connections(1);
76  connections[0]=(tsrc->periodic_cells[2*pf]);
77  sys.matasm().createBlockRow(tsrc->periodic_cells[2*pf+1], connections, DofPerCell);
78  connections[0]=(tsrc->periodic_cells[2*pf+1]);
79  sys.matasm().createBlockRow(tsrc->periodic_cells[2*pf], connections, DofPerCell);
80  }
81  */
82  sys.matasm().finalizeStructure();
83  }
84 
85  template <class ReservoirState,
86  class Grid ,
87  class SourceTerms ,
88  class JacobianSystem>
89  void
90  assemble(const ReservoirState& state,
91  const Grid& g ,
92  const SourceTerms* src ,
93  const double dt ,
94  JacobianSystem& sys ) {
95 
96  for (int c = 0; c < g.number_of_cells; ++c) {
97  this->computeCellContrib(state, g, dt, c);
98  this->assembleCellContrib(g, c, sys);
99  }
100 
101  if (src != 0) {
102  this->assembleSourceContrib(g, src, dt, sys);
103  }
104  }
105 
106  private:
107  template <class Grid>
108  int
109  countConnections(const Grid& g, int c) const {
110  int i, f, c1, c2, n;
111 
112  for (i = g.cell_facepos[c + 0], n = 0;
113  i < g.cell_facepos[c + 1]; ++i) {
114  f = g.cell_faces[i];
115  c1 = g.face_cells[2*f + 0];
116  c2 = g.face_cells[2*f + 1];
117 
118  n += (c1 >= 0) && (c2 >= 0);
119  }
120 
121  return n;
122  }
123 
124  template <class Grid>
125  int
126  countConnections(const Grid& g) const {
127  int n = 0;
128 
129  for (int c = 0; c < g.number_of_cells; ++c) {
130  n += countConnections(g, c);
131  }
132 
133  return n;
134  }
135 
136  template <class Grid, class System>
137  void
138  createRowStructure(const Grid& g ,
139  const int c ,
140  System& sys) const {
141 
142  int nconn = countConnections(g, c);
143 
144  ::std::vector<int> connections;
145  connections.reserve (nconn + 1);
146  connections.push_back(c);
147 
148  for (int i = g.cell_facepos[c + 0];
149  i < g.cell_facepos[c + 1]; ++i) {
150  int f = g.cell_faces[i];
151  int c1 = g.face_cells[2*f + 0];
152  int c2 = g.face_cells[2*f + 1];
153 
154  if ((c1 >= 0) && (c2 >= 0)) {
155  connections.push_back((c1 == c) ? c2 : c1);
156  }
157  }
158  sys.matasm().createBlockRow(c, connections, DofPerCell);
159 
160  }
161 
162  template <class ReservoirState, class Grid>
163  int
164  computeCellContrib(const ReservoirState& state,
165  const Grid& g ,
166  const double dt ,
167  const int c ) {
168  const int ndof = DofPerCell;
169  const int ndof2 = ndof * ndof;
170 
171  nconn_ = countConnections(g, c);
172 
173  asm_buffer_.resize((2*nconn_ + 1)*ndof2 + (nconn_ + 2)*ndof);
174  std::fill(asm_buffer_.begin(), asm_buffer_.end(), 0.0);
175 
176  double* F = &asm_buffer_[(2*nconn_ + 1) * ndof2];
177  double* J1 = &asm_buffer_[(0*nconn_ + 1) * ndof2];
178  double* J2 = J1 + (1*nconn_ + 0) * ndof2 ;
179 
180  model_.initResidual(c, F);
181  F += ndof;
182 
183  for (int i = g.cell_facepos[c + 0];
184  i < g.cell_facepos[c + 1]; ++i) {
185  int f = g.cell_faces[i];
186  int c1 = g.face_cells[2*f + 0];
187  int c2 = g.face_cells[2*f + 1];
188 
189  if ((c1 >= 0) && (c2 >= 0)) {
190  model_.fluxConnection(state, g, dt, c, f, J1, J2, F);
191  J1 += ndof2; J2 += ndof2; F += ndof;
192  }
193  }
194 
195  model_.accumulation(g, c, &asm_buffer_[0], F);
196 
197  return 1;
198  }
199 
200  template <class Grid, class System>
201  void
202  assembleCellContrib(const Grid& g ,
203  const int c ,
204  System& sys) const {
205  const int ndof = DofPerCell;
206  const int ndof2 = ndof * ndof;
207 
208  const double* J1 = &asm_buffer_[0];
209  const double* J2 = J1 + ((1*nconn_ + 1) * ndof2);
210 
211  // Assemble contributions from accumulation term
212  sys.matasm().assembleBlock(ndof, c, c, J1); J1 += ndof2;
213 
214  // Assemble connection contributions.
215  for (int i = g.cell_facepos[c + 0];
216  i < g.cell_facepos[c + 1]; ++i) {
217  int f = g.cell_faces[i];
218  int c1 = g.face_cells[2*f + 0];
219  int c2 = g.face_cells[2*f + 1];
220 
221  c2 = (c1 == c) ? c2 : c1;
222 
223  if (c2 >= 0) {
224  sys.matasm().assembleBlock(ndof, c, c , J1);
225  sys.matasm().assembleBlock(ndof, c, c2, J2);
226 
227  J1 += ndof2;
228  J2 += ndof2;
229  }
230  }
231  // Assemble residual
232  const double* F = &asm_buffer_[(2*nconn_ + 1) * ndof2];
233  for (int conn = 0; conn < nconn_ + 2; ++conn, F += ndof) {
234  sys.vector().assembleBlock(ndof, c, F);
235  }
236  }
237 
238  template <class Grid, class SourceTerms, class System>
239  void
240  assembleSourceContrib(const Grid& g,
241  const SourceTerms* src,
242  const double dt,
243  System& sys) {
244  const int ndof = DofPerCell;
245  const int ndof2 = ndof * ndof;
246 
247  for (int i = 0; i < src->nsrc; ++i) {
248  std::fill_n(asm_buffer_.begin(), ndof2 + ndof, 0.0);
249 
250  double *J = &asm_buffer_[0];
251  double *F = J + ndof2;
252 
253  model_.sourceTerms(g, src, i, dt, J, F);
254 
255  const int c = src->cell[i];
256 
257  sys.matasm().assembleBlock(ndof, c, c, J);
258  sys.vector().assembleBlock(ndof, c, F);
259  }
260  }
261 
262  Model& model_ ;
263  int nconn_ ;
264  std::vector<double> asm_buffer_;
265  };
266 }
267 #endif /* OPM_IMPLICITASSEMBLY_HPP_HEADER */
Opm
Definition: AnisotropicEikonal.hpp:43
Opm::ImplicitAssembly::ImplicitAssembly
ImplicitAssembly(Model &model)
Definition: ImplicitAssembly.hpp:49
Opm::ImplicitAssembly::assemble
void assemble(const ReservoirState &state, const Grid &g, const SourceTerms *src, const double dt, JacobianSystem &sys)
Definition: ImplicitAssembly.hpp:90
Opm::ImplicitAssembly::createSystem
void createSystem(const Grid &g, JacobianSystem &sys) const
Definition: ImplicitAssembly.hpp:58
Opm::ImplicitAssembly
Definition: ImplicitAssembly.hpp:45