FYS3150/Project-4
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Update code

Browse Source
This commit is contained in:
Cory Balaton 2023-11-21 11:24:34 +01:00
parent 816e38e9e4
commit f07fb8829b
No known key found for this signature in database
GPG Key ID: 3E5FCEBFD80F432B
12 changed files with 522 additions and 343 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
include/IsingModel.hpp
View File

@ -21,7 +21,6 @@
#include <random>
#include <unordered_map>
// Faster modulo
#define INDEX(I, N) (I + N) % N
// Indeces for the neighbor matrix.
@ -50,7 +49,7 @@ private:
* for the left and upper neighbor, and we can use the same column for the
* right and lower neighbor.
* */
arma::Mat<uint> neighbors;
arma::Mat<int> neighbors;
/** @brief A hash map containing all possible energy changes.
* */
@ -62,15 +61,15 @@ private:
/** @brief Size of the lattice.
* */
uint L;
int L;
/** @brief The current energy state. unit: \f$ J \f$.
* */
double E;
int E;
/** @brief The current magnetic strength. unit: Unitless.
* */
double M;
int M;
/** @brief Initialize the lattice with a random distribution of 1s and
* -1s.
@ -103,7 +102,7 @@ public:
* @param L The size of the lattice.
* @param T The temperature for the system.
* */
IsingModel(uint L, double T);
IsingModel(int L, double T);
/** @brief Constructor for the Ising model.
*
@ -111,23 +110,23 @@ public:
* @param T The temperature for the system.
* @param val The value to set for all spins.
* */
IsingModel(uint L, double T, int val);
IsingModel(int L, double T, int val);
/** @brief The Metropolis algorithm.
* */
data_t Metropolis(std::mt19937 &engine);
data_t Metropolis();
/** @brief Get the current energy.
*
* @return double
* */
double get_E();
int get_E();
/** @brief Get the current magnetization.
*
* @return double
* */
double get_M();
int get_M();
};
#endif

157
include/data_type.hpp
View File

