Describe the diffusion problem

project/src/Field.h

 // -*-c++-*-
 
+#pragma once
+
 #include <ostream>
 
 class RowView {
@@ -13,6 +15,10 @@ public:
     return row[i];
   }
 
+  double operator[](unsigned i) const {
+    return row[i];
+  }
+
   double* raw() {
     return row;
   }
@@ -41,6 +47,10 @@ public:
     return RowView(field + i * n);
   }
 
+  RowView operator[](unsigned i) const {
+    return RowView(field + i * n);
+  }
+
   RowView firstRow() {
     return operator[](1);
   }

project/src/Main.cpp

+#include <mpi.h>
+#include "Field.h"
 
-int main(int argc, char **argv) {
+void addGates(Field &f, int nGates, double gateSize) {
 
+}
 
+int main(int argc, char **argv) {
+  MPI::Init();
 
+  MPI::Finalize();
   return 0;
 }

project/src/Makefile

@@ -4,12 +4,16 @@ CPP_SOURCE_FILES = $(filter-out $(ENTRY_POINTS), $(wildcard *.cpp))
 C_SOURCE_FILES = $(wildcard *.c)
 HEADER_FILES = $(wildcard *.h)
 
+MPIFLAGS = $(shell mpic++ -showme:compile)
+MPI_LD_FLAGS = $(shell mpic++ -showme:link)
+
 CXX=g++
 MPICXX=mpic++
-CXX_FLAGS=-O3
+CXX_FLAGS=-O3 -I/usr/lib/x86_64-linux-gnu/openmpi/include/
+
 
 main : Main.cpp $(CPP_SOURCE_FILES) $(HEADER_FILES)
-	$(MPICXX) -O3 Main.cpp $(CPP_SOURCE_FILES) -o main
+	$(CXX) $(MPIFLAGS) -O3 Main.cpp $(CPP_SOURCE_FILES) $(MPI_LD_FLAGS) -o main
 
 test : Test.cpp $(CPP_SOURCE_FILES) $(HEADER_FILES)
 	$(CXX) -g -O0 Test.cpp $(CPP_SOURCE_FILES) -o test

project/src/Problem.cpp

+#include "Problem.h"
+
+void Problem::init() {
+  for (int i = 0; i < ydim + 2; ++i) {
+    for (int j = 0; j < xdim + 2; ++j) {
+      f0[i][j] = 0;
+    }
+  }
+}
+
+void Problem::addGates(unsigned xGates, unsigned yGates, double value) {
+  const double GATE_WIDTH = 0.1;
+
+  if (iField == 0) {
+    // add the xgates
+    for (unsigned i = 0; i < xGates; ++i) {
+      unsigned gateCentre = (1 + 2*i) * xdim / (2 * xGates);
+      unsigned gateStart = gateCentre - 0.5 * GATE_WIDTH * xdim;
+      unsigned gateEnd = gateCentre + 0.5 * GATE_WIDTH * xdim;
+
+      auto row =  f0.firstHalo();
+
+      for (unsigned j = gateStart; j < gateEnd; ++j) {
+        row[j] = value;
+      }
+    }
+  }
+
+  unsigned starty = ydim * iField;
+  unsigned endy = starty + ydim;
+
+  // add the ygates
+  for (unsigned i = 0; i < yGates; ++i) {
+    unsigned gateCentre = (1 + 2*i) * ydim / (2 * yGates);
+    unsigned gateStart = gateCentre - 0.5 * GATE_WIDTH * ydim;
+    unsigned gateEnd = gateCentre + 0.5 * GATE_WIDTH * ydim;
+
+    if (gateStart >= starty && gateStart < endy) {
+      for (unsigned j = gateStart; j < std::min(gateEnd, endy); j++) {
+        f0[j - starty + 1][xdim] = value;
+      }
+    }
+  }
+}
+
+void Problem::step(double cx, double cy) {
+  auto &nextf = *this->nextf;
+  auto &f = *this->currentf;
+
+  for (unsigned y = 1; y < ydim + 1; y++) {
+    for (unsigned x = 1; x < xdim + 1; x++) {
+      nextf[y][x] = 0.25 * ((1.0 + 0.5 * cx) * f[y][x+1] + (1.0 - 0.5 * cx) * f[y][x-1]
+                            + (1.0 + 0.5 * cy) * f[y+1][x] + (1.0 - 0.5 * cy) * f[y-1][x]);
+
+    }
+  }
+
+  std::swap(this->nextf, this->currentf);
+}
+
+double Problem::residual(double cx, double cy) const {
+  double res_sqr = 0.0;
+  const auto &solution = *this->currentf;
+
+  for(unsigned iy = 1 ; iy < ydim + 1; iy++) {
+    for(unsigned ix = 1; ix < xdim + 1; ix++) {
+        double res_i = (1.0+0.5*cy) * solution[iy+1][ix]
+          + (1.0-0.5*cy) * solution[iy-1][ix]
+          + (1.0+0.5*cx) * solution[iy][ix+1]
+          + (1.0-0.5*cx) * solution[iy][ix-1]
+          - 4.0 * solution[iy][ix];
+
+        res_sqr += res_i*res_i;
+    }
+  }
+
+  return res_sqr;
+}

project/src/Problem.h

+#pragma once
+
+#include "Field.h"
+
+class Problem {
+  unsigned nFields;
+  unsigned iField;
+
+  unsigned xdim;
+  unsigned ydim;
+
+  Field f0;
+  Field f1;
+
+  Field *currentf;
+  Field *nextf;
+
+  void init();
+  void addGates(unsigned xgates, unsigned ygates, double value);
+
+public:
+  Problem(unsigned nFields, unsigned iField, unsigned xdim, unsigned ydim,
+          unsigned xgates, unsigned ygates, double gateVal) :
+    nFields(nFields),
+    iField(iField),
+    xdim(xdim),
+    ydim(ydim),
+    f0(xdim, ydim),
+    f1(xdim, ydim),
+    currentf(&f0),
+    nextf(&f1) {
+
+    init();
+    addGates(xgates, ygates, gateVal);
+  }
+
+  void step(double cx, double cy);
+
+  double residual(double cx, double cy) const;
+
+};