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

Make changes

Browse Source
This commit is contained in:
Cory Balaton
2023-12-03 13:33:48 +01:00
parent 84692040d3
commit 237bd37184
13 changed files with 377 additions and 207 deletions
Download Patch File Download Diff File Expand all files Collapse all files

22
src/IsingModel.cpp
View File

@@ -10,18 +10,18 @@
* @bug No known bugs
* */
#include "IsingModel.hpp"
#include <cmath>
#include <random>
IsingModel::IsingModel()
{
this->initialize_engine();
}
IsingModel::IsingModel(int L, double T)
{
this->L = L;
this->T = T;
this->initialize_engine();
this->initialize_lattice();
this->initialize_neighbors();
this->initialize_energy_diff();
@@ -33,6 +33,7 @@ IsingModel::IsingModel(int L, double T, int val)
{
this->L = L;
this->T = T;
this->initialize_engine();
this->lattice.set_size(this->L, this->L);
this->lattice.fill(val);
this->initialize_neighbors();
@@ -41,16 +42,20 @@ IsingModel::IsingModel(int L, double T, int val)
this->initialize_energy();
}
void IsingModel::initialize_engine()
{
std::random_device rd{};
this->engine = std::mt19937{rd()};
}
void IsingModel::initialize_lattice()
{
this->lattice.set_size(this->L, this->L);
std::random_device rd{};
std::mt19937 engine{rd()};
std::uniform_int_distribution<> coin_flip(0, 1);
for (size_t i = 0; i < this->lattice.n_elem; i++)
this->lattice(i) = 2 * coin_flip(engine) - 1;
this->lattice(i) = 2 * coin_flip(this->engine) - 1;
}
void IsingModel::initialize_neighbors()
@@ -67,7 +72,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[i+8] = std::exp(-(double)i / this->T);
}
}
@@ -95,9 +100,6 @@ void IsingModel::initialize_energy()
data_t IsingModel::Metropolis()
{
std::random_device rd{};
std::mt19937_64 engine{rd()};
int ri, rj;
int dE;
@@ -120,7 +122,7 @@ data_t IsingModel::Metropolis()
+ this->lattice(this->neighbors(ri, DOWN), rj));
// Choose whether or not to accept the new configuration
if (random_number(engine) <= this->energy_diff[dE]) {
if (random_number(engine) <= this->energy_diff[dE+8]) {
// Update if the configuration is accepted
this->lattice(ri, rj) *= -1;
this->M += 2 * this->lattice(ri, rj);

23
src/Makefile
View File

@@ -3,9 +3,11 @@ CC=mpic++
LIBSRCS=utils.cpp testlib.cpp data_type.cpp
LIBOBJS=$(LIBSRCS:.cpp=.o)
LIBPROFOBJS=$(addprefix prof/, $(LIBOBJS))
CLASSSRCS=IsingModel.cpp monte_carlo.cpp
CLASSOBJS=$(CLASSSRCS:.cpp=.o)
CLASSPROFOBJS=$(addprefix prof/, $(CLASSOBJS))
INCLUDE=../include
@@ -29,16 +31,14 @@ else
PROFFLAG=
endif
.PHONY: clean
.PHONY: clean instrument
all: main phase_transition_mpi test_suite
all: main phase_transition_mpi test_suite time
#all: main
# Instrumentation using scorep for parallel analysis
instrument:
scorep $(CC) -c utils.cpp -o utils.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
scorep $(CC) -c main.cpp -o main.o $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
scorep $(CC) $(LIBOBJS) $(CLASSOBJS) main.o -o main $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
instrument: prof/phase_transition_mpi.o $(LIBPROFOBJS) $(CLASSPROFOBJS)
scorep $(CC) $(LIBPROFOBJS) $(CLASSPROFOBJS) $< -o phase_transition_mpi_prof $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
# Rules for executables
main: main.o $(LIBOBJS) $(CLASSOBJS)
@@ -50,10 +50,17 @@ phase_transition_mpi: phase_transition_mpi.o $(LIBOBJS) $(CLASSOBJS)
test_suite: test_suite.o $(LIBOBJS) $(CLASSOBJS)
$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
time: time.o $(LIBOBJS) $(CLASSOBJS)
$(CC) $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
# Rule for object files
%.o: %.cpp
$(CC) -c $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
# Rule for instrumented object files
prof/%.o: %.cpp
scorep $(CC) -c $^ -o $@ $(CFLAGS) $(DBGFLAG) $(PROFFLAG) -I$(INCLUDE) $(OPENMP)
clean:
rm *.o
rm test_suite main
find . -maxdepth 2 -name "*.o" -type f -delete
rm test_suite main phase_transition_mpi phase_transition_mpi_prof time

76
src/main.cpp
View File

@@ -13,31 +13,43 @@
#include "monte_carlo.hpp"
#include "utils.hpp"
#include <csignal>
#include <cstdlib>
#include <iostream>
#include <omp.h>
/** @brief Create the data for the burn-in time for temperatures 1.0 and 2.4
* for both unordered and ordered initial states.
* */
void create_burn_in_time_data()
{
// Test burn-in time
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");
montecarlo::progression(1.0, 20, 20000,
"../output/burn_in_time/unordered_1_0.txt");
montecarlo::progression(1.0, 20, 20000, 1,
"../output/burn_in_time/ordered_1_0.txt");
montecarlo::progression(2.4, 20, 20000,
"../output/burn_in_time/unordered_2_4.txt");
montecarlo::progression(2.4, 20, 20000, 1,
"../output/burn_in_time/ordered_2_4.txt");
}
/** @brief Create the data used to estimate the probability distribution
* for tempratures 1.0 anbd 2.4.
* */
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");
montecarlo::pd_estimate(1.0, 20, 1000000,
"../output/pd_estimate/estimate_1_0.txt");
montecarlo::pd_estimate(2.4, 20, 1000000,
"../output/pd_estimate/estimate_2_4.txt");
}
void test_burn_in_time()
{
montecarlo::phase_transition(100, 2.1, 2.4, 40, 1e5, montecarlo::mcmc_serial,
"../output/test_burn_in_time/no_burn_in.txt", 0);
montecarlo::phase_transition(100, 2.1, 2.4, 40, 1e5, montecarlo::mcmc_serial,
"../output/test_burn_in_time/burn_in.txt", 5000);
}
/** @brief Test how much Openmp speeds up.*/
void test_parallel_speedup()
{
// Test the openmp speedup
@@ -46,10 +58,10 @@ void test_parallel_speedup()
int tries = 5;
t0 = omp_get_wtime();
for (size_t i = 0; i < tries; i++)
monte_carlo_serial(20, 1.0, 10000);
montecarlo::mcmc_serial(20, 1.0, 10000);
t1 = omp_get_wtime();
for (size_t i = 0; i < tries; i++)
monte_carlo_parallel(20, 1.0, 10000);
montecarlo::mcmc_parallel(20, 1.0, 10000);
t2 = omp_get_wtime();
std::cout << "Time serial : " << (t1 - t0) / tries << " seconds"
@@ -59,22 +71,29 @@ void test_parallel_speedup()
std::cout << "Speedup parallel: " << (t1 - t0) / (t2 - t1) << '\n';
}
/** @brief Create data for studying phase transition.
* */
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");
montecarlo::phase_transition(20, 2.1, 2.4, 40, 1e4,
montecarlo::mcmc_parallel,
"../output/phase_transition/size_20.txt");
montecarlo::phase_transition(40, 2.1, 2.4, 40, 1e4,
montecarlo::mcmc_parallel,
"../output/phase_transition/size_40.txt");
montecarlo::phase_transition(60, 2.1, 2.4, 40, 1e4,
montecarlo::mcmc_parallel,
"../output/phase_transition/size_60.txt");
montecarlo::phase_transition(80, 2.1, 2.4, 40, 1e4,
montecarlo::mcmc_parallel,
"../output/phase_transition/size_80.txt");
montecarlo::phase_transition(100, 2.1, 2.4, 40, 1e4,
montecarlo::mcmc_parallel,
"../output/phase_transition/size_100.txt");
t1 = omp_get_wtime();
std::cout << "Time: " << t1 - t0 << std::endl;
@@ -104,6 +123,9 @@ int main(int argc, char **argv)
case 4:
create_phase_transition_data();
break;
case 5:
test_burn_in_time();
break;
default:
std::cout << "Not a valid option!" << std::endl;
abort();

29
src/monte_carlo.cpp
View File

@@ -11,11 +11,8 @@
* */
#include "monte_carlo.hpp"
#include <cmath>
#include <cstdint>
void monte_carlo_progression(double T, int L, int cycles,
const std::string filename)
namespace montecarlo {
void progression(double T, int L, int cycles, const std::string filename)
{
// Set some variables
data_t data, tmp;
@@ -48,8 +45,8 @@ void monte_carlo_progression(double T, int L, int cycles,
ofile.close();
}
void monte_carlo_progression(double T, int L, int cycles, int value,
const std::string filename)
void progression(double T, int L, int cycles, int value,
const std::string filename)
{
// Set some variables
data_t data, tmp;
@@ -108,7 +105,7 @@ void pd_estimate(double T, int L, int cycles, const std::string filename)
}
// Code for seeing phase transitions.
data_t monte_carlo_serial(int L, double T, int cycles)
data_t mcmc_serial(int L, double T, int cycles, int burn_in_time)
{
data_t data;
IsingModel model(L, T);
@@ -122,10 +119,12 @@ data_t monte_carlo_serial(int L, double T, int cycles)
data += model.Metropolis();
}
return data;
double norm = 1. / (double)cycles;
return data * norm;
}
data_t monte_carlo_parallel(int L, double T, int cycles)
data_t mcmc_parallel(int L, double T, int cycles, int burn_in_time)
{
data_t data;
#pragma omp parallel
@@ -153,12 +152,11 @@ data_t monte_carlo_parallel(int L, double T, int cycles)
return data * norm;
}
void phase_transition(int L, double start, double end, int points,
std::function<data_t(int, double, int)> monte_carlo,
std::string outfile)
void phase_transition(int L, double start, double end, int points, int cycles,
std::function<data_t(int, double, int, int)> monte_carlo,
std::string outfile, int burn_in_time)
{
double dt = (end - start) / (double)points;
int cycles = 10000;
int N = L * L;
std::ofstream ofile;
@@ -172,7 +170,7 @@ void phase_transition(int L, double start, double end, int points,
using utils::scientific_format;
for (size_t i = 0; i < points; i++) {
temp = start + dt * i;
data = monte_carlo(L, temp, cycles);
data = monte_carlo(L, temp, cycles, burn_in_time);
E_var = (data.E2 - data.E * data.E) / (double)N;
M_var = (data.M2 - data.M_abs * data.M_abs) / (double)N;
@@ -184,3 +182,4 @@ void phase_transition(int L, double start, double end, int points,
}
ofile.close();
}
} // namespace montecarlo

141
src/phase_transition_mpi.cpp
View File

@@ -7,38 +7,42 @@
*
* @brief Sweep over different temperatures and generate data.
*
* @details This program takes in 4 arguments: the start temperature,
* the end temperature, the amount of temperature points to simulate, and
* the amount of monte carlo samples to collect, in that order.
*
* @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>
#include <sstream>
/** @brief The main function*/
/** @brief The main function
*
* */
int main(int argc, char **argv)
{
if (argc < 5) {
std::cout << "You need at least 4 arguments" << std::endl;
// Check that the number of arguments is at least 4.
if (argc < 7) {
std::cout << "You need at least 6 arguments" << std::endl;
abort();
}
// Timing variables
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;
int cycles = atoi(argv[4]);
std::ofstream ofile;
// Define/initialize variables
double start = atof(argv[1]), end = atof(argv[2]);
int points = atoi(argv[3]), cycles = atoi(argv[5]), L = atoi(argv[4]),
burn_in_time = atoi(argv[6]), N = L * L;
double dt = (end - start) / points;
std::ofstream ofile;
data_t data[points];
// MPI stuff
// MPI specific variables
int rank, cluster_size;
// Initialize MPI
@@ -49,8 +53,9 @@ int main(int argc, char **argv)
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
int remainder = points % cluster_size;
double i_start;
int i_points;
double i_start; // What temperature to start from
int i_points; // How many points to simulate
// Distribute temperature points
if (rank < remainder) {
i_points = points / cluster_size + 1;
@@ -61,61 +66,67 @@ int main(int argc, char **argv)
i_start = start + dt * (i_points * rank + remainder);
}
// Initialize array to contains data for each temperature point
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++) {
i_data[i] = monte_carlo_parallel(L, i_start + dt * i, cycles);
}
// Simulate and save data to array
for (size_t i = 0; i < i_points; i++) {
i_data[i] = montecarlo::mcmc_parallel(L, i_start + dt * i, cycles,
burn_in_time);
}
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);
}
// Rank 0 collects all the data and copies it to the "master"
// data array.
if (rank == 0) {
// Copy its own i_data to the data array
std::copy_n(i_data, i_points, data);
// Collect i_data from other ranks in order and copy to 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);
}
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)
/ ((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';
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);
}
ofile.close();
}
else {
MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
MPI_COMM_WORLD);
// Write everything from data to file
std::stringstream outfile;
outfile << "../output/phase_transition/mpi/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)
/ ((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';
}
ofile.close();
}
// For all other ranks, send the data to rank 0
else {
MPI_Send(i_data, i_points * sizeof(data_t), MPI_CHAR, 0, rank,
MPI_COMM_WORLD);
}
t1 = MPI_Wtime();

