summaryrefslogtreecommitdiff
path: root/src/core/cell.h
blob: 003b1e8cd305ffd8d0646f4f24b7b3d3542bbb63 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2006-2007 Jonas Latt,
 *                2015-2019 Mathias J. Krause, Adrian Kummerlaender
 *  E-mail contact: info@openlb.net
 *  The most recent release of OpenLB can be downloaded at
 *  <http://www.openlb.net/>
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public
 *  License along with this program; if not, write to the Free
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/

/** \file
 * Definition of a LB cell -- header file.
 */
#ifndef CELL_H
#define CELL_H

#include "olbDebug.h"
#include "serializer.h"
#include "dynamics/latticeDescriptors.h"
#include "dynamics/dynamics.h"

namespace olb {

template<typename T, typename DESCRIPTORBASE>
class CellBase {
protected:
  /// The lattice populations and fields are stored as a C-array.
  T data[DESCRIPTORBASE::size()]; ///< distribution functions and additional fields

public:
  /// Read-write access to distribution functions.
  /**
   * \param iPop index of the accessed distribution function
   *
   * Note that for legacy-purposes this method allows access to all data of the cell
   * i.e. not just its distribution but also all fields that are declared by the
   * descriptor. Use of e.g. Cell::getFieldPointer or Cell::defineField is strongly
   * recommended when writing new code.
   **/
  T& operator[](int const& iPop)
  {
    OLB_PRECONDITION( iPop < DESCRIPTORBASE::size() );
    return data[iPop];
  }

  /// Read-only access to distribution functions.
  /**
   * \param iPop index of the accessed distribution function
   **/
  T const& operator[](int const& iPop) const
  {
    OLB_PRECONDITION( iPop < DESCRIPTORBASE::size() );
    return data[iPop];
  }
};


/// A LB lattice cell.
/**
 * A cell contains the q values of the distribution functions f on one lattice
 * point, additional "external" fields as well as a pointer to the dynamics of
 * the cell. Thanks to this pointer, one can have a space dependend definition
 * of the dynamics. This mechanism is useful e.g. for the implementation of
 * boundary conditions, or an inhomogeneous body force.
 *
 * The dynamics object is not owned by the class and as such not destructed in
 * the Cell destructor.
 *
 * This class is not intended to be derived from.
 */
template<typename T, typename DESCRIPTOR>
class Cell :
  public CellBase<T, typename DESCRIPTOR::BaseDescriptor>,
  public Serializable {
private:
  Dynamics<T,DESCRIPTOR>* dynamics;  ///< local LB dynamics

public:
  /// Default constructor.
  Cell();
  /// Constructor, to be used whenever possible.
  Cell(Dynamics<T,DESCRIPTOR>* dynamics_);

  /// Return pointer to FIELD of cell
  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<X::template provides<FIELD>(), T*>
  getFieldPointer()
  {
    const int offset = DESCRIPTOR::template index<FIELD>();
    return &(this->data[offset]);
  }

  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<!X::template provides<FIELD>(), T*>
  getFieldPointer()
  {
    throw std::invalid_argument("DESCRIPTOR does not provide FIELD.");
    return nullptr;
  }

  /// Return read-only pointer to FIELD of cell
  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<X::template provides<FIELD>(), const T*>
  getFieldPointer() const
  {
    const int offset = DESCRIPTOR::template index<FIELD>();
    return &(this->data[offset]);
  }

  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<!X::template provides<FIELD>(), const T*>
  getFieldPointer() const
  {
    throw std::invalid_argument("DESCRIPTOR does not provide FIELD.");
    return nullptr;
  }

  /// Return copy of FIELD as a vector
  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<(X::template size<FIELD>() > 1), Vector<T,X::template size<FIELD>()>>
  getField() const
  {
    return Vector<T,DESCRIPTOR::template size<FIELD>()>(
      getFieldPointer<FIELD>()
    );
  }

  /// Return copy of FIELD as a scalar
  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<(X::template size<FIELD>() == 1), T>
  getField() const
  {
    return getFieldPointer<FIELD>()[0];
  }

  /// Set value of FIELD from a vector
  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<(X::template size<FIELD>() > 1), void>
  setField(const Vector<T,DESCRIPTOR::template size<FIELD>()>& field)
  {
    std::copy_n(
      field.data,
      DESCRIPTOR::template size<FIELD>(),
      getFieldPointer<FIELD>());
  }

  /// Set value of FIELD from a scalar
  template <typename FIELD, typename X = DESCRIPTOR>
  utilities::meta::enable_if_t<(X::template size<FIELD>() == 1), void>
  setField(T value)
  {
    getFieldPointer<FIELD>()[0] = value;
  }

  /// Copy FIELD content to given memory location
  template <typename FIELD>
  void computeField(T* ext) const
  {
    const T* field = getFieldPointer<FIELD>();
    for (int iExt=0; iExt < DESCRIPTOR::template size<FIELD>(); ++iExt) {
      ext[iExt] = field[iExt];
    }
  }

  /// Set FIELD value from given memory location
  template <typename FIELD>
  void defineField(const T* ext)
  {
    T* field = getFieldPointer<FIELD>();
    for (int iExt=0; iExt < DESCRIPTOR::template size<FIELD>(); ++iExt) {
      field[iExt] = ext[iExt];
    }
  }

  /// Add to FIELD from given memory location
  /**
   * Similar to defineField(),but instead of replacing existing values
   * the data at ext is added onto the existing values.
   **/
  template <typename FIELD>
  inline void addField(const T* ext)
  {
    T* field = getFieldPointer<FIELD>();
    for (int iExt=0; iExt < DESCRIPTOR::template size<FIELD>(); ++iExt) {
      field[iExt] += ext[iExt];
    }
  }

  /// Multiply FIELD with values at given memory location
  /**
   * Similar to defineField(), but instead of replacing existing values
   * the data at ext is multiplied to the existing values.
   **/
  template <typename FIELD>
  inline void multiplyField(const T* ext)
  {
    T* field = getFieldPointer<FIELD>();
    for (int iExt=0; iExt < DESCRIPTOR::template size<FIELD>(); ++iExt) {
      field[iExt] *= ext[iExt];
    }
  }

  /// Define or re-define dynamics of the cell.
  /**
   * \param dynamics_ a pointer to the dynamics object, whos memory management
   *                  falls under the responsibility of the user
   **/
  void defineDynamics(Dynamics<T,DESCRIPTOR>* dynamics_);
  /// Get a non-modifiable pointer to the dynamics
  Dynamics<T,DESCRIPTOR> const* getDynamics() const;
  /// Get a non-modifiable pointer to the dynamics
  Dynamics<T,DESCRIPTOR>* getDynamics();

  // The following helper functions forward the function call
  // to the Dynamics object
public:
  /// Apply LB collision to the cell according to local dynamics.
  void collide(LatticeStatistics<T>& statistics)
  {
    OLB_PRECONDITION( dynamics );
    dynamics->collide(*this, statistics);
  }

  /// Compute particle density on the cell.
  /** \return particle density
   */
  T computeRho() const
  {
    OLB_PRECONDITION( dynamics );
    return dynamics->computeRho(*this);
  }
  /// Compute fluid velocity on the cell.
  /** \param u fluid velocity
   */
  void computeU(T u[descriptors::d<DESCRIPTOR>()]) const
  {
    OLB_PRECONDITION( dynamics );
    dynamics->computeU(*this, u);
  }
  /// Compute fluid momentum (j = rho * u) on the cell.
  /** \param j fluid momentum
   */
  void computeJ(T j[descriptors::d<DESCRIPTOR>()]) const
  {
    OLB_PRECONDITION( dynamics );
    dynamics->computeJ(*this, j);
  }
  /// Compute components of the stress tensor on the cell.
  /** \param pi stress tensor */
  void computeStress (
    T pi[util::TensorVal<DESCRIPTOR >::n]) const
  {
    OLB_PRECONDITION( dynamics );
    T rho, u[descriptors::d<DESCRIPTOR>()];
    dynamics->computeRhoU(*this, rho, u);
    dynamics->computeStress(*this, rho, u, pi);
  }
  /// Compute fluid velocity and particle density on the cell.
  /** \param rho particle density
   *  \param u fluid velocity
   */
  void computeRhoU(T& rho, T u[descriptors::d<DESCRIPTOR>()]) const
  {
    OLB_PRECONDITION( dynamics );
    dynamics->computeRhoU(