From 090510175b45c68fc7a96060fcdccd3fce8bb78c Mon Sep 17 00:00:00 2001 From: Cory Date: Wed, 8 Nov 2023 17:31:14 +0100 Subject: [PATCH] Implementations of different problems --- include/IsingModel.hpp | 34 +++---- include/data_type.hpp | 78 +++++++++++++++ include/monte_carlo.hpp | 84 +++++++++++++++- include/testlib.hpp | 34 +++++-- include/utils.hpp | 56 +++++++---- src/IsingModel.cpp | 44 +++++---- src/data_type.cpp | 72 ++++++++++++++ src/main.cpp | 44 +++++++-- src/monte_carlo.cpp | 183 +++++++++++++++++++++++++++++++++++ src/phase_transition_mpi.cpp | 108 +++++++++++++++++++++ src/test_suite.cpp | 28 ++++-- src/testlib.cpp | 24 ++--- src/utils.cpp | 13 ++- 13 files changed, 704 insertions(+), 98 deletions(-) create mode 100644 include/data_type.hpp create mode 100644 src/data_type.cpp create mode 100644 src/phase_transition_mpi.cpp diff --git a/include/IsingModel.hpp b/include/IsingModel.hpp index d4adc9b..9f71bb7 100644 --- a/include/IsingModel.hpp +++ b/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 #include @@ -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 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 - * \f$ L \cross 2 \f$, is that we can see that we can use the same column + * @details The reason why it's \f$ L \times 2 \f$ instead of + * \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; - - /** @brief Energy per spin. unit: \f$ J \f$. - * */ - double eps; - - /** @brief Magnetization per spin. unit: Unitless. - * */ - double m; + 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 diff --git a/include/data_type.hpp b/include/data_type.hpp new file mode 100644 index 0000000..7382db2 --- /dev/null +++ b/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 +#include + +/** @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 +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 +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 diff --git a/include/monte_carlo.hpp b/include/monte_carlo.hpp index 2ad1738..7be8f4f 100644 --- a/include/monte_carlo.hpp +++ b/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 +#include + +#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 monte_carlo, + std::string outfile); #endif diff --git a/include/testlib.hpp b/include/testlib.hpp index b44ec3f..452088c 100644 --- a/include/testlib.hpp +++ b/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 &a, arma::Mat &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 ::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,10 +114,11 @@ static bool is_equal(arma::Mat &a, arma::Mat &b) } return true; } + /** @brief Test that all elements fulfill the condition. * * @param expr The boolean expression to apply to each element - * @param M The matrix/vector to iterate over + * @param M The matrix/vector to iterate over * * @return bool * */ @@ -111,5 +133,5 @@ static bool assert_each(std::function expr, arma::Mat &M) } return true; } - +} // namespace testlib #endif diff --git a/include/utils.hpp b/include/utils.hpp index 9283735..8cd0554 100644 --- a/include/utils.hpp +++ b/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 &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 diff --git a/src/IsingModel.cpp b/src/IsingModel.cpp index dda98b9..f61aa91 100644 --- a/src/IsingModel.cpp +++ b/src/IsingModel.cpp @@ -9,14 +9,8 @@ * * @bug No known bugs * */ - #include "IsingModel.hpp" -#include "constants.hpp" - -#include -#include - 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 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) - this->energy_diff.insert({i, std::exp(-((double)i) / this->T)}); + for (int i = -8; i <= 8; i += 4) { + 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 random_index(0, this->L); + std::uniform_int_distribution 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->E += dE; + this->M += 2. * (double)this->lattice(ri, rj); + 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; } diff --git a/src/data_type.cpp b/src/data_type.cpp new file mode 100644 index 0000000..e602fe8 --- /dev/null +++ b/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 +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 +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; +} diff --git a/src/main.cpp b/src/main.cpp index 221aaee..dc4a72d 100644 --- a/src/main.cpp +++ b/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 #include -typedef struct data { - double E = 0.; - double E_squared = 0.; - double M = 0.; - double M_squared = 0.; - double M_abs = 0.; -} data_t; - +/** @brief The main function.*/ 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; } diff --git a/src/monte_carlo.cpp b/src/monte_carlo.cpp index 7bd596d..b3bb867 100644 --- a/src/monte_carlo.cpp +++ b/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 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 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(); +} diff --git a/src/phase_transition_mpi.cpp b/src/phase_transition_mpi.cpp new file mode 100644 index 0000000..ac925a0 --- /dev/null +++ b/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 +#include +#include +#include +#include + +/** @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(); +} diff --git a/src/test_suite.cpp b/src/test_suite.cpp index 4d3437d..59ec82c 100644 --- a/src/test_suite.cpp +++ b/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 -#include - +/** @brief Test class for the Ising model + * */ 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 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 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; diff --git a/src/testlib.cpp b/src/testlib.cpp index 7cce040..1027b24 100644 --- a/src/testlib.cpp +++ b/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) -{ - if (msg.size() > 0) { - std::cout << "message: " << msg << "\n\n"; - } - else { - std::cout << "\n"; - } -} - +namespace details { void m_assert(bool expr, std::string expr_str, std::string f, std::string file, int line, std::string msg) { + std::function print_message = + [](const std::string &msg) { + if (msg.size() > 0) { + std::cout << "message: " << msg << "\n\n"; + } + else { + std::cout << "\n"; + } + }; 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 diff --git a/src/utils.cpp b/src/utils.cpp index bf602fc..b51fa72 100644 --- a/src/utils.cpp +++ b/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