61
src/test_suite.cpp
View File

@@ -12,8 +12,6 @@
#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))
@@ -29,7 +27,7 @@
* */
class IsingModelTest {
public:
/** @brief Test That initializing works as intended.
/** @brief Test that initializing works as intended.
* */
void test_init_functions()
{
@@ -84,15 +82,9 @@ public:
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)
@@ -102,44 +94,6 @@ public:
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;
}
};
@@ -150,18 +104,23 @@ int main()
IsingModelTest test;
test.test_init_functions();
int res = 0;
int tmp;
for (size_t i=0; i < 1000; i++) {
int iterations = 10000;
int accepted_values = 0;
// Run through the test multiple times to get a better estimate.
for (size_t i=0; i < iterations; i++) {
tmp = test.test_2x2_lattice(1e-2, 1e5);
if (tmp == 0) {
std::cout << "not enough cycles\n";
break;
continue;
}
accepted_values++;
res += tmp;
}
std::cout << "Res: " << res / 1000 << std::endl;
std::cout << "Res: " << res / accepted_values << std::endl;
return 0;
}

68
src/time.cpp Normal file
View File

@@ -0,0 +1,68 @@
/** @file time.cpp
*
* @author Cory Alexander Balaton (coryab)
* @author Janita Ovidie Sandtrøen Willumsen (janitaws)
*
* @version 0.1
*
* @brief Timing various things
*
* @bug No known bugs
* */
#include "data_type.hpp"
#include "monte_carlo.hpp"
#include "utils.hpp"
#include <omp.h>
#include <ostream>
void time_lattice_sizes()
{
std::string outfile = "output/timing/lattice_sizes.txt";
std::ofstream ofile;
int lattice_sizes[] = {20, 40, 60, 80, 100};
utils::mkpath(utils::dirname(outfile));
ofile.open(outfile);
double t0, t1;
for (int L : lattice_sizes) {
t0 = omp_get_wtime();
montecarlo::phase_transition(L, 2.1, 2.4, 40, 100000,
montecarlo::mcmc_parallel,
"output/garbage/null.txt");
t1 = omp_get_wtime();
ofile << utils::scientific_format(L) << ','
<< utils::scientific_format(t1 - t0) << '\n';
}
ofile.close();
}
void time_sample_sizes()
{
std::string outfile = "output/timing/sample_sizes.txt";
std::ofstream ofile;
int sample_sizes[] = {1000, 10000, 100000};
utils::mkpath(utils::dirname(outfile));
ofile.open(outfile);
double t0, t1;
for (int samples : sample_sizes) {
t0 = omp_get_wtime();
montecarlo::phase_transition(20, 2.1, 2.4, 40, samples,
montecarlo::mcmc_parallel,
"output/garbage/null.txt");
t1 = omp_get_wtime();
ofile << utils::scientific_format(samples) << ','
<< utils::scientific_format(t1 - t0) << '\n';
}
ofile.close();
}
int main()
{
time_lattice_sizes();
time_sample_sizes();
return 0;
}