Implementations of different problems

include/IsingModel.hpp

@@ -12,7 +12,10 @@
 #ifndef __ISING_MODEL__
 #define __ISING_MODEL__
 
+#include "constants.hpp"
+#include "data_type.hpp"
 #include "typedefs.hpp"
+#include "utils.hpp"
 
 #include <armadillo>
 #include <random>
@@ -27,7 +30,6 @@
 #define DOWN 1
 #define RIGHT 1
 
-
 /** @brief The Ising model in 2 dimensions.
  *
  *  @details None of the methods are parallelized, as there is very little
@@ -36,15 +38,15 @@
 class IsingModel {
 private:
     friend class IsingModelTest;
-    /** @brief \f$ L \cross L \f$ matrix where element $ x \in {-1, 1}$.
+    /** @brief \f$ L \times L \f$ matrix where element \f$ x \in {-1, 1}\f$.
      * */
     arma::Mat<int> lattice;
 
-    /** @brief \f$ L \cross 2 \f$ matrix with the neighbors of each element 
+    /** @brief \f$ L \times 2 \f$ matrix with the neighbors of each element
      *  \f$ x_i \f$.
      *
-     *  @details The reason why it's \f$ L \cross 2 \f$ instead of 
+     *  @details The reason why it's \f$ L \times 2 \f$ instead of
-     *  \f$ L \cross 2 \f$, is that we can see that we can use the same column
+     *  \f$ L \times 2 \f$, is that we can see that we can use the same column
      *  for the left and upper neighbor, and we can use the same column for the
      *  right and lower neighbor.
      * */
@@ -64,19 +66,11 @@ private:
 
     /** @brief The current energy state. unit: \f$ J \f$.
      * */
-    int E;
+    double E;
 
     /** @brief The current magnetic strength. unit: Unitless.
      * */
-    int M;
+    double M;
-
-    /** @brief Energy per spin. unit: \f$ J \f$.
-     * */
-    double eps;
-
-    /** @brief Magnetization per spin. unit: Unitless.
-     * */
-    double m;
 
     /** @brief Initialize the lattice with a random distribution of 1s and
      *  -1s.
@@ -121,19 +115,19 @@ public:
 
     /** @brief The Metropolis algorithm.
      * */
-    void Metropolis(std::mt19937 engine);
+    data_t Metropolis(std::mt19937 &engine);
 
     /** @brief Get the current energy.
      *
-     *  @return int
+     *  @return double
      * */
-    int get_E();
+    double get_E();
 
     /** @brief Get the current magnetization.
      *
-     *  @return int
+     *  @return double
      * */
-    int get_M();
+    double get_M();
 };
 
 #endif

include/data_type.hpp

@@ -0,0 +1,78 @@
+/** @file data_type.hpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 1.0
+ *
+ *  @brief Header for the data_t type.
+ *
+ *  @bug No known bugs
+ * */
+#ifndef __DATA_TYPE__
+#define __DATA_TYPE__
+
+#include <sys/types.h>
+#include <type_traits>
+
+/** @brief Data structure that contains the data needed for the project*/
+struct data_t {
+    double E = 0.; ///< The expected energy
+    double M = 0.; ///< The expected magnetization
+    double E2 = 0.; ///< The expected variance of the energy
+    double M2 = 0.; ///< The expected variance of magnetization
+    double M_abs = 0.; ///< The expected absolute magnetization
+};
+
+/** @brief Define dividing data_t by a number
+ *
+ *  @param data The data to divide
+ *  @param num The number to divide data by
+ *
+ *  @return data_t
+ * */
+template <class T>
+data_t operator/(const data_t &data, T num);
+
+// Explicit instantiation
+extern template data_t operator/(const data_t &, uint);
+extern template data_t operator/(const data_t &, ulong);
+extern template data_t operator/(const data_t &,int);
+extern template data_t operator/(const data_t &,double);
+
+
+/** @brief Define /= on data_t by a number.
+ *
+ *  @param data The data to divide
+ *  @param num The number to divide data by
+ *
+ *  @return data_t
+ * */
+template <class T>
+data_t& operator/=(data_t &data, T num);
+
+// Explicit instantiation
+extern template data_t& operator/=(data_t &, uint);
+extern template data_t& operator/=(data_t &, ulong);
+extern template data_t& operator/=(data_t &,int);
+extern template data_t& operator/=(data_t &,double);
+
+/** @brief Define + on data_t by a data_t.
+ *
+ *  @param a The left side
+ *  @param b The right side
+ *
+ *  @return data_t
+ * */
+data_t operator+(const data_t &a, const data_t &b);
+
+/** @brief Define += on data_t by a data_t.
+ *
+ *  @param a The left side
+ *  @param b The right side
+ *
+ *  @return data_t
+ * */
+data_t& operator+=(data_t &a, const data_t &b);
+
+#endif