@ -15,64 +15,105 @@
#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
class data_t {
public:
double E, M, E2, M2, M_abs;
data_t()
{
this->E = 0.;
this->E2 = 0.;
this->M = 0.;
this->M2 = 0.;
this->M_abs = 0.;
}
data_t(double E, double E2, double M, double M2, double M_abs)
{
this->E = E;
this->E2 = E2;
this->M = M;
this->M2 = M2;
this->M_abs = M_abs;
}
template <class T> data_t operator/(T num)
{
data_t res;
res.E = this->E / (double)num;
res.E2 = this->E2 / (double)num;
res.M = this->M / (double)num;
res.M2 = this->M2 / (double)num;
res.M_abs = this->M_abs / (double)num;
return res;
}
template <class T> data_t& operator/=(T num)
{
this->E /= (double)num;
this->E2 /= (double)num;
this->M /= (double)num;
this->M2 /= (double)num;
this->M_abs /= (double)num;
return *this;
}
template <class T> data_t operator*(T num)
{
data_t res;
res.E = this->E * (double)num;
res.E2 = this->E2 * (double)num;
res.M = this->M * (double)num;
res.M2 = this->M2 * (double)num;
res.M_abs = this->M_abs * (double)num;
return res;
}
template <class T> data_t& operator*=(T num)
{
this->E *= (double)num;
this->E2 *= (double)num;
this->M *= (double)num;
this->M2 *= (double)num;
this->M_abs *= (double)num;
return *this;
}
data_t operator+(const data_t &b)
{
data_t res;
res.E = this->E + b.E;
res.E2 = this->E2 + b.E2;
res.M = this->M + b.M;
res.M2 = this->M2 + b.M2;
res.M_abs = this->M_abs + b.M_abs;
return res;
}
data_t& operator+=(const data_t &b)
{
this->E += b.E;
this->E2 += b.E2;
this->M += b.M;
this->M2 += b.M2;
this->M_abs += b.M_abs;
return *this;
}
template <class T> void operator=(T num)
{
this->E = (double)num;
this->E2 = (double)num;
this->M = (double)num;
this->M2 = (double)num;
this->M_abs = (double)num;
}
};
/** @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

54
include/monte_carlo.hpp
View File

@ -18,31 +18,21 @@
#include <functional>
#include <string>
#include <omp.h>
#define BURN_IN_TIME 1000
//#define BURN_IN_TIME 12500
#define BURN_IN_TIME 5000
#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))
#pragma omp declare reduction(+: data_t: omp_out += omp_in)
/** @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
* return int
* */
uint test_2x2_lattice(double tol, uint max_cycles);
int test_2x2_lattice(double tol, int max_cycles);
/** @brief Write the expected values for each Monte Carlo cycles to file.
*
@ -51,7 +41,18 @@ uint test_2x2_lattice(double tol, uint max_cycles);
* @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,
void monte_carlo_progression(double T, int L, int cycles,
const std::string filename);
/** @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 value The value to set the elements in the lattice
* @param filename The file to write to
* */
void monte_carlo_progression(double T, int L, int cycles, int value,
const std::string filename);
/** @brief Estimate the probability distribution for the energy.
@ -61,26 +62,29 @@ void monte_carlo_progression(double T, uint L, uint cycles,
* @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);
void pd_estimate(double T, int L, int 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);
* @param cycles The amount of Monte Carlo cycles to do
*
* @return data_t
* */
data_t monte_carlo_serial(int L, double T, int 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
*
* @return data_t
* */
void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles);
data_t monte_carlo_parallel(int L, double T, int cycles);
/** @brief Perform the MCMC algorithm using a range of temperatures.
*
@ -92,8 +96,8 @@ void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles);
* @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,
int L, double start_T, double end_T, int points_T,
std::function<data_t(int, double, int)> monte_carlo,
std::string outfile);
#endif

5
include/testlib.hpp
View File

@ -8,7 +8,8 @@
* @brief A small test library.
*
* @details This a small testing library that is tailored for the needs of the
* project.
* project. Anything that is in the details namespace should not be used
* directly, or else it might cause undefined behavior if not used correctly.
*
* @bug No known bugs
* */
@ -90,7 +91,7 @@ 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;
return std::fabs(a - b) < tol;
}
/** @brief Test if two armadillo matrices/vectors are equal.

4
include/utils.hpp
View File

@ -8,7 +8,9 @@
* @brief Function prototypes and macros that are useful.
*
* These utility function are mainly for convenience and aren't directly
* related to the project.
* related to the project. Anything that is in the details namespace should
* not be used directly, or else it might cause undefined behavior if not used
* correctly.
*
* @bug No known bugs
* */

47
src/IsingModel.cpp
View File

@ -11,11 +11,14 @@
* */
#include "IsingModel.hpp"
#include <cmath>
#include <random>
IsingModel::IsingModel()
{
}
IsingModel::IsingModel(uint L, double T)
IsingModel::IsingModel(int L, double T)
{
this->L = L;
this->T = T;
@ -26,7 +29,7 @@ IsingModel::IsingModel(uint L, double T)
this->initialize_energy();
}
IsingModel::IsingModel(uint L, double T, int val)
IsingModel::IsingModel(int L, double T, int val)
{
this->L = L;
this->T = T;
@ -41,13 +44,13 @@ 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::mt19937 engine(rd());
std::random_device rd{};
std::mt19937 engine{rd()};
std::uniform_int_distribution<int> coin_flip(0, 1);
std::uniform_int_distribution<> coin_flip(0, 1);
for (size_t i = 0; i < this->lattice.n_elem; i++)
this->lattice(i) = coin_flip(engine) == 1 ? 1 : -1;
this->lattice(i) = 2 * coin_flip(engine) - 1;
}
void IsingModel::initialize_neighbors()
@ -64,7 +67,7 @@ void IsingModel::initialize_neighbors()
void IsingModel::initialize_energy_diff()
{
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)});
}
}
@ -72,7 +75,7 @@ void IsingModel::initialize_magnetization()
{
this->M = 0.;
for (size_t i = 0; i < this->lattice.n_elem; i++) {
this->M += (double)this->lattice(i);
this->M += this->lattice(i);
}
}
@ -83,21 +86,23 @@ void IsingModel::initialize_energy()
// Loop through the matrix
for (size_t j = 0; j < this->L; j++) {
for (size_t i = 0; i < this->L; i++) {
this->E -= (double)this->lattice(i, j)
this->E -= this->lattice(i, j)
* (this->lattice(i, this->neighbors(j, RIGHT))
+ this->lattice(this->neighbors(i, DOWN), j));
}
}
}
data_t IsingModel::Metropolis(std::mt19937 &engine)
data_t IsingModel::Metropolis()
{
uint ri, rj;
std::random_device rd{};
std::mt19937_64 engine{rd()};
int ri, rj;
int dE;
data_t res;
// Create random distribution for indeces
std::uniform_int_distribution<uint> random_index(0, this->L - 1);
std::uniform_int_distribution<> random_index(0, this->L - 1);
// Create random distribution for acceptance
std::uniform_real_distribution<> random_number(0., 1.);
@ -118,25 +123,21 @@ data_t 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. * (double)this->lattice(ri, rj);
this->E += (double)dE;
this->M += 2 * this->lattice(ri, rj);
this->E += 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;
return data_t((double)this->E, (double)(this->E * this->E), (double)this->M,
(double)(this->M * this->M), std::fabs((double)this->M));
}
double IsingModel::get_E()
int IsingModel::get_E()
{
return this->E;
}
double IsingModel::get_M()
int IsingModel::get_M()
{
return this->M;
}

