data_access/infinitely_elementary.cpp at main · timspainNERSC/data_access · GitHub

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
/*
 * Every function is a member of ElementData and can read-write any
 * data member
 */
#include <iostream>
#include <vector>
#include <string>

class ElementData {
public:
  ElementData() : m_w(0.), m_x(0.), m_y(0.), m_z(0.) { };
  ElementData(double w, double x, double y, double z) : m_w(w), m_x(x), m_y(y), m_z(z) { };

  void level1Function();
  void level2Function();

  double m_w;
  double m_x;
  double m_y;
  double m_z;

private:
  double m_dx_dt;

};

void ElementData::level1Function() {
  const double beta = 8./3.;
  const double rho = 28.;

  const double dt = 0.1;

  double dy_dt;
  double dz_dt;

  dy_dt = m_x * (rho - m_z) - m_y;
  dz_dt = m_x * m_y - beta * m_z;

  m_x += m_dx_dt * dt;
  m_y += dy_dt * dt;
  m_z += dz_dt * dt;
}

void ElementData::level2Function() {
  const double sigma = 10.;

  m_dx_dt = sigma * (m_y - m_x);
}

void scalar( ) {
  ElementData e(2., 1., 1., 1.);
  std::cout << "x = " << e.m_x << ", y = " << e.m_y << ", z = " << e.m_z << std::endl;

  e.level1Function();

  std::cout << "Scalar run" << std::endl;
  std::cout << "x = " << e.m_x << ", y = " << e.m_y << ", z = " << e.m_z << std::endl;
}

void iterated_array(int nxyz, int timesteps) {
  double lorenz_range = 10.;
  double lorenz_start = -lorenz_range/2;
  double lorenz_step = lorenz_range/nxyz;

  std::vector<ElementData> modelData = std::vector<ElementData>(nxyz*nxyz*nxyz);

  // Fill the arrays
  for (int i = 0; i < nxyz; i++) {
    for (int j = 0; j < nxyz; j++) {
      for (int k = 0; k < nxyz; k++) {
	modelData[k + nxyz * (j + nxyz * i)] =
	  ElementData(
	   i * 10000. + j * 100. + k,
	   -lorenz_start + lorenz_step * i,
	   -lorenz_start + lorenz_step * j,
	   -lorenz_start + lorenz_step * k
	   );
      }
    }
  }

  // Loop over time and all elements
  for (int t = 0; t < timesteps; t++) {
    for (auto iter = modelData.begin(); iter != modelData.end(); iter++) {
	iter->level1Function();
    }
  }

  // Sum all the elements, so that the compiler cannot omptimize away the main loop as not doing anything
  double sumx = 0.;
  double sumy = 0.;
  double sumz = 0.;
  for (auto iter = modelData.begin(); iter != modelData.end(); iter++) {
    sumx += iter->m_x;
    sumy += iter->m_y;
    sumz += iter->m_z;
  }

  std::cout << "Variable totals: x = " << sumx << ", y = " << sumy << ", z = " << sumz << std::endl;

}

int main(int argc, char* argv[]) {

  scalar();

  int nxyz;
  int ntimesteps;
  if (argc > 1) {
    nxyz = std::stoi(std::string(argv[1]));
  } else {
    nxyz = 100;
  }

  if (argc > 2) {
    ntimesteps = std::stoi(std::string(argv[2]));
  } else {
    ntimesteps = 100;
  }

  iterated_array(nxyz, ntimesteps);

  return 0;
}