include/monte_carlo.hpp

@@ -12,8 +12,88 @@
 #ifndef __MONTE_CARLO__
 #define __MONTE_CARLO__
 
-void burn_in_time();
+#include "IsingModel.hpp"
+#include "data_type.hpp"
+#include "utils.hpp"
 
-void pd_estimate();
+#include <functional>
+#include <string>
+
+#define BURN_IN_TIME 1000
+
+#define EPS_2 (-2 * std::sinh(8.)) / (std::cosh(8.) + 3)
+
+#define MAG_2 (std::exp(8.) + 1) / (2 * cosh(8.) + 3)
+
+#define CV_2                                                                   \
+    16                                                                         \
+        * (3 * std::cosh(8.) + std::cosh(8.) * std::cosh(8.)                   \
+           - std::sinh(8.) * std::sinh(8.))                                    \
+        / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))
+
+#define X_2                                                                    \
+    (3 * std::exp(8.) + std::exp(-8.) + 3)                                     \
+        / ((std::cosh(8.) + 3) * (std::cosh(8.) + 3))
+
+/** @brief Test numerical data with analytical data.
+ *
+ *  @param tol The tolerance between the analytical and numerical solution.
+ *  @param max_cycles The max number of Monte Carlo cycles.
+ *
+ *  return uint
+ * */
+uint test_2x2_lattice(double tol, uint max_cycles);
+
+/** @brief Write the expected values for each Monte Carlo cycles to file.
+ *
+ *  @param T Temperature
+ *  @param L The size of the lattice
+ *  @param cycles The amount of Monte Carlo cycles to do
+ *  @param filename The file to write to
+ * */
+void monte_carlo_progression(double T, uint L, uint cycles,
+                             const std::string filename);
+
+/** @brief Estimate the probability distribution for the energy.
+ *
+ *  @param T The temperature of the Ising model
+ *  @param L The size of the lattice
+ *  @param cycles The amount of Monte Carlo cycles to do
+ *  @param filename The file to write to
+ * */
+void pd_estimate(double T, uint L, uint cycles, const std::string filename);
+
+/** @brief Execute the Metropolis algorithm for a certain amount of Monte 
+ *  Carlo cycles.
+ *
+ *  @param data The data to store the results
+ *  @param L The size of the lattice
+ *  @param T The Temperature for the Ising model
+ *  @param cycles The amount of Monte Carlo cycles to do*/
+void monte_carlo_serial(data_t &data, uint L, double T, uint cycles);
+
+/** @brief Execute the Metropolis algorithm for a certain amount of Monte 
+ *  Carlo cycles in parallel.
+ *
+ *  @param data The data to store the results
+ *  @param L The size of the lattice
+ *  @param T The Temperature for the Ising model
+ *  @param cycles The amount of Monte Carlo cycles to do
+ * */
+void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles);
+
+/** @brief Perform the MCMC algorithm using a range of temperatures.
+ *
+ *  @param L The size of the lattice
+ *  @param start_T The start temperature
+ *  @param end_T The end temperature
+ *  @param point_T The amount of point to measure
+ *  @param monte_carlo Which Monte Carlo implementation to use
+ *  @param outfile The file to write the data to
+ * */
+void phase_transition(
+    uint L, double start_T, double end_T, uint points_T,
+    std::function<void(data_t &, uint, double, uint)> monte_carlo,
+    std::string outfile);
 
 #endif

include/testlib.hpp

@@ -5,9 +5,10 @@
  *
  *  @version 1.0
  *
- *  @brief Function prototypes and macros that are useful.
+ *  @brief A small test library.
  *