2
src/Makefile
View File

@ -16,7 +16,7 @@ OPENMP=-fopenmp
# Add a debug flag when compiling (For the DEBUG macro in utils.hpp)
DEBUG ?= 0
ifeq ($(DEBUG), 1)
DBGFLAG=-DDBG
DBGFLAG=-DDBG -g
else
DBGFLAG=
endif

60
src/data_type.cpp
View File

@ -10,63 +10,3 @@
* @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;
}

96
src/main.cpp
View File

@ -9,35 +9,105 @@
*
* @bug No known bugs
* */
#include "data_type.hpp"
#include "monte_carlo.hpp"
#include "utils.hpp"
#include <csignal>
#include <cstdlib>
#include <iostream>
#include <omp.h>
/** @brief The main function.*/
int main()
void create_burn_in_time_data()
{
// 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");
monte_carlo_progression(1.0, 20, 20000,
"output/burn_in_time/unordered_1_0.txt");
monte_carlo_progression(1.0, 20, 20000, 1,
"output/burn_in_time/ordered_1_0.txt");
monte_carlo_progression(2.4, 20, 20000,
"output/burn_in_time/unordered_2_4.txt");
monte_carlo_progression(2.4, 20, 20000, 1,
"output/burn_in_time/ordered_2_4.txt");
}
void create_pd_estimate_data()
{
// Estimate pd
pd_estimate(1.0, 20, 1000000, "output/pd_estimate/estimate_1_0.txt");
pd_estimate(2.4, 20, 1000000, "output/pd_estimate/estimate_2_4.txt");
}
void test_parallel_speedup()
{
// Test the openmp speedup
data_t data;
double t0, t1, t2;
int tries = 5;
t0 = omp_get_wtime();
phase_transition(20, 2.1, 2.4, 1000, monte_carlo_serial,
"output/phase_transition/size_20.txt");
for (size_t i = 0; i < tries; i++)
monte_carlo_serial(20, 1.0, 10000);
t1 = omp_get_wtime();
phase_transition(20, 2.1, 2.4, 1000, monte_carlo_parallel,
"output/phase_transition/size_20.txt");
for (size_t i = 0; i < tries; i++)
monte_carlo_parallel(20, 1.0, 10000);
t2 = omp_get_wtime();
std::cout << "Time serial : " << t1 - t0 << " seconds" << '\n';
std::cout << "Time parallel : " << t2 - t1 << " seconds" << '\n';
std::cout << "Time serial : " << (t1 - t0) / tries << " seconds"
<< '\n';
std::cout << "Time parallel : " << (t2 - t1) / tries << " seconds"
<< '\n';
std::cout << "Speedup parallel: " << (t1 - t0) / (t2 - t1) << '\n';
}
void create_phase_transition_data()
{
double t0, t1;
t0 = omp_get_wtime();
// Phase transition
phase_transition(20, 2.1, 2.4, 40, monte_carlo_parallel,
"output/phase_transition/size_20.txt");
phase_transition(40, 2.1, 2.4, 40, monte_carlo_parallel,
"output/phase_transition/size_40.txt");
phase_transition(60, 2.1, 2.4, 40, monte_carlo_parallel,
"output/phase_transition/size_60.txt");
phase_transition(80, 2.1, 2.4, 40, monte_carlo_parallel,
"output/phase_transition/size_80.txt");
phase_transition(100, 2.1, 2.4, 40, monte_carlo_parallel,
"output/phase_transition/size_100.txt");
t1 = omp_get_wtime();
std::cout << "Time: " << t1 - t0 << std::endl;
}
/** @brief The main function.*/
int main(int argc, char **argv)
{
if (argc < 2) {
std::cout << "Need at least 1 argument, got " << argc - 1
<< " arguments." << std::endl;
abort();
}
int arg = atoi(argv[1]);
switch (arg) {
case 1:
create_burn_in_time_data();
break;
case 2:
create_pd_estimate_data();
break;
case 3:
test_parallel_speedup();
break;
case 4:
create_phase_transition_data();
break;
default:
std::cout << "Not a valid option!" << std::endl;
abort();
}
return 0;
}

