Learn » Examples » Example Integral Reduction Program

The full code for this example can be found here: examples/collection/reduce_integral.cc.

//
// This code computes a composite trapezoidal rule for the integral
// \int_{0}^{1} f(x) dx
//
// The function 'f' is defined in the code.
//
// The interval [0, 1] is broken into uniform sub-intervals of length
//
// H = 1.0 / (# of objects)
//
// Each sub-interval [xi, xi + H] is broken into smaller parts of length
//
// h = H / (# of parts per object)
//


#include <vt/transport.h>
#include <cmath>

static constexpr std::size_t const default_nparts_object = 8;
static constexpr std::size_t const default_num_objs = 4;
static constexpr std::size_t const verbose = 1;
static constexpr vt::NodeType const reduce_root_node = 0;

static bool root_reduce_finished = false;
static double exactIntegral = 0.0;


//
// Function 'f' to integrate over [0, 1]
//


double f(double x) {
  //----
  //exactIntegral = 1.0 / 3.0;
  //return x*x;
  //----
  exactIntegral = M_2_PI;
  return sin(M_PI * x);
}

using ReduceMsg = vt::collective::ReduceTMsg<double>;

struct Integration1D : vt::Collection<Integration1D, vt::Index1D> {

private:

  size_t numObjs_ = default_num_objs;
  size_t numPartsPerObject_ = default_nparts_object;

public:

  explicit Integration1D()
    : numObjs_(default_num_objs),
      numPartsPerObject_(default_nparts_object)
  { }

  void checkIntegral(double val) {
    fmt::print(" >> The integral over [0, 1] is {}\n", val);
    fmt::print(
      " >> The absolute error is {}\n",
      std::fabs(val - exactIntegral)
    );
    //
    // Set the 'root_reduce_finished' variable to true.
    //
    root_reduce_finished = true;
  }

  struct InitMsg : vt::CollectionMessage<Integration1D> {

    size_t numObjects = 0;
    size_t nIntervalPerObject = 0;

    InitMsg(const size_t nobjs, const size_t nint)
      : numObjects(nobjs), nIntervalPerObject(nint)
    { }

  };

  void compute(InitMsg *msg) {

    numObjs_ = msg->numObjects;
    numPartsPerObject_ = msg->nIntervalPerObject;

    //
    // Compute the integral with the trapezoidal rule
    // over the interval [a, b]
    //
    // a = numPartsPerObject_ * (IndexID) * h
    //
    // b = a + numPartsPerObject_ * h
    //
    // where h = 1.0 / (numPartsPerObject_ * numObjs_ )
    //
    // Since 0 <= (IndexID) < numObjs_, the union of these intervals
    // covers [0, 1]
    //

    double h = 1.0 / (numPartsPerObject_ * numObjs_ );
    double quadsum = 0.0;

    double a = numPartsPerObject_ * getIndex().x() * h;

    //
    // Apply composite trapezoidal rule over
    //
    // [a, a + numPartsPerObject_ * h]
    //
    for (size_t ii = 0; ii < numPartsPerObject_; ++ii) {
      double x0 = a + ii * h;
      /* --- Trapeze Quadrature Rule over the interval [x0, x0+h] */
      quadsum += 0.5 * h * ( f(x0) + f(x0+h) );
    }

    if (verbose > 0) {
      double b = a + numPartsPerObject_ * h;
      fmt::print(
        " Interval [{}, {}], on node {} & object {}, "
        "has integral {}.\n", a, b, vt::theContext()->getNode(),
        getIndex(), quadsum
      );
    }

    //
    // Reduce the partial sum to get the integral over [0, 1]
    //

    auto proxy = this->getCollectionProxy();
    proxy.reduce<&Integration1D::checkIntegral, vt::collective::PlusOp>(
      proxy[0], quadsum
    );
  }

};


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

  size_t num_objs = default_num_objs;
  size_t numIntPerObject = default_nparts_object;

  std::string name(argv[0]);

  vt::initialize(argc, argv);

  vt::NodeType this_node = vt::theContext()->getNode();
  vt::NodeType num_nodes = vt::theContext()->getNumNodes();

  if (argc == 1) {
    if (this_node == 0) {
      fmt::print(
        stderr, "{}: using default arguments since none provided\n", name
      );
    }
    num_objs = default_num_objs * num_nodes;
  } else {
    if (argc == 2) {
      num_objs = (size_t) strtol(argv[1], nullptr, 10);
    }
    else if (argc == 3) {
      num_objs = (size_t) strtol(argv[1], nullptr, 10);
      numIntPerObject = (size_t) strtol(argv[2], nullptr, 10);
    }
    else {
      fmt::print(
        stderr,
        "usage: {} <num-objects> <num-interval-per-object>\n", name
      );
      return 1;
    }
  }

  //
  // Create the interval decomposition into objects
  //
  using BaseIndexType = typename vt::Index1D::DenseIndexType;
  auto range = vt::Index1D(static_cast<BaseIndexType>(num_objs));

  auto proxy = vt::makeCollection<Integration1D>("examples_reduce_integral")
    .bounds(range)
    .bulkInsert()
    .wait();

  vt::runInEpochCollective([proxy, num_objs, numIntPerObject]{
    proxy.broadcastCollective<Integration1D::InitMsg, &Integration1D::compute>(
      num_objs, numIntPerObject
    );
  });

  // Add something like this to validate the reduction.
  // Create the variable root_reduce_finished as a static variable,
  // which is only checked on one node.
  if (vt::theContext()->getNode() == reduce_root_node) {
    vtAssertExpr(root_reduce_finished == true);
  }

  vt::finalize();

  return 0;

}