Main part of code for project 4

2025-11-12 20:45:39 +01:00
parent 1dde2c2594
commit 0d5b6e2267
9 changed files with 2760 additions and 1 deletions
@@ -291,4 +291,10 @@ src/project3/results/*/*/*.red
 *.tar.gz
 *_timing.txt
 run_stats.json
-src/project3/python/plots/
+src/project3/python/plots/
+
+
+# Project 4 executables
+src/project4/phase_transition/
+src/project4/system_evolutions/
+src/project4/main
@@ -0,0 +1,62 @@
+#include <vector>
+#include <random>
+
+#ifndef ISING_HPP
+#define ISING_HPP
+
+
+const double kB = 1.0; // Boltzmann constant
+
+class IsingModel {
+public:
+    // Construction & Initialization
+    IsingModel(int size, double temperature, bool randomInit = true, int seed = 42);
+    void initializeLattice(bool randomInit = true);
+
+    // Simulation Control
+    void monteCarloCycle();
+    void run(int cycles);
+
+    // Measurements & Observables
+    double currentEnergy(bool normalized = true) const;
+    double currentMagnetization(bool normalized = true) const;
+
+    double calculateEnergy(bool normalized = true) const;
+    double calculateMagnetization(bool normalized = true) const;
+
+    double calculateEnergyVariance(bool normalized = true) const;
+    double calculateMagnetizationVariance(bool normalized = true) const;
+
+    double calculateHeatCapacity(bool normalized = true) const;
+    double calculateSusceptibility(bool normalized = true) const;
+
+    // Data Output
+    void saveResults(const std::string& filename) const;
+    void saveEvolution(const std::string& filename) const;
+
+private:
+    // Internal State
+    int latticeSize;
+    double temp; // Units here are [T] = J/kB
+    std::mt19937 rng;
+    std::vector<std::vector<int>> lattice;
+    std::vector<std::vector<double>> evolutionData;
+
+    std::vector<double> boltzmannFactors;
+
+    int n_steps = 0;
+    double sumEnergy = 0.0;
+    double sumSquaredEnergy = 0.0;
+    double sumMagnetization = 0.0;
+    double sumSquaredMagnetization = 0.0;
+
+    // Internal Helpers
+    int getNeighborIndex(int i, int j, int dir, int offset = 1) const;
+    int getNeighborSum(int i, int j) const;
+    double getBoltzmannFactor(int neighborSum, int spin) const;
+    void monteCarloStep();
+    void updateAverages();
+    void recordEvolution();
+};
+#endif
+
@@ -0,0 +1,46 @@
+#include <vector>
+#include <string>
+#include <random>
+#include "ising.hpp"
+
+#ifndef SCHEDULING_HPP
+#define SCHEDULING_HPP
+
+class Scheduler {
+public:
+    // Construction & Initialization
+    Scheduler(int numThreads, std::string outPutDir, bool recordEvolution = true);
+    void addRun(int size, double temperature, int cycles, bool randomInit, int N_duplicates = 1);
+
+    // Execution Control
+    void start();
+
+    // Stochastic Configuration
+    void setMasterSeed(int seed) { masterSeed = seed; seedGen.seed(masterSeed); }
+
+
+private:
+    // Internal State
+    int numThreads;
+    std::string outPutDir;
+    bool recordEvolution;
+
+    int masterSeed = 42;
+    std::mt19937 seedGen;
+
+    struct RunConfig {
+        int size;
+        double temperature;
+        int cycles;
+        bool randomInit;
+        int seed;
+        int duplicateID;
+    };
+
+    std::vector<RunConfig> runQueue;
+
+    // Internal Helpers
+    std::string getOutputFileName(const RunConfig& config, bool isEvolution) const;
+    void executeRun(const RunConfig& config);
+};
+#endif
@@ -0,0 +1,26 @@
+#include "scheduling.hpp"
+
+int main() {
+    // Example usage of Scheduler
+    Scheduler scheduler(12, "system_evolutions", true);
+
+    // Add runs with different parameters
+    scheduler.addRun(20, 1.0, 2000, true, 25);
+    scheduler.addRun(20, 2.4, 2000, true, 25);
+    scheduler.addRun(20, 1.0, 2000, false, 25);
+    scheduler.addRun(20, 2.4, 2000, false, 25);
+
+    // Start the scheduled runs
+    scheduler.start();
+
+    Scheduler scheduler2(12, "phase_transition", false);
+    scheduler2.setMasterSeed(69); // Different master seed for different runs
+    for (int L : {40, 60, 80, 100}) {
+        for (double T = 2.1; T <= 2.4; T += 0.0005) {
+            scheduler2.addRun(L, T, 50000, false);
+        }
+    }
+    scheduler2.start();
+
+    return 0;
+}
@@ -0,0 +1,34 @@
+import matplotlib.pyplot as plt
+
+
+def get_rc_params():
+    colors = ["FF220C", "70D6FF", "8AAA79", "666370", "1C1F33"]
+    rcParams = plt.rcParams
+    # Use LaTeX for rendering
+    # Setup fonts
+    rcParams["text.usetex"] = True
+    rcParams["font.family"] = "serif"
+    rcParams["font.size"] = 10
+    rcParams["axes.labelsize"] = 10
+    rcParams["axes.titlesize"] = 10
+    rcParams["legend.fontsize"] = 8
+    rcParams["xtick.labelsize"] = 8
+    rcParams["ytick.labelsize"] = 8
+    # Figure size and resolution
+    rcParams["figure.figsize"] = (4.5, 3)
+    rcParams["figure.dpi"] = 300
+    # Use colors from the palette
+    rcParams["axes.prop_cycle"] = plt.cycler(color=[f"#{color}" for color in colors])
+    # Grid
+    rcParams["axes.grid"] = True
+    rcParams["grid.alpha"] = 0.5
+    rcParams["grid.linestyle"] = "--"
+    # Point ticks to the inside of the axes
+    rcParams["xtick.direction"] = "in"
+    rcParams["ytick.direction"] = "in"
+    rcParams["xtick.top"] = True
+    rcParams["ytick.right"] = True
+    return rcParams
+
+
+plt.rcParams.update(get_rc_params())
@@ -0,0 +1,193 @@
+#include <vector>
+#include <random>
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include "ising.hpp"
+
+
+IsingModel::IsingModel(int size, double temperature, bool randomInit, int seed)
+    : latticeSize(size), temp(temperature), lattice(size, std::vector<int>(size, 1)), rng(seed) {
+    // Precompute Boltzmann factors for efficiency
+    boltzmannFactors.resize(5);
+    for (int deltaE = -8; deltaE <= 8; deltaE += 4) {
+        boltzmannFactors[(deltaE + 8) / 4] = exp(-deltaE / kB / temp);
+        //std::cout << "Boltzmann factor for ΔE=" << deltaE << ": " << boltzmannFactors[(deltaE + 8) / 4] << "\n";
+    }
+    initializeLattice(randomInit);
+}
+
+void IsingModel::initializeLattice(bool randomInit) {
+    if (!randomInit) {
+        for (int i = 0; i < latticeSize; ++i) {
+            for (int j = 0; j < latticeSize; ++j) {
+                lattice[i][j] = 1; // All spins up
+            }
+        }
+        return;
+    }
+    std::random_device rd;
+    std::mt19937 gen(rd());
+    std::uniform_int_distribution<> dis(0, 1);
+
+    for (int i = 0; i < latticeSize; ++i) {
+        for (int j = 0; j < latticeSize; ++j) {
+            lattice[i][j] = dis(gen) == 0 ? -1 : 1; // Randomly assign -1 or +1
+        }
+    }
+}
+
+int IsingModel::getNeighborIndex(int i, int j, int dir, int offset) const {
+    switch (dir) {
+        case 0: return (i + offset) % latticeSize; // Down
+        case 1: return (j + offset) % latticeSize; // Right
+        case 2: return (i - offset + latticeSize) % latticeSize; // Up
+        case 3: return (j - offset + latticeSize) % latticeSize; // Left
+        default: return -1; // Invalid direction
+    }
+}
+int IsingModel::getNeighborSum(int i, int j) const {
+    return lattice[getNeighborIndex(i, j, 0)][j] + // Down
+           lattice[i][getNeighborIndex(i, j, 1)] + // Right
+           lattice[getNeighborIndex(i, j, 2)][j] + // Up
+           lattice[i][getNeighborIndex(i, j, 3)];   // Left
+}
+
+double IsingModel::getBoltzmannFactor(int neighborSum, int spin) const {
+    int deltaE = 2 * spin * neighborSum;
+    return boltzmannFactors[(deltaE + 8) / 4];
+}
+
+void IsingModel::monteCarloStep() {
+    std::uniform_int_distribution<> dist(0, latticeSize - 1);
+    int i = dist(rng);
+    int j = dist(rng);
+    int neighborSum = getNeighborSum(i, j);
+    int spin = lattice[i][j];
+    double acceptanceProb = getBoltzmannFactor(neighborSum, spin);
+
+    std::uniform_real_distribution<> probDist(0.0, 1.0);
+    if (probDist(rng) < acceptanceProb) {
+        lattice[i][j] *= -1; // Flip spin
+    }
+}
+
+void IsingModel::monteCarloCycle() {
+    for (int step = 0; step < latticeSize * latticeSize; ++step) {
+        monteCarloStep();
+    }
+    updateAverages();
+    recordEvolution();
+}
+
+void IsingModel::run(int cycles) {
+    for (int cycle = 0; cycle < cycles; ++cycle) {
+        monteCarloCycle();
+    }
+    //std::cout << "Accepted Moves: " << acceptedMoves << "\n";
+    //std::cout << "Declined Moves: " << declinedMoves << "\n";
+}
+
+void IsingModel::updateAverages() {
+    double mag = currentMagnetization(false);
+    double en = currentEnergy(false);
+
+    n_steps++;
+    sumMagnetization += mag;
+    sumSquaredMagnetization += mag * mag;
+    sumEnergy += en;
+    sumSquaredEnergy += en * en;
+}
+
+double IsingModel::currentMagnetization(bool normalized) const {
+    double totalMagnetization = 0.0;
+    for (int i = 0; i < latticeSize; ++i) {
+        for (int j = 0; j < latticeSize; ++j) {
+            totalMagnetization += lattice[i][j];
+        }
+    }
+    return normalized? abs(totalMagnetization) / (latticeSize * latticeSize) : abs(totalMagnetization);
+}
+
+double IsingModel::currentEnergy(bool normalized) const {
+    double totalEnergy = 0.0;
+    for (int i = 0; i < latticeSize; ++i) {
+        for (int j = 0; j < latticeSize; ++j) {
+            int spin = lattice[i][j];
+            int neighborSum = getNeighborSum(i, j);
+            totalEnergy -= spin * neighborSum / 2.0; // Each pair counted twice
+        }
+    }
+    return normalized ? totalEnergy / (latticeSize * latticeSize) : totalEnergy;
+}
+
+double IsingModel::calculateMagnetization(bool normalized) const {
+    double avgMagnetization = sumMagnetization / n_steps;
+    return normalized ? avgMagnetization / (latticeSize * latticeSize) : avgMagnetization;
+}
+
+double IsingModel::calculateEnergy(bool normalized) const {
+    double avgEnergy = sumEnergy / n_steps;
+    return normalized ? avgEnergy / (latticeSize * latticeSize) : avgEnergy;
+}
+
+double IsingModel::calculateMagnetizationVariance(bool normalized) const {
+    double variance = (sumSquaredMagnetization / n_steps) - sumMagnetization * sumMagnetization / (n_steps * n_steps);
+    return normalized ? variance / (latticeSize * latticeSize * latticeSize * latticeSize) : variance;
+}
+
+double IsingModel::calculateEnergyVariance(bool normalized) const {
+    double variance = (sumSquaredEnergy / n_steps) - sumEnergy * sumEnergy / (n_steps * n_steps);
+    return normalized ? variance / (latticeSize * latticeSize * latticeSize * latticeSize) : variance;
+}
+
+double IsingModel::calculateHeatCapacity(bool normalized) const {
+    double meanEnergyVar = calculateEnergyVariance(true);
+    double heatCapacity = meanEnergyVar / (kB * temp * temp);
+    return normalized ? heatCapacity / (latticeSize * latticeSize) : heatCapacity;
+}
+
+double IsingModel::calculateSusceptibility(bool normalized) const {
+    double meanMagVar = calculateMagnetizationVariance(true);
+    double susceptibility = meanMagVar / (kB * temp);
+    return normalized ? susceptibility / (latticeSize * latticeSize) : susceptibility;
+}
+
+void IsingModel::saveResults(const std::string& filename) const {
+    std::ofstream outFile(filename);
+    if (!outFile) {
+        throw std::runtime_error("Could not open file for writing: " + filename);
+    }
+
+    outFile << "Lattice Size: " << latticeSize << "\n";
+    outFile << "Temperature: " << temp << "\n";
+    outFile << "Number of Steps: " << n_steps << "\n";
+    outFile << "Average Energy: " << calculateEnergy() << "\n";
+    outFile << "Average Magnetization: " << calculateMagnetization() << "\n";
+    outFile << "Heat Capacity: " << calculateHeatCapacity() << "\n";
+    outFile << "Susceptibility: " << calculateSusceptibility() << "\n";
+
+    outFile.close();
+}
+
+void IsingModel::recordEvolution() {
+    double en = currentEnergy();
+    double mag = currentMagnetization();
+    double En = calculateEnergy();
+    double Mag = calculateMagnetization();
+    evolutionData.push_back({static_cast<double>(n_steps), en, mag, En, Mag});
+}
+
+void IsingModel::saveEvolution(const std::string& filename) const {
+    std::ofstream outFile(filename);
+    if (!outFile) {
+        throw std::runtime_error("Could not open file for writing: " + filename);
+    }
+
+    outFile << "Step,CurrentEnergy,CurrentMagnetization,AvgEnergy,AvgMagnetization\n";
+    for (const auto& record : evolutionData) {
+        outFile << record[0] << "," << record[1] << "," << record[2] << "," << record[3] << "," << record[4] << "\n";
+    }
+
+    outFile.close();
+}
@@ -0,0 +1,58 @@
+#include "scheduling.hpp"
+#include "omp.h"
+#include <iostream>
+#include <filesystem>
+
+Scheduler::Scheduler(int numThreads, std::string outPutDir, bool recordEvolution) : numThreads(numThreads), outPutDir(outPutDir), recordEvolution(recordEvolution) {
+    // Create output directory if it doesn't exist
+    std::filesystem::create_directories(outPutDir);
+}
+
+void Scheduler::addRun(int size, double temperature, int cycles, bool randomInit, int N_duplicates) {
+    std::uniform_int_distribution<> dist(0, 1000000);
+
+    for (int i = 0; i < N_duplicates; ++i) {
+        runQueue.push_back({size, temperature, cycles, randomInit, dist(seedGen), i});
+    }
+}
+
+void Scheduler::executeRun(const RunConfig& config) {
+    IsingModel model(config.size, config.temperature, config.randomInit, config.seed);
+    model.run(config.cycles);
+
+    double finalEnergy = model.calculateEnergy();
+    double finalMagnetization = model.calculateMagnetization();
+
+    std::string outputFileName = getOutputFileName(config, false);
+    model.saveResults(outputFileName);
+    if (recordEvolution) {
+        std::string evolutionFileName = getOutputFileName(config, true);
+        model.saveEvolution(evolutionFileName);
+    }
+
+    // Output results (could be saved to file or processed further)
+    #pragma omp critical
+    {
+        std::cout << "Run (Size: " << config.size
+                  << ", Temp: " << config.temperature
+                  << ", Cycles: " << config.cycles
+                  << ", Seed: " << config.seed << ") -> "
+                  << "Final Energy: " << finalEnergy
+                  << ", Final Magnetization: " << finalMagnetization << std::endl;
+    }
+}
+
+std::string Scheduler::getOutputFileName(const RunConfig& config, bool isEvolution) const {
+    return outPutDir + "/ising_L" + std::to_string(config.size) +
+           (config.randomInit ? "_R" : "_O") +
+           "_T" + std::to_string(config.temperature) +
+           "_C" + std::to_string(config.cycles) +
+           "_D" + std::to_string(config.duplicateID) + (isEvolution ? "_evolution.txt" : ".txt");
+}
+
+void Scheduler::start() {
+    #pragma omp parallel for num_threads(numThreads)
+    for (size_t i = 0; i < runQueue.size(); ++i) {
+        executeRun(runQueue[i]);
+    }
+}