173
src/monte_carlo.cpp
View File

@ -11,50 +11,15 @@
* */
#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;
#include <cmath>
#include <cstdint>
// 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,
void monte_carlo_progression(double T, int L, int cycles,
const std::string filename)
{
// Set some variables
data_t data, tmp;
uint n_spins = L * L;
int n_spins = L * L;
// File stuff
std::string directory = utils::dirname(filename);
@ -62,7 +27,10 @@ void monte_carlo_progression(double T, uint L, uint cycles,
// Create random engine using the mersenne twister
std::random_device rd;
std::mt19937 engine(rd());
uint32_t rd_sample = rd();
std::mt19937 engine(rd_sample);
std::cout << "Seed: " << rd_sample << std::endl;
IsingModel ising(L, T);
@ -72,7 +40,7 @@ void monte_carlo_progression(double T, uint L, uint cycles,
// Loop through cycles
for (size_t i = 1; i <= cycles; i++) {
data += ising.Metropolis(engine);
data += ising.Metropolis();
tmp = data / (i * n_spins);
ofile << i << ',' << tmp.E << ',' << tmp.E2 << ',' << tmp.M << ','
<< tmp.M2 << ',' << tmp.M_abs << '\n';
@ -80,11 +48,12 @@ void monte_carlo_progression(double T, uint L, uint cycles,
ofile.close();
}
void pd_estimate(double T, uint L, uint cycles, const std::string filename)
void monte_carlo_progression(double T, int L, int cycles, int value,
const std::string filename)
{
// Set some variables
data_t data;
uint n_spins = L * L;
data_t data, tmp;
int n_spins = L * L;
// File stuff
std::string directory = utils::dirname(filename);
@ -92,104 +61,126 @@ void pd_estimate(double T, uint L, uint cycles, const std::string filename)
// Create random engine using the mersenne twister
std::random_device rd;
uint32_t rd_sample = rd();
std::mt19937 engine(rd());
std::cout << "Seed: " << rd_sample << std::endl;
IsingModel ising(L, T, value);
// 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();
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, int L, int cycles, const std::string filename)
{
// Set some variables
data_t data, tmp;
int n_spins = L * L;
// File stuff
std::string directory = utils::dirname(filename);
std::ofstream ofile;
IsingModel ising(L, T);
// Create path and open file
utils::mkpath(directory);
ofile.open(filename);
double E, M;
// Figure out bin widths and such
// Loop through cycles
for (size_t i = 1; i <= cycles; i++) {
data = ising.Metropolis() / n_spins;
ofile << data.E << ',' << data.E2 << ',' << data.M << ',' << data.M2
<< ',' << data.M_abs << '\n';
}
ofile.close();
}
// Code for seeing phase transitions.
void monte_carlo_serial(data_t &data, uint L, double T, uint cycles)
data_t monte_carlo_serial(int L, double T, int cycles)
{
data_t data;
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);
model.Metropolis();
}
for (size_t i = 0; i < cycles; i++) {
data += model.Metropolis(engine);
double E, M;
for (size_t i = 0; i < (uint)cycles; i++) {
data += model.Metropolis();
}
data /= cycles;
return data;
}
void monte_carlo_parallel(data_t &data, uint L, double T, uint cycles)
data_t monte_carlo_parallel(int L, double T, int cycles)
{
data_t data;
#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);
model.Metropolis();
}
// 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;
#pragma omp for schedule(static) reduction(+ : data)
for (size_t i = 0; i < (uint)cycles; i++) {
data += model.Metropolis();
}
}
data /= cycles;
double norm = 1. / (double)cycles;
return data * norm;
}
void phase_transition(
uint L, double start_T, double end_T, uint points_T,
std::function<void(data_t &, uint, double, uint)> monte_carlo,
void phase_transition(int L, double start, double end, int points,
std::function<data_t(int, double, int)> monte_carlo,
std::string outfile)
{
double dt_T = (end_T - start_T) / points_T;
uint cycles = 10000;
uint N = L * L;
double dt = (end - start) / (double)points;
int cycles = 10000;
int 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);
}
data_t data;
utils::mkpath(utils::dirname(outfile));
ofile.open(outfile);
double temp, CV, X;
double temp, CV, X, E_var, M_var;
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);
for (size_t i = 0; i < points; i++) {
temp = start + dt * i;
data = monte_carlo(L, temp, cycles);
E_var = (data.E2 - data.E * data.E) / (double)N;
M_var = (data.M2 - data.M_abs * data.M_abs) / (double)N;
ofile << scientific_format(temp) << ','
<< scientific_format(data[i].E / N) << ','
<< scientific_format(data[i].M_abs / N) << ','
<< scientific_format(CV) << ',' << scientific_format(X) << '\n';
<< scientific_format(data.E / (double)N) << ','
<< scientific_format(data.M_abs / N) << ','
<< scientific_format(E_var / (temp * temp)) << ','
<< scientific_format(M_var / temp) << '\n';
}
ofile.close();
}

