functor.cpp

/*
 * All the data members of ElementData are public, allowing the
 * argument class direct access. The auxiliary classes are functors.
 */
#include <iostream>
#include <vector>
#include <string>

// Forward definitions
class Level1;
class Level2;

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

class Level1 {
public:
  Level1(ElementData& data) : x(data.x), y(data.y), z(data.z), dx_dt(0.) { };
  // Only a limited set of variables from ElementData
  double& x;
  double& y;
  double& z;

  double dx_dt;
  
  void operator()();
};

class Level2 {
public:
  Level2(Level1& data) : x(data.x), y(data.y), dx_dt(data.dx_dt) { };
  // A limited and const qualified set of references from level 1
  const double& x;
  const double& y;
  double& dx_dt;

  void operator()();
};

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

  const double dt = 0.1;
  
  double dy_dt;
  double dz_dt;

  Level2 l2(*this);
  l2();

  dy_dt = x * (rho - z) - y;
  dz_dt = x * y - beta * z;

  x += dx_dt * dt;
  y += dy_dt * dt;
  z += dz_dt * dt;
}  
  
void Level2::operator()() {
  const double sigma = 10.;
  
  dx_dt = sigma * (y - x);
}

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

  Level1 l1(e);
  l1();

  std::cout << "Scalar run" << std::endl;
  std::cout << "x = " << e.x << ", y = " << e.y << ", z = " << e.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++) {
        Level1 l1(*iter);
	l1();
    }
  }

  // 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->x;
    sumy += iter->y;
    sumz += iter->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;
}