- *  @details This a small testing library that is tailored for the needs of the project.
+ *  @details This a small testing library that is tailored for the needs of the
+ * project.
  *
  *  @bug No known bugs
  * */
@@ -27,8 +28,9 @@
  *  assertion function than the regular assert function from cassert.
  * */
 #define ASSERT(expr, msg)                                                      \
-    m_assert(expr, #expr, __METHOD_NAME__, __FILE__, __LINE__, msg)
+    details::m_assert(expr, #expr, __METHOD_NAME__, __FILE__, __LINE__, msg)
 
+namespace details {
 /** @brief Test an expression, confirm that test is ok, or abort execution.
  *
  *  @details This function takes in an expression and prints an OK message if
@@ -43,7 +45,9 @@
  * */
 void m_assert(bool expr, std::string expr_str, std::string func,
               std::string file, int line, std::string msg);
+} // namespace details
 
+namespace testlib {
 /** @brief Test if two armadillo matrices/vectors are close to each other.
  *
  *  @details This function takes in 2 matrices/vectors and checks if they are
@@ -64,12 +68,29 @@ static bool close_to(arma::Mat<T> &a, arma::Mat<T> &b, double tol = 1e-8)
     }
 
     for (size_t i = 0; i < a.n_elem; i++) {
-        if (std::abs(a(i) - b(i)) >= tol) {
+        if (!close_to(a(i), b(i))) {
             return false;
         }
     }
     return true;
+}
 
+/** @brief Test if two numbers are close to each other.
+ *
+ *  @details This function takes in 2 matrices/vectors and checks if they are
+ *  approximately equal to each other given a tolerance.
+ *
+ *  @param a Matrix/vector a
+ *  @param b Matrix/vector b
+ *  @param tol The tolerance
+ *
+ *  @return bool
+ * */
+template <class T,
+          class = typename std::enable_if<std::is_arithmetic<T>::value>::type>
+static bool close_to(T a, T b, double tol = 1e-8)
+{
+    return std::abs(a - b) < tol;
 }
 
 /** @brief Test if two armadillo matrices/vectors are equal.
@@ -93,6 +114,7 @@ static bool is_equal(arma::Mat<T> &a, arma::Mat<T> &b)
     }
     return true;
 }
+
 /** @brief Test that all elements fulfill the condition.
  *
  *  @param expr The boolean expression to apply to each element
@@ -111,5 +133,5 @@ static bool assert_each(std::function<bool(T)> expr, arma::Mat<T> &M)
     }
     return true;
 }
-
+} // namespace testlib
 #endif

include/utils.hpp

@@ -40,7 +40,31 @@
 /** @def __METHOD_NAME__
  *  @brief Get the name of the current method/function without the return type.
  * */
-#define __METHOD_NAME__ methodName(__PRETTY_FUNCTION__)
+#define __METHOD_NAME__ details::methodName(__PRETTY_FUNCTION__)
+
+namespace details {
+/** @brief Takes in the __PRETTY_FUNCTION__ string and removes the return type.
+ *
+ *  @details This function should only be used for the __METHOD_NAME__ macro,
+ *  since it takes the output from __PRETTY_FUNCTION__ and strips the return
+ *  type.
+ *
+ *  @param pretty_function The string from __PRETTY_FUNCTION__
+ *
+ *  @return std::string
+ * */
+inline std::string methodName(const std::string &pretty_function)
+{
+    size_t colons = pretty_function.find("::");
+    size_t begin = pretty_function.substr(0, colons).rfind(" ") + 1;
+    size_t end = pretty_function.rfind("(") - begin;
+
+    return pretty_function.substr(begin, end) + "()";
+}
+
+} // namespace details
+
+namespace utils {
 
 /** @brief Turns a double into a string written in scientific format.
  *
@@ -68,26 +92,6 @@ std::string scientific_format(double d, int width = 20, int prec = 10);
 std::string scientific_format(const std::vector<double> &v, int width = 20,
                               int prec = 10);
 
-
-/** @brief Takes in the __PRETTY_FUNCTION__ string and removes the return type.
- *
- *  @details This function should only be used for the __METHOD_NAME__ macro,
- *  since it takes the output from __PRETTY_FUNCTION__ and strips the return
- *  type.
- *
- *  @param pretty_function The string from __PRETTY_FUNCTION__
- *
- *  @return std::string
- * */
-static inline std::string methodName(const std::string &pretty_function)
-{
-    size_t colons = pretty_function.find("::");
-    size_t begin = pretty_function.substr(0, colons).rfind(" ") + 1;
-    size_t end = pretty_function.rfind("(") - begin;
-
-    return pretty_function.substr(begin, end) + "()";
-}
-
 /** @brief Make path given.
  *
  *  @details This tries to be the equivalent to "mkdir -p" and creates a new
@@ -100,4 +104,14 @@ static inline std::string methodName(const std::string &pretty_function)
  * */
 bool mkpath(std::string path, int mode = 0777);
 
+/** @brief Get the directory name of the path
+ *
+ *  @param path The path to use.
+ *
+ *  @return string
+ * */
+std::string dirname(const std::string &path);
+
+} // namespace utils
+
 #endif

src/IsingModel.cpp

@@ -9,14 +9,8 @@
  *
  *  @bug No known bugs
  * */
-
 #include "IsingModel.hpp"
 
-#include "constants.hpp"
-
-#include <cmath>
-#include <random>
-
 IsingModel::IsingModel()
 {
 }
@@ -47,7 +41,7 @@ IsingModel::IsingModel(uint L, double T, int val)
 void IsingModel::initialize_lattice()
 {
     this->lattice.set_size(this->L, this->L);
-    std::random_device rd{};
+    std::random_device rd;
     std::mt19937 engine(rd());
 
     std::uniform_int_distribution<int> coin_flip(0, 1);
@@ -69,36 +63,41 @@ void IsingModel::initialize_neighbors()
 
 void IsingModel::initialize_energy_diff()
 {
-    for (size_t i = -8; i <= 8; i += 4)
+    for (int i = -8; i <= 8; i += 4) {
-        this->energy_diff.insert({i, std::exp(-((double)i) / this->T)});
+        this->energy_diff.insert({i, std::exp(-((double)i / this->T))});
+    }
 }
 
 void IsingModel::initialize_magnetization()
 {
-    this->M = arma::accu(this->lattice);
+    this->M = 0.;
+    for (size_t i = 0; i < this->lattice.n_elem; i++) {
+        this->M += (double)this->lattice(i);
+    }
 }
 
 void IsingModel::initialize_energy()
 {
-    this->E = 0;
+    this->E = 0.;
 
     // Loop through the matrix
     for (size_t j = 0; j < this->L; j++) {
         for (size_t i = 0; i < this->L; i++) {
-            this->E -= this->lattice(i, j)
+            this->E -= (double)this->lattice(i, j)
                        * (this->lattice(i, this->neighbors(j, RIGHT))
                           + this->lattice(this->neighbors(i, DOWN), j));
         }
     }
 }
 
-void IsingModel::Metropolis(std::mt19937 engine)
+data_t IsingModel::Metropolis(std::mt19937 &engine)
 {
     uint ri, rj;
     int dE;
+    data_t res;
 
     // Create random distribution for indeces
-    std::uniform_int_distribution<uint> random_index(0, this->L);
+    std::uniform_int_distribution<uint> random_index(0, this->L - 1);
     // Create random distribution for acceptance
     std::uniform_real_distribution<> random_number(0., 1.);
 
@@ -109,7 +108,7 @@ void IsingModel::Metropolis(std::mt19937 engine)
         rj = random_index(engine);
 
         // Calculate the difference in energy
-        dE = 2. * this->lattice(ri, rj)
+        dE = 2 * this->lattice(ri, rj)
              * (this->lattice(ri, this->neighbors(rj, LEFT))
                 + this->lattice(ri, this->neighbors(rj, RIGHT))
                 + this->lattice(this->neighbors(ri, UP), rj)
@@ -119,18 +118,25 @@ void IsingModel::Metropolis(std::mt19937 engine)
         if (random_number(engine) <= this->energy_diff[dE]) {
             // Update if the configuration is accepted
             this->lattice(ri, rj) *= -1;
-            this->M += 2 * this->lattice(ri, rj);
+            this->M += 2. * (double)this->lattice(ri, rj);
-            this->E += dE;
+            this->E += (double)dE;
         }
     }
+    res.E = this->E;
+    res.E2 = this->E * this->E;
+    res.M = this->M;
+    res.M2 = this->M * this->M;
+    res.M_abs = std::abs(this->M);
+
+    return res;
 }
 
-int IsingModel::get_E()
+double IsingModel::get_E()
 {
     return this->E;
 }
 
-int IsingModel::get_M()
+double IsingModel::get_M()
 {
     return this->M;
 }

src/data_type.cpp

@@ -0,0 +1,72 @@
+/** @file data_type.cpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 1.0
+ *
+ *  @brief Implementation for the data_t type.
+ *
+ *  @bug No known bugs
+ * */
+#include "data_type.hpp"
+
+template <class T>
+data_t operator/(const data_t &data, T num)
+{
+    data_t res = data;
+    res.E /= num;
+    res.E2 /= num;
+    res.M /= num;
+    res.M2 /= num;
+    res.M_abs /= num;
+
+    return res;
+}
+
+// Explicit instantiation
+template data_t operator/(const data_t &, uint);
+template data_t operator/(const data_t &, ulong);
+template data_t operator/(const data_t &,int);
+template data_t operator/(const data_t &,double);
+
+template <class T>
+data_t& operator/=(data_t &data, T num)
+{
+    data.E /= num;
+    data.E2 /= num;
+    data.M /= num;
+    data.M2 /= num;
+    data.M_abs /= num;
+
+    return data;
+}
+
+// Explicit instantiation
+template data_t& operator/=(data_t &, uint);
+template data_t& operator/=(data_t &, ulong);
+template data_t& operator/=(data_t &,int);
+template data_t& operator/=(data_t &,double);
+
+data_t operator+(const data_t &a, const data_t &b)
+{
+    data_t res = a;
+    res.E += b.E;
+    res.E2 += b.E2;
+    res.M += b.M;
+    res.M2 += b.M2;
+    res.M_abs += b.M_abs;
+
+    return res;
+}
+
+data_t& operator+=(data_t &a, const data_t &b)
+{
+    a.E += b.E;
+    a.E2 += b.E2;
+    a.M += b.M;
+    a.M2 += b.M2;
+    a.M_abs += b.M_abs;
+
+    return a;
+}

src/main.cpp

@@ -1,17 +1,43 @@
-#include "IsingModel.hpp"
+/** @file main.cpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 0.1
+ *
+ *  @brief The main program
+ *
+ *  @bug No known bugs
+ * */
+#include "monte_carlo.hpp"
+#include "utils.hpp"
 
 #include <iostream>
 #include <omp.h>
 
-typedef struct data {
+/** @brief The main function.*/
-    double E = 0.;
-    double E_squared = 0.;
-    double M = 0.;
-    double M_squared = 0.;
-    double M_abs = 0.;
-} data_t;
-
 int main()
 {
+    // uint test_cycles = test_2x2_lattice(1e-7, 10000);
+    // monte_carlo(1.0, 2, 10000, "output/2_lattice_test.txt");
+
+    // Test burn-in time
+    monte_carlo_progression(1.0, 20, 10000, "output/burn_in_time_1_0.txt");
+    monte_carlo_progression(2.4, 20, 10000, "output/burn_in_time_2_4.txt");
+
+    // Test the openmp speedup
+    double t0, t1, t2;
+    t0 = omp_get_wtime();
+    phase_transition(20, 2.1, 2.4, 1000, monte_carlo_serial,
+                     "output/phase_transition/size_20.txt");
+    t1 = omp_get_wtime();
+    phase_transition(20, 2.1, 2.4, 1000, monte_carlo_parallel,
+                     "output/phase_transition/size_20.txt");
+    t2 = omp_get_wtime();
+
+    std::cout << "Time serial     : " << t1 - t0 << " seconds" << '\n';
+    std::cout << "Time parallel   : " << t2 - t1 << " seconds" << '\n';
+    std::cout << "Speedup parallel: " << (t1 - t0) / (t2 - t1) << '\n';
+
     return 0;
 }

src/monte_carlo.cpp

@@ -10,3 +10,186 @@
  *  @bug No known bugs
  * */
 #include "monte_carlo.hpp"
+
+uint test_2x2_lattice(double tol, uint max_cycles)
+{
+    data_t data, tmp;
+    size_t L = 2;
+    size_t n_spins = L * L;
+    double T = 1.;
+    size_t cycles = 0;
+
+    // Create random engine using the mersenne twister
+    std::random_device rd;
+    std::mt19937 engine(rd());
+
+    IsingModel test(L, T);
+
+    double E, M;
+
+    std::function<bool(double, double)> is_close = [tol](double a, double b) {
+        return std::abs(a - b) < tol;
+    };
+
+    // Loop through cycles
+    while (cycles++ < max_cycles) {
+        data += test.Metropolis(engine);
+        tmp = data / (cycles * n_spins);
+        // if (close(EPS_2, tmp.E)
+        //&& close(MAG_2, tmp.M)
+        //&& close(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T))
+        //&& close(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T)) {
+        // return cycles;
+        //}
+        if (is_close(EPS_2, tmp.E) && is_close(MAG_2, tmp.M)) {
+            return cycles;
+        }
+    }
+    std::cout << "hello" << std::endl;
+    return 0;
+}
+
+void monte_carlo_progression(double T, uint L, uint cycles,
+                             const std::string filename)
+{
+    // Set some variables
+    data_t data, tmp;
+    uint n_spins = L * L;
+
+    // File stuff
+    std::string directory = utils::dirname(filename);
+    std::ofstream ofile;
+
+    // Create random engine using the mersenne twister
+    std::random_device rd;
+    std::mt19937 engine(rd());
+
+    IsingModel ising(L, T);
+
+    // Create path and open file
+    utils::mkpath(directory);
+    ofile.open(filename);
+
+    // Loop through cycles
+    for (size_t i = 1; i <= cycles; i++) {
+        data += ising.Metropolis(engine);
+        tmp = data / (i * n_spins);
+        ofile << i << ',' << tmp.E << ',' << tmp.E2 << ',' << tmp.M << ','
+              << tmp.M2 << ',' << tmp.M_abs << '\n';
+    }
+    ofile.close();
+}
+
+void pd_estimate(double T, uint L, uint cycles, const std::string filename)
+{
+    // Set some variables
+    data_t data;
+    uint n_spins = L * L;
+
+    // File stuff
+    std::string directory = utils::dirname(filename);
+    std::ofstream ofile;
+
+    // Create random engine using the mersenne twister
+    std::random_device rd;
+    std::mt19937 engine(rd());
+
+    IsingModel ising(L, T);
+
+    // Create path and open file
+    utils::mkpath(directory);
+    ofile.open(filename);
+
+    double E, M;
+
+    // Figure out bin widths and such
+}
+
+// Code for seeing phase transitions.
+void monte_carlo_serial(data_t &data, uint L, double T, uint cycles)
+{
+    IsingModel model(L, T);
+    // Create random engine using the mersenne twister
+    std::random_device rd;
+    std::mt19937 engine(rd());
+
+    for (size_t i = 0; i < BURN_IN_TIME; i++) {
+        model.Metropolis(engine);
+    }
+
+    for (size_t i = 0; i < cycles; i++) {
+        data += model.Metropolis(engine);
+    }
+    data /= cycles;
+}
+
+void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles)
+{
+#pragma omp parallel
+    {
+        // Each thread creates an instance of IsingModel.
+        IsingModel model(L, T);
+        // Each thread creates an instance of the mersenne twister
+        std::random_device rd;
+        std::mt19937 engine(rd());
+
+        data_t tmp;
+
+        // Each thread runs the Metropolis algorithm before starting to collect
+        // samples
+        for (size_t i = 0; i < BURN_IN_TIME; i++) {
+            model.Metropolis(engine);
+        }
+
+        // Now each thread work on one loop together, but using their own
+        // instances of things, but the total of cycles add up.
+        // static ensure that each thread gets the same amount of iterations
+#pragma omp for schedule(static)
+        for (size_t i = 0; i < cycles; i++) {
+            tmp = tmp + model.Metropolis(engine);
+        }
+
+        // Combine all the data.
+#pragma omp critical
+        {
+            data += tmp;
+        }
+    }
+
+    data /= cycles;
+}
+
+void phase_transition(
+    uint L, double start_T, double end_T, uint points_T,
+    std::function<void(data_t &, uint, double, uint)> monte_carlo,
+    std::string outfile)
+{
+    double dt_T = (end_T - start_T) / points_T;
+    uint cycles = 10000;
+    uint N = L * L;
+    std::ofstream ofile;
+
+    data_t data[points_T];
+
+    for (size_t i = 0; i < points_T; i++) {
+        monte_carlo(data[i], L, start_T + dt_T * i, cycles);
+    }
+
+    utils::mkpath(utils::dirname(outfile));
+    ofile.open(outfile);
+
+    double temp, CV, X;
+
+    using utils::scientific_format;
+    for (size_t i = 0; i < points_T; i++) {
+        temp = start_T + dt_T * i;
+        CV = (data[i].E2 - data[i].E * data[i].E) / (N * temp * temp);
+        X = (data[i].M2 - data[i].M_abs * data[i].M_abs) / (N * temp);
+
+        ofile << scientific_format(temp) << ','
+              << scientific_format(data[i].E / N) << ','
+              << scientific_format(data[i].M_abs / N) << ','
+              << scientific_format(CV) << ',' << scientific_format(X) << '\n';
+    }
+    ofile.close();
+}

src/phase_transition_mpi.cpp

@@ -0,0 +1,108 @@
+/** @file phase_transition_mpi.cpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 1.0
+ *
+ *  @brief Sweep over different temperatures and generate data.
+ *
+ *  @bug No known bugs
+ * */
+#include "data_type.hpp"
+#include "monte_carlo.hpp"
+#include "utils.hpp"
+
+#include <algorithm>
+#include <fstream>
+#include <iostream>
+#include <iterator>
+#include <mpi.h>
+
+/** @brief The main function*/
+int main(int argc, char **argv)
+{
+    double start = 1., end = 3.;
+    uint points = 1000, L = 20, N;
+    double dt = (end - start) / points;
+    uint cycles = 10000;
+    N = L * L;
+    std::ofstream ofile;
+
+    data_t data[points];
+
+    // MPI stuff
+    int rank, cluster_size;
+
+    // Initialize MPI
+    MPI_Init(&argc, &argv);
+
+    // Get the cluster size and rank
+    MPI_Comm_size(MPI_COMM_WORLD, &cluster_size);
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+    uint remainder = points % cluster_size;
+    double i_start;
+    uint i_points;
+    // The last
+    if (rank < remainder) {
+        i_points = points / cluster_size + 1;
+        i_start = start + dt * i_points * rank;
+    }
+    else {
+        i_points = points / cluster_size;
+        i_start = start + dt * (i_points * rank + remainder);
+    }
+
+    data_t i_data[i_points];
+    std::cout << "Rank " << rank << ": " << i_points << ',' << i_start << '\n';
+
+    for (size_t i = 0; i < i_points; i++) {
+        monte_carlo_serial(i_data[i], L, i_start + dt * i, cycles);
+    }
+
+    if (rank == 0) {
+        std::copy_n(i_data, i_points, data);
+        for (size_t i = 1; i < cluster_size; i++) {
+            if (rank < remainder) {
+                MPI_Recv((void *)i_data,
+                         sizeof(data_t) * (points / cluster_size + 1), MPI_CHAR,
+                         i, MPI_ANY_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                std::copy_n(i_data, points / cluster_size + 1,
+                            data + (points / cluster_size) * i);
+            }
+            else {
+                MPI_Recv((void *)i_data,
+                         sizeof(data_t) * (points / cluster_size), MPI_CHAR, i,
+                         MPI_ANY_TAG, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+                std::copy_n(i_data, points / cluster_size,
+                            data + (points / cluster_size) * i + remainder);
+            }
+        }
+    }
+    else {
+        MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
+                 MPI_COMM_WORLD);
+    }
+
+    MPI_Finalize();
+
+    std::string outfile = "output/phase_transition/size_20.txt";
+    utils::mkpath(utils::dirname(outfile));
+    ofile.open(outfile);
+
+    double temp, CV, X;
+
+    using utils::scientific_format;
+    for (size_t i = 0; i < points; i++) {
+        temp = start + dt * i;
+        CV = (data[i].E2 - data[i].E * data[i].E) / (N * temp * temp);
+        X = (data[i].M2 - data[i].M_abs * data[i].M_abs) / (N * temp);
+
+        ofile << scientific_format(temp) << ','
+              << scientific_format(data[i].E / N) << ','
+              << scientific_format(data[i].M_abs / N) << ','
+              << scientific_format(CV) << ',' << scientific_format(X) << '\n';
+    }
+    ofile.close();
+}

src/test_suite.cpp

@@ -1,11 +1,23 @@
+/** @file test_suite.cpp
+ *
+ *  @author Cory Alexander Balaton (coryab)
+ *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
+ *
+ *  @version 1.0
+ *
+ *  @brief Sweep over different temperatures and generate data.
+ *
+ *  @bug No known bugs
+ * */
 #include "IsingModel.hpp"
 #include "testlib.hpp"
 
-#include <cxxabi.h>
+/** @brief Test class for the Ising model
-#include <typeinfo>
+ * */
-
 class IsingModelTest {
 public:
+    /** @brief Test That initializing works as intended.
+     * */
     void test_init_functions()
     {
         IsingModel test;
@@ -15,11 +27,12 @@ public:
         // Test that initializing the lattice only yields 1s and -1s.
         test.initialize_lattice();
         std::function<bool(int)> f = [](int x) { return x == 1 || x == -1; };
-        ASSERT(assert_each(f, test.lattice), "Test lattice initialization.");
+        ASSERT(testlib::assert_each(f, test.lattice),
+               "Test lattice initialization.");
 
         test.initialize_neighbors();
         arma::Mat<uint> neighbor_matrix("2, 1 ; 0, 2 ; 1, 0");
-        ASSERT(is_equal(neighbor_matrix, test.neighbors),
+        ASSERT(testlib::is_equal(neighbor_matrix, test.neighbors),
                "Test neighbor matrix.");
 
         // Fill the lattice with 1s to be able to test the next functions.
@@ -27,14 +40,15 @@ public:
 
         // Test the initial magnetization.
         test.initialize_magnetization();
-        ASSERT(test.M == 9, "Test intial magnetization");
+        ASSERT(std::abs(test.M - 9.) < 1e-8, "Test intial magnetization");
 
         // Test that the initial energy is correct
         test.initialize_energy();
-        ASSERT(test.E == -18, "Test initial energy.");
+        ASSERT(std::abs(test.E - (-18)) < 1e-8, "Test initial energy.");
     }
 };
 
+/** @brief The main function.*/
 int main()
 {
     IsingModelTest test;

src/testlib.cpp

@@ -1,4 +1,4 @@
-/** @file utils.cpp
+/** @file testlib.cpp
  *
  *  @author Cory Alexander Balaton (coryab)
  *  @author Janita Ovidie Sandtrøen Willumsen (janitaws)
@@ -9,22 +9,21 @@
  *
  *  @bug No known bugs
  * */
-
 #include "testlib.hpp"
 
-static void print_message(std::string msg)
+namespace details {
+void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
+              int line, std::string msg)
 {
+    std::function<void(const std::string &)> print_message =
+        [](const std::string &msg) {
             if (msg.size() > 0) {
                 std::cout << "message: " << msg << "\n\n";
             }
             else {
                 std::cout << "\n";
             }
-}
+        };
-
-void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
-              int line, std::string msg)
-{
     std::string new_assert(f.size() + (expr ? 4 : 6), '-');
     std::cout << "\x1B[36m" << new_assert << "\033[0m\n";
     std::cout << f << ": ";
@@ -40,3 +39,4 @@ void m_assert(bool expr, std::string expr_str, std::string f, std::string file,
         abort();
     }
 }
+} // namespace details

src/utils.cpp

@@ -9,9 +9,10 @@
  *
  *  @bug No known bugs
  * */
-
 #include "utils.hpp"
 
+namespace utils {
+
 std::string scientific_format(double d, int width, int prec)
 {
     std::stringstream ss;
@@ -46,9 +47,17 @@ bool mkpath(std::string path, int mode)
                 && stat(cur_dir.c_str(), &buf) != 0) {
                 return -1;
             }
-        } else {
+        }
+        else {
             break;
         }
     }
     return 0;
 }
+
+std::string dirname(const std::string &path)
+{
+    return path.substr(0, path.find_last_of("/"));
+}
+
+} // namespace utils