72
src/phase_transition_mpi.cpp
View File

@ -18,15 +18,22 @@
#include <iostream>
#include <iterator>
#include <mpi.h>
#include <sstream>
/** @brief The main function*/
int main(int argc, char **argv)
{
double start = 1., end = 3.;
uint points = 1000, L = 20, N;
if (argc < 5) {
std::cout << "You need at least 4 arguments" << std::endl;
abort();
}
double t0, t1;
t0 = MPI_Wtime();
double start = atof(argv[1]), end = atof(argv[2]);
int points = atoi(argv[3]), N;
int lattice_sizes[] = {20, 40, 60, 80, 100};
double dt = (end - start) / points;
uint cycles = 10000;
N = L * L;
int cycles = atoi(argv[4]);
std::ofstream ofile;
data_t data[points];
@ -41,10 +48,10 @@ int main(int argc, char **argv)
MPI_Comm_size(MPI_COMM_WORLD, &cluster_size);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
uint remainder = points % cluster_size;
int remainder = points % cluster_size;
double i_start;
uint i_points;
// The last
int i_points;
// Distribute temperature points
if (rank < remainder) {
i_points = points / cluster_size + 1;
i_start = start + dt * i_points * rank;
@ -57,8 +64,10 @@ int main(int argc, char **argv)
data_t i_data[i_points];
std::cout << "Rank " << rank << ": " << i_points << ',' << i_start << '\n';
for (int L : lattice_sizes) {
N = L * L;
for (size_t i = 0; i < i_points; i++) {
monte_carlo_serial(i_data[i], L, i_start + dt * i, cycles);
i_data[i] = monte_carlo_parallel(L, i_start + dt * i, cycles);
}
if (rank == 0) {
@ -66,43 +75,54 @@ int main(int argc, char **argv)
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);
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);
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);
std::stringstream outfile;
outfile << "output/phase_transition/size_" << L << ".txt";
utils::mkpath(utils::dirname(outfile.str()));
ofile.open(outfile.str());
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);
CV = (data[i].E2 - data[i].E * data[i].E)
/ ((double)N * temp * temp);
X = (data[i].M2 - data[i].M_abs * data[i].M_abs)
/ ((double)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';
<< scientific_format(CV) << ',' << scientific_format(X)
<< '\n';
}
ofile.close();
}
else {
MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
MPI_COMM_WORLD);
}
}
t1 = MPI_Wtime();
if (rank == 0) {
std::cout << "Time: " << t1 - t0 << " seconds\n";
}
MPI_Finalize();
}

