Athena++/Atmosphere
Planetary Atmosphere Simulator
athena_arrays.hpp
Go to the documentation of this file.
1 #ifndef ATHENA_ARRAYS_HPP_
2 #define ATHENA_ARRAYS_HPP_
3 //========================================================================================
4 // Athena++ astrophysical MHD code
5 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
6 // Licensed under the 3-clause BSD License, see LICENSE file for details
7 //========================================================================================
9 // \brief provides array classes valid in 1D to 5D.
10 //
11 // The operator() is overloaded, e.g. elements of a 4D array of size [N4xN3xN2xN1]
12 // are accessed as: A(n,k,j,i) = A[i + N1*(j + N2*(k + N3*n))]
13 // NOTE THE TRAILING INDEX INSIDE THE PARENTHESES IS INDEXED FASTEST
14 
15 // C headers
16 
17 // C++ headers
18 #include <cstddef> // size_t
19 #include <cstring> // memset()
20 #include <utility> // swap()
21 
22 // Athena++ headers
23 #include "stride_iterator.hpp"
24 
25 template <typename T>
26 class AthenaArray {
27  public:
28  enum class DataStatus {empty, shallow_slice, allocated}; // formerly, "bool scopy_"
29  // ctors
30  // default ctor: simply set null AthenaArray
31  AthenaArray() : pdata_(nullptr), nx1_(0), nx2_(0), nx3_(0),
32  nx4_(0), nx5_(0), nx6_(0), state_(DataStatus::empty) {}
33  // ctor overloads: set expected size of unallocated container, maybe allocate (default)
35  pdata_(nullptr), nx1_(nx1), nx2_(1), nx3_(1), nx4_(1), nx5_(1), nx6_(1),
36  state_(init) { AllocateData(); }
38  pdata_(nullptr), nx1_(nx1), nx2_(nx2), nx3_(1), nx4_(1), nx5_(1), nx6_(1),
39  state_(init) { AllocateData(); }
41  pdata_(nullptr), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(1), nx5_(1), nx6_(1),
42  state_(init) { AllocateData(); }
44  pdata_(nullptr), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(nx4), nx5_(1), nx6_(1),
45  state_(init) { AllocateData(); }
46  AthenaArray(int nx5, int nx4, int nx3, int nx2, int nx1,
48  pdata_(nullptr), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(nx4), nx5_(nx5), nx6_(1),
49  state_(init) { AllocateData(); }
50  AthenaArray(int nx6, int nx5, int nx4, int nx3, int nx2, int nx1,
52  pdata_(nullptr), nx1_(nx1), nx2_(nx2), nx3_(nx3), nx4_(nx4), nx5_(nx5), nx6_(nx6),
53  state_(init) { AllocateData(); }
54  // still allowing delayed-initialization (after constructor) via array.NewAthenaArray()
55  // or array.InitWithShallowSlice() (only used in outputs.cpp + 3x other files)
56  // TODO(felker): replace InitWithShallowSlice with ??? and remove shallow_copy enum val
57  // TODO(felker): replace raw pointer with std::vector + reshape (if performance is same)
58 
59  // user-provided dtor, "rule of five" applies:
61  // define copy constructor and overload assignment operator so both do deep copies.
63  __attribute__((nothrow)) AthenaArray<T> &operator= (const AthenaArray<T> &t);
64  // define move constructor and overload assignment operator to transfer ownership
66  __attribute__((nothrow)) AthenaArray<T> &operator= (AthenaArray<T> &&t);
67 
68  // public functions to allocate/deallocate memory for 1D-5D data
69  __attribute__((nothrow)) void NewAthenaArray(int nx1);
70  __attribute__((nothrow)) void NewAthenaArray(int nx2, int nx1);
71  __attribute__((nothrow)) void NewAthenaArray(int nx3, int nx2, int nx1);
72  __attribute__((nothrow)) void NewAthenaArray(int nx4, int nx3, int nx2, int nx1);
73  __attribute__((nothrow)) void NewAthenaArray(int nx5, int nx4, int nx3, int nx2,
74  int nx1);
75  __attribute__((nothrow)) void NewAthenaArray(int nx6, int nx5, int nx4, int nx3,
76  int nx2, int nx1);
78 
79  // public function to swap underlying data pointers of two equally-sized arrays
81  void ZeroClear();
82 
83  // functions to get array dimensions
84  int GetDim1() const { return nx1_; }
85  int GetDim2() const { return nx2_; }
86  int GetDim3() const { return nx3_; }
87  int GetDim4() const { return nx4_; }
88  int GetDim5() const { return nx5_; }
89  int GetDim6() const { return nx6_; }
90 
91  // functions to set array dimensions
92  void SetDim1(int nx1) { nx1_ = nx1; }
93  void SetDim2(int nx2) { nx2_ = nx2; }
94  void SetDim3(int nx3) { nx3_ = nx3; }
95  void SetDim4(int nx4) { nx4_ = nx4; }
96  void SetDim5(int nx5) { nx5_ = nx5; }
97  void SetDim6(int nx6) { nx6_ = nx6; }
98 
99  // a function to get the total size of the array
100  int GetSize() const {
101  if (state_ == DataStatus::empty)
102  return 0;
103  else
104  return nx1_*nx2_*nx3_*nx4_*nx5_*nx6_;
105  }
106  std::size_t GetSizeInBytes() const {
107  if (state_ == DataStatus::empty)
108  return 0;
109  else
110  return nx1_*nx2_*nx3_*nx4_*nx5_*nx6_*sizeof(T);
111  }
112 
114  bool IsEmpty() { return (state_ == DataStatus::empty); }
115  bool IsAllocated() { return (state_ == DataStatus::allocated); }
116  // "getter" function to access private data member
117  // TODO(felker): Replace this unrestricted "getter" with a limited, safer alternative.
118  // TODO(felker): Rename function. Conflicts with "AthenaArray<> data" OutputData member.
119  T *data() { return pdata_; }
120  const T *data() const { return pdata_; }
121 
122  // overload "function call" operator() to access 1d-5d data
123  // provides Fortran-like syntax for multidimensional arrays vs. "subscript" operator[]
124 
125  // "non-const variants" called for "AthenaArray<T>()" provide read/write access via
126  // returning by reference, enabling assignment on returned l-value, e.g.: a(3) = 3.0;
127  T &operator() (const int n) {
128  return pdata_[n]; }
129  // "const variants" called for "const AthenaArray<T>" returns T by value, since T is
130  // typically a built-in type (versus "const T &" to avoid copying for general types)
131  T operator() (const int n) const {
132  return pdata_[n]; }
133 
134  T &operator() (const int n, const int i) {
135  return pdata_[i + nx1_*n]; }
136  T operator() (const int n, const int i) const {
137  return pdata_[i + nx1_*n]; }
138 
139  T &operator() (const int n, const int j, const int i) {
140  return pdata_[i + nx1_*(j + nx2_*n)]; }
141  T operator() (const int n, const int j, const int i) const {
142  return pdata_[i + nx1_*(j + nx2_*n)]; }
143 
144  T &operator() (const int n, const int k, const int j, const int i) {
145  return pdata_[i + nx1_*(j + nx2_*(k + nx3_*n))]; }
146  T operator() (const int n, const int k, const int j, const int i) const {
147  return pdata_[i + nx1_*(j + nx2_*(k + nx3_*n))]; }
148 
149  T &operator() (const int m, const int n, const int k, const int j, const int i) {
150  return pdata_[i + nx1_*(j + nx2_*(k + nx3_*(n + nx4_*m)))]; }
151  T operator() (const int m, const int n, const int k, const int j, const int i) const {
152  return pdata_[i + nx1_*(j + nx2_*(k + nx3_*(n + nx4_*m)))]; }
153 
154  // int l?, int o?
155  T &operator() (const int p, const int m, const int n, const int k, const int j,
156  const int i) {
157  return pdata_[i + nx1_*(j + nx2_*(k + nx3_*(n + nx4_*(m + nx5_*p))))]; }
158  T operator() (const int p, const int m, const int n, const int k, const int j,
159  const int i) const {
160  return pdata_[i + nx1_*(j + nx2_*(k + nx3_*(n + nx4_*(m + nx5_*p))))]; }
161 
162  // (deferred) initialize an array with slice from another array
163  void InitWithShallowSlice(AthenaArray<T> &src, const int dim, const int indx,
164  const int nvar);
165 
166  // access stride iterator
167  StrideIterator<T*> at(int i) const {
168  return StrideIterator<T*>(pdata_ + i, nx1_);
169  }
170 
171  StrideIterator<T*> at(int j, int i) const {
172  return StrideIterator<T*>(pdata_ + i + nx1_*j, nx1_*nx2_);
173  }
174 
175  StrideIterator<T*> at(int k, int j, int i) const {
176  return StrideIterator<T*>(pdata_ + i + nx1_*j + nx1_*nx2_*k, nx1_*nx2_*nx3_);
177  }
178 
179  private:
180  T *pdata_;
182  DataStatus state_; // describe what "pdata_" points to and ownership of allocated data
183 
184  void AllocateData();
185 };
186 
187 
188 // destructor
189 
190 template<typename T>
192  DeleteAthenaArray();
193 }
194 
195 // copy constructor (does a deep copy)
196 
197 template<typename T>
199  nx1_ = src.nx1_;
200  nx2_ = src.nx2_;
201  nx3_ = src.nx3_;
202  nx4_ = src.nx4_;
203  nx5_ = src.nx5_;
204  nx6_ = src.nx6_;
205  if (src.pdata_) {
206  std::size_t size = (src.nx1_)*(src.nx2_)*(src.nx3_)*(src.nx4_)*(src.nx5_);
207  pdata_ = new T[size]; // allocate memory for array data
208  for (std::size_t i=0; i<size; ++i) {
209  pdata_[i] = src.pdata_[i]; // copy data (not just addresses!) into new memory
210  }
211  state_ = DataStatus::allocated;
212  }
213 }
214 
215 // copy assignment operator (does a deep copy). Does not allocate memory for destination.
216 // THIS REQUIRES THAT THE DESTINATION ARRAY IS ALREADY ALLOCATED & THE SAME SIZE AS SOURCE
217 
218 template<typename T>
219 __attribute__((nothrow))
221  if (this != &src) {
222  // setting nxN_ is redundant given the above (unenforced) constraint on allowed usage
223  nx1_ = src.nx1_;
224  nx2_ = src.nx2_;
225  nx3_ = src.nx3_;
226  nx4_ = src.nx4_;
227  nx5_ = src.nx5_;
228  nx6_ = src.nx6_;
229  std::size_t size = (src.nx1_)*(src.nx2_)*(src.nx3_)*(src.nx4_)*(src.nx5_)*(src.nx6_);
230  for (std::size_t i=0; i<size; ++i) {
231  this->pdata_[i] = src.pdata_[i]; // copy data (not just addresses!)
232  }
233  state_ = DataStatus::allocated;
234  }
235  return *this;
236 }
237 
238 // move constructor
239 template<typename T>
241  nx1_ = src.nx1_;
242  nx2_ = src.nx2_;
243  nx3_ = src.nx3_;
244  nx4_ = src.nx4_;
245  nx5_ = src.nx5_;
246  nx6_ = src.nx6_;
247  if (src.pdata_) {
248  // && (src.state_ != DataStatus::allocated){ // (if forbidden to move shallow slices)
249  // ---- >state_ = DataStatus::allocated;
250 
251  // Allowing src shallow-sliced AthenaArray to serve as move constructor argument
252  state_ = src.state_;
253  pdata_ = src.pdata_;
254  // remove ownership of data from src to prevent it from free'ing the resources
255  src.pdata_ = nullptr;
256  src.state_ = DataStatus::empty;
257  src.nx1_ = 0;
258  src.nx2_ = 0;
259  src.nx3_ = 0;
260  src.nx4_ = 0;
261  src.nx5_ = 0;
262  src.nx6_ = 0;
263  }
264 }
265 
266 // move assignment operator
267 template<typename T>
268 __attribute__((nothrow))
270  if (this != &src) {
271  // free the target AthenaArray to prepare to receive src pdata_
272  DeleteAthenaArray();
273  if (src.pdata_) {
274  nx1_ = src.nx1_;
275  nx2_ = src.nx2_;
276  nx3_ = src.nx3_;
277  nx4_ = src.nx4_;
278  nx5_ = src.nx5_;
279  nx6_ = src.nx6_;
280  state_ = src.state_;
281  pdata_ = src.pdata_;
282 
283  src.pdata_ = nullptr;
284  src.state_ = DataStatus::empty;
285  src.nx1_ = 0;
286  src.nx2_ = 0;
287  src.nx3_ = 0;
288  src.nx4_ = 0;
289  src.nx5_ = 0;
290  src.nx6_ = 0;
291  }
292  }
293  return *this;
294 }
295 
296 
297 //----------------------------------------------------------------------------------------
299 // \brief shallow copy of nvar elements in dimension dim of an array, starting at
300 // index=indx. Copies pointer to data, but not data itself.
301 
302 // Shallow slice is only able to address the "nvar" range in "dim", and all entries of
303 // the src array for d<dim (cannot access any nx4=2, etc. entries if dim=3 for example)
304 
305 template<typename T>
307  const int indx, const int nvar) {
308  pdata_ = src.pdata_;
309  if (dim == 6) {
310  nx6_ = nvar;
311  nx5_ = src.nx5_;
312  nx4_ = src.nx4_;
313  nx3_ = src.nx3_;
314  nx2_ = src.nx2_;
315  nx1_ = src.nx1_;
316  pdata_ += indx*(nx1_*nx2_*nx3_*nx4_*nx5_);
317  } else if (dim == 5) {
318  nx6_ = 1;
319  nx5_ = nvar;
320  nx4_ = src.nx4_;
321  nx3_ = src.nx3_;
322  nx2_ = src.nx2_;
323  nx1_ = src.nx1_;
324  pdata_ += indx*(nx1_*nx2_*nx3_*nx4_);
325  } else if (dim == 4) {
326  nx6_ = 1;
327  nx5_ = 1;
328  nx4_ = nvar;
329  nx3_ = src.nx3_;
330  nx2_ = src.nx2_;
331  nx1_ = src.nx1_;
332  pdata_ += indx*(nx1_*nx2_*nx3_);
333  } else if (dim == 3) {
334  nx6_ = 1;
335  nx5_ = 1;
336  nx4_ = 1;
337  nx3_ = nvar;
338  nx2_ = src.nx2_;
339  nx1_ = src.nx1_;
340  pdata_ += indx*(nx1_*nx2_);
341  } else if (dim == 2) {
342  nx6_ = 1;
343  nx5_ = 1;
344  nx4_ = 1;
345  nx3_ = 1;
346  nx2_ = nvar;
347  nx1_ = src.nx1_;
348  pdata_ += indx*(nx1_);
349  } else if (dim == 1) {
350  nx6_ = 1;
351  nx5_ = 1;
352  nx4_ = 1;
353  nx3_ = 1;
354  nx2_ = 1;
355  nx1_ = nvar;
356  pdata_ += indx;
357  }
358  state_ = DataStatus::shallow_slice;
359  return;
360 }
361 
362 //----------------------------------------------------------------------------------------
364 // \brief allocate new 1D array with elements initialized to zero.
365 
366 template<typename T>
367 __attribute__((nothrow)) void AthenaArray<T>::NewAthenaArray(int nx1) {
368  state_ = DataStatus::allocated;
369  nx1_ = nx1;
370  nx2_ = 1;
371  nx3_ = 1;
372  nx4_ = 1;
373  nx5_ = 1;
374  nx6_ = 1;
375  pdata_ = new T[nx1](); // allocate memory and initialize to zero
376 }
377 
378 //----------------------------------------------------------------------------------------
380 // \brief 2d data allocation
381 
382 template<typename T>
383 __attribute__((nothrow)) void AthenaArray<T>::NewAthenaArray(int nx2, int nx1) {
384  state_ = DataStatus::allocated;
385  nx1_ = nx1;
386  nx2_ = nx2;
387  nx3_ = 1;
388  nx4_ = 1;
389  nx5_ = 1;
390  nx6_ = 1;
391  pdata_ = new T[nx1*nx2](); // allocate memory and initialize to zero
392 }
393 
394 //----------------------------------------------------------------------------------------
396 // \brief 3d data allocation
397 
398 template<typename T>
399 __attribute__((nothrow)) void AthenaArray<T>::NewAthenaArray(int nx3, int nx2, int nx1) {
400  state_ = DataStatus::allocated;
401  nx1_ = nx1;
402  nx2_ = nx2;
403  nx3_ = nx3;
404  nx4_ = 1;
405  nx5_ = 1;
406  nx6_ = 1;
407  pdata_ = new T[nx1*nx2*nx3](); // allocate memory and initialize to zero
408 }
409 
410 //----------------------------------------------------------------------------------------
412 // \brief 4d data allocation
413 
414 template<typename T>
415 __attribute__((nothrow)) void AthenaArray<T>::NewAthenaArray(int nx4, int nx3, int nx2,
416  int nx1) {
417  state_ = DataStatus::allocated;
418  nx1_ = nx1;
419  nx2_ = nx2;
420  nx3_ = nx3;
421  nx4_ = nx4;
422  nx5_ = 1;
423  nx6_ = 1;
424  pdata_ = new T[nx1*nx2*nx3*nx4](); // allocate memory and initialize to zero
425 }
426 
427 //----------------------------------------------------------------------------------------
429 // \brief 5d data allocation
430 
431 template<typename T>
432 __attribute__((nothrow)) void AthenaArray<T>::NewAthenaArray(int nx5, int nx4, int nx3,
433  int nx2, int nx1) {
434  state_ = DataStatus::allocated;
435  nx1_ = nx1;
436  nx2_ = nx2;
437  nx3_ = nx3;
438  nx4_ = nx4;
439  nx5_ = nx5;
440  nx6_ = 1;
441  pdata_ = new T[nx1*nx2*nx3*nx4*nx5](); // allocate memory and initialize to zero
442 }
443 
444 //----------------------------------------------------------------------------------------
446 // \brief 6d data allocation
447 
448 template<typename T>
449 __attribute__((nothrow)) void AthenaArray<T>::NewAthenaArray(int nx6, int nx5, int nx4,
450  int nx3, int nx2, int nx1) {
451  state_ = DataStatus::allocated;
452  nx1_ = nx1;
453  nx2_ = nx2;
454  nx3_ = nx3;
455  nx4_ = nx4;
456  nx5_ = nx5;
457  nx6_ = nx6;
458  pdata_ = new T[nx1*nx2*nx3*nx4*nx5*nx6](); // allocate memory and initialize to zero
459 }
460 
461 //----------------------------------------------------------------------------------------
463 // \brief free memory allocated for data array
464 
465 template<typename T>
467  // state_ is tracked partly for correctness of delete[] operation in DeleteAthenaArray()
468  switch (state_) {
469  case DataStatus::empty:
470  case DataStatus::shallow_slice:
471  pdata_ = nullptr;
472  break;
473  case DataStatus::allocated:
474  delete[] pdata_;
475  pdata_ = nullptr;
476  state_ = DataStatus::empty;
477  break;
478  }
479 }
480 
481 //----------------------------------------------------------------------------------------
483 // \brief swap pdata_ pointers of two equally sized AthenaArrays (shallow swap)
484 // Does not allocate memory for either AthenArray
485 // THIS REQUIRES THAT THE DESTINATION AND SOURCE ARRAYS BE ALREADY ALLOCATED (state_ !=
486 // empty) AND HAVE THE SAME SIZES (does not explicitly check either condition)
487 
488 template<typename T>
490  std::swap(pdata_, array2.pdata_);
491  return;
492 }
493 
494 //----------------------------------------------------------------------------------------
496 // \brief fill the array with zero
497 
498 template<typename T>
500  switch (state_) {
501  case DataStatus::empty:
502  break;
503  case DataStatus::shallow_slice:
504  case DataStatus::allocated:
505  // allocate memory and initialize to zero
506  std::memset(pdata_, 0, GetSizeInBytes());
507  break;
508  }
509 }
510 
511 //----------------------------------------------------------------------------------------
513 // \brief to be called in non-default ctors, if immediate memory allocation is requested
514 // (could replace all "new[]" calls in NewAthenaArray function overloads)
515 
516 template<typename T>
518  switch (state_) {
519  case DataStatus::empty:
520  case DataStatus::shallow_slice: // init=shallow_slice should never be passed to ctor
521  break;
522  case DataStatus::allocated:
523  // allocate memory and initialize to zero
524  pdata_ = new T[nx1_*nx2_*nx3_*nx4_*nx5_*nx6_]();
525  break;
526  }
527 }
528 
529 #endif // ATHENA_ARRAYS_HPP_
__attribute__((nothrow)) AthenaArray< T >
int GetDim3() const
__attribute__((nothrow)) AthenaArray(AthenaArray< T > &&t)
int GetDim2() const
bool IsAllocated()
std::size_t GetSizeInBytes() const
DataStatus state_
__attribute__((nothrow)) void NewAthenaArray(int nx3
__attribute__((nothrow)) void NewAthenaArray(int nx4
__attribute__((nothrow)) void NewAthenaArray(int nx6
void InitWithShallowSlice(AthenaArray< T > &src, const int dim, const int indx, const int nvar)
void SetDim6(int nx6)
int GetDim4() const
StrideIterator< T * > at(int j, int i) const
T & operator()(const int n)
StrideIterator< T * > at(int i) const
StrideIterator< T * > at(int k, int j, int i) const
const T * data() const
__attribute__((nothrow)) AthenaArray< T > &operator
__attribute__((nothrow)) AthenaArray(const AthenaArray< T > &t)
AthenaArray(int nx5, int nx4, int nx3, int nx2, int nx1, DataStatus init=DataStatus::allocated)
AthenaArray(int nx2, int nx1, DataStatus init=DataStatus::allocated)
int GetDim6() const
__attribute__((nothrow)) void NewAthenaArray(int nx2
int GetSize() const
void SwapAthenaArray(AthenaArray< T > &array2)
void SetDim1(int nx1)
void DeleteAthenaArray()
void SetDim2(int nx2)
void SetDim3(int nx3)
AthenaArray(int nx6, int nx5, int nx4, int nx3, int nx2, int nx1, DataStatus init=DataStatus::allocated)
void AllocateData()
AthenaArray(int nx1, DataStatus init=DataStatus::allocated)
void SetDim4(int nx4)
int GetDim1() const
int GetDim5() const
bool IsShallowSlice()
void SetDim5(int nx5)
__attribute__((nothrow)) void NewAthenaArray(int nx1)
AthenaArray(int nx3, int nx2, int nx1, DataStatus init=DataStatus::allocated)
AthenaArray(int nx4, int nx3, int nx2, int nx1, DataStatus init=DataStatus::allocated)
__attribute__((nothrow)) void NewAthenaArray(int nx5