Fix issue in calculation of heat cap and susc.. Add performance test and analytical comp

.gitignore

@@ -297,4 +297,9 @@ src/project3/python/plots/
 # Project 4 executables
 src/project4/phase_transition/
 src/project4/system_evolutions/
-src/project4/main
+src/project4/results/
+src/project4/analytical_comparison/
+src/project4/main
+src/project4/analytical
+src/project4/single_threaded_performance
+src/project4/multi_threaded_performance

src/project4/analytical.cpp

@@ -0,0 +1,10 @@
+#include "scheduling.hpp"
+
+int main() {
+    // Comparison runs against analytical results
+    Scheduler scheduler3(12, "analytical_comparison", true);
+    scheduler3.setMasterSeed(12345);
+    scheduler3.addRun(2, 1.0, 10000, true, 10);
+    scheduler3.start();
+
+}

src/project4/generate_results.sh

@@ -0,0 +1,39 @@
+#!/bin/bash
+
+# Compile the C++ programs
+rm -f analytical main single_threaded_performance multi_threaded_performance
+echo "Compiling C++ programs..."
+g++ -I include -fopenmp -O3 src/* analytical.cpp -o analytical
+g++ -I include -fopenmp -O3 src/* main.cpp -o main
+g++ -I include -fopenmp -O3 src/* performance_test.cpp -o multi_threaded_performance
+g++ -I include -O3 src/* performance_test.cpp -o single_threaded_performance
+echo "Compilation done."
+
+# Deleting old results
+rm -rf analytical_comparison
+rm -rf phase_transition
+rm -rf system_evolutions
+rm -rf results
+
+# Creating new results
+echo "Generating new results..."
+./analytical
+echo "Analytical results generated."
+./main
+echo "Main results generated."
+{ time ./single_threaded_performance; } 2> single_threaded_time.txt
+{ time ./multi_threaded_performance; } 2> multi_threaded_time.txt
+echo "Performance test results generated."
+echo "Simulation results generated."
+
+# Backing up results
+echo "Backing up results..."
+timestamp=$(date +"%Y%m%d_%H%M%S")
+backup_dir="/home/lars/Programming/Backup_$timestamp"
+mkdir -p "$backup_dir"
+for dir in analytical_comparison phase_transition system_evolutions results; do
+    if [ -d "$dir" ]; then
+        tar -czf "$backup_dir/${dir}_$timestamp.tar.gz" "$dir"
+    fi
+done
+echo "Backup completed at $backup_dir."

src/project4/include/ising.hpp

@@ -24,12 +24,12 @@ public:
     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;
+    // The following are always normalized
+    double calculateEnergyVariance() const; 
+    double calculateMagnetizationVariance() const;
 
+    double calculateHeatCapacity() const;
+    double calculateSusceptibility() const;
     // Data Output
     void saveResults(const std::string& filename) const;
     void saveEvolution(const std::string& filename) const;

src/project4/main.cpp

@@ -16,7 +16,7 @@ int main() {
     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) {
+        for (double T = 2.1; T <= 2.4; T += 0.001) {
             scheduler2.addRun(L, T, 50000, false);
         }
     }

src/project4/performance_test.cpp

@@ -0,0 +1,9 @@
+#include "scheduling.hpp"
+
+int main() {
+    // Just a dummy run
+    Scheduler scheduler(12, "results", true);
+    scheduler.setMasterSeed(2025);
+    scheduler.addRun(100, 2.3, 100000, true, 100);
+    scheduler.start();
+}

src/project4/python/01_systemEvolution.ipynb

@@ -1459,7 +1459,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.13.9"
+   "version": "3.13.7"
   }
  },
  "nbformat": 4,

src/project4/src/ising.cpp

@@ -89,8 +89,8 @@ void IsingModel::run(int cycles) {
 }
 
 void IsingModel::updateAverages() {
-    double mag = currentMagnetization(false);
-    double en = currentEnergy(false);
+    double mag = currentMagnetization(true);
+    double en = currentEnergy(true);
 
     n_steps++;
     sumMagnetization += mag;
@@ -106,7 +106,7 @@ double IsingModel::currentMagnetization(bool normalized) const {
             totalMagnetization += lattice[i][j];
         }
     }
-    return normalized? abs(totalMagnetization) / (latticeSize * latticeSize) : abs(totalMagnetization);
+    return normalized? static_cast<double>(abs(totalMagnetization)) / (latticeSize * latticeSize) : static_cast<double>(abs(totalMagnetization));
 }
 
 double IsingModel::currentEnergy(bool normalized) const {
@@ -118,39 +118,41 @@ double IsingModel::currentEnergy(bool normalized) const {
             totalEnergy -= spin * neighborSum / 2.0; // Each pair counted twice
         }
     }
-    return normalized ? totalEnergy / (latticeSize * latticeSize) : totalEnergy;
+    return normalized ? static_cast<double>(totalEnergy) / (latticeSize * latticeSize) : static_cast<double>(totalEnergy);
 }
 
 double IsingModel::calculateMagnetization(bool normalized) const {
     double avgMagnetization = sumMagnetization / n_steps;
-    return normalized ? avgMagnetization / (latticeSize * latticeSize) : avgMagnetization;
+    return normalized ? avgMagnetization : avgMagnetization * (latticeSize * latticeSize);
 }
 
 double IsingModel::calculateEnergy(bool normalized) const {
     double avgEnergy = sumEnergy / n_steps;
-    return normalized ? avgEnergy / (latticeSize * latticeSize) : avgEnergy;
+    return normalized ? avgEnergy : avgEnergy * (latticeSize * latticeSize);
 }
 
-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::calculateMagnetizationVariance() const {
+    double normalization = 1.0 / (n_steps);
+    double variance = (sumSquaredMagnetization * normalization) - (sumMagnetization * normalization) * (sumMagnetization * normalization);
+    return 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::calculateEnergyVariance() const {
+    double normalization = 1.0 / (n_steps);
+    double variance = (sumSquaredEnergy * normalization) - (sumEnergy * normalization) * (sumEnergy * normalization);
+    return 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::calculateHeatCapacity() const {
+    double meanEnergyVar = calculateEnergyVariance();
+    double heatCapacity = meanEnergyVar / (kB * temp * temp) * (latticeSize * latticeSize);
+    return heatCapacity;
 }
 
-double IsingModel::calculateSusceptibility(bool normalized) const {
-    double meanMagVar = calculateMagnetizationVariance(true);
-    double susceptibility = meanMagVar / (kB * temp);
-    return normalized ? susceptibility / (latticeSize * latticeSize) : susceptibility;
+double IsingModel::calculateSusceptibility() const {
+    double meanMagVar = calculateMagnetizationVariance();
+    double susceptibility = meanMagVar / (kB * temp) * (latticeSize * latticeSize);
+    return susceptibility;
 }
 
 void IsingModel::saveResults(const std::string& filename) const {

src/project4/src/scheduling.cpp

@@ -51,8 +51,17 @@ std::string Scheduler::getOutputFileName(const RunConfig& config, bool isEvoluti
 }
 
 void Scheduler::start() {
-    #pragma omp parallel for num_threads(numThreads)
+#ifdef _OPENMP
+    std::cout << "Running with OpenMP using " << numThreads << " threads.\n";
+    #pragma omp parallel for schedule(dynamic, 4) num_threads(numThreads)
     for (size_t i = 0; i < runQueue.size(); ++i) {
         executeRun(runQueue[i]);
     }
-}
+#else
+    // Fallback: single-threaded execution
+    std::cout << "Running without OpenMP (single-threaded).\n";
+    for (size_t i = 0; i < runQueue.size(); ++i) {
+        executeRun(runQueue[i]);
+    }
+#endif
+}