116
src/test_suite.cpp
View File

@ -12,6 +12,19 @@
#include "IsingModel.hpp"
#include "testlib.hpp"
#include <fstream>
#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.) + 1) / ((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 class for the Ising model
* */
class IsingModelTest {
@ -31,7 +44,7 @@ public:
"Test lattice initialization.");
test.initialize_neighbors();
arma::Mat<uint> neighbor_matrix("2, 1 ; 0, 2 ; 1, 0");
arma::Mat<int> neighbor_matrix("2, 1 ; 0, 2 ; 1, 0");
ASSERT(testlib::is_equal(neighbor_matrix, test.neighbors),
"Test neighbor matrix.");
@ -40,11 +53,94 @@ public:
// Test the initial magnetization.
test.initialize_magnetization();
ASSERT(std::abs(test.M - 9.) < 1e-8, "Test intial magnetization");
ASSERT(std::fabs(test.M - 9.) < 1e-8, "Test intial magnetization");
// Test that the initial energy is correct
test.initialize_energy();
ASSERT(std::abs(test.E - (-18)) < 1e-8, "Test initial energy.");
ASSERT(std::fabs(test.E - (-18)) < 1e-8, "Test initial energy.");
}
/** @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 int
* */
int test_2x2_lattice(double tol, int 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);
int arr[]{0, 0, 0, 0};
// Loop through cycles
//std::ofstream ofile;
//ofile.open("output/test_2x2.txt");
while (cycles++ < max_cycles) {
data += test.Metropolis();
tmp = data / cycles;
//ofile << cycles << ',' << tmp.E / n_spins << ','
//<< tmp.M_abs / n_spins << ','
//<< (tmp.E2 - tmp.E * tmp.E) / (T * T) / n_spins << ','
//<< (tmp.M2 - tmp.M_abs * tmp.M_abs) / T / n_spins << '\n';
if (testlib::close_to(EPS_2, tmp.E / n_spins, tol)
&& testlib::close_to(MAG_2, tmp.M_abs / n_spins, tol)
&& testlib::close_to(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T)
/ n_spins, tol)
&& testlib::close_to(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T
/ n_spins, tol)) {
return cycles;
}
}
//std::cout << EPS_2 << ',' << MAG_2 << ',' << CV_2 << ',' << X_2
//<< std::endl;
//ofile.close();
// cycles = 0;
// data = 0;
// IsingModel test_mag(L, T);
// while (cycles++ < max_cycles) {
// data += test.Metropolis();
// tmp = data / (cycles * n_spins);
// if (testlib::close_to(MAG_2, tmp.M, tol)) {
// arr[1] = cycles;
// break;
//}
//}
// cycles = 0;
// data = 0;
// IsingModel test_CV(L, T);
// while (cycles++ < max_cycles) {
// data += test.Metropolis();
// tmp = data / (cycles * n_spins);
// if (testlib::close_to(CV_2, (tmp.E2 - tmp.E * tmp.E) / (T * T),
// tol)) {
// arr[2] = cycles;
// break;
//}
//}
// cycles = 0;
// data = 0;
// IsingModel test_X(L, T);
// while (cycles++ < max_cycles) {
// data += test.Metropolis();
// tmp = data / (cycles * n_spins);
// if (testlib::close_to(X_2, (tmp.M2 - tmp.M_abs * tmp.M_abs) / T,
// tol)) {
// arr[3] = cycles;
// break;
//}
//}
return 0;
}
};
@ -54,4 +150,18 @@ int main()
IsingModelTest test;
test.test_init_functions();
int res = 0;
int tmp;
for (size_t i=0; i < 1000; i++) {
tmp = test.test_2x2_lattice(1e-2, 1e5);
if (tmp == 0) {
std::cout << "not enough cycles\n";
break;
}
res += tmp;
}
std::cout << "Res: " << res / 1000 << std::endl;
return 0;
}