Assignment3

CMakeLists.txt

@@ -66,6 +66,8 @@ file(GLOB_RECURSE TEST_SOURCES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/objects/*.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/calculators/*.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/boundaries/*.cpp"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/utils/*.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/utils/*.cpp"
 )
 
 file(GLOB_RECURSE BENCHMARK_SOURCES
@@ -77,6 +79,8 @@ file(GLOB_RECURSE BENCHMARK_SOURCES
     "${CMAKE_CURRENT_SOURCE_DIR}/src/io/in/*.h"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/calculators/*.cpp"
     "${CMAKE_CURRENT_SOURCE_DIR}/src/calculators/*.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/utils/*.h"
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/utils/*.cpp"
 )
 
 #------------------------------------------LOGGING---------------------------------------------------------

src/Main.cpp

@@ -10,6 +10,8 @@
 #include "MolSim.h"
 #include <iostream>
 #include <memory>
+#include <map>
+#include "boundaries/BoundaryController.h"
 
 /**
  * @brief The main function of the program.
@@ -37,18 +39,19 @@ int main(int argc, char *argsv[]) {
     double end_time;
     std::unique_ptr<outputWriters::FileWriter> outputWriter;
     std::unique_ptr<calculators::Calculator> calculator;
+    std::map<boundaries::BoundaryDirection, boundaries::BoundaryType> boundaryMap;
 
-    if (!MolSim::processArguments(argc, argsv, inputFilePath, delta_t, end_time, outputWriter, calculator)) {
+    if (!MolSim::processArguments(argc, argsv, inputFilePath, delta_t, end_time, outputWriter, calculator, boundaryMap)) {
         return 1;
     }
 
     ParticleContainer particleContainer = FileReader::readFile(inputFilePath);
 
     SPDLOG_INFO("Starting simulation with delta_t: {}, end_time: {}", delta_t, end_time);
 
-    MolSim::performSimulation(particleContainer, delta_t, end_time, outputWriter, calculator);
+    MolSim::performSimulation(particleContainer, delta_t, end_time, outputWriter, calculator, boundaryMap);
 
     SPDLOG_INFO("Simulation completed.");
 
     return 0;
-}
+}

src/MolSim.h

@@ -70,7 +70,8 @@ class MolSim {
     static bool processArguments(int argc, char *argv[], std::string &inputFilePath,
                                  double &delta_t, double &end_time,
                                  std::unique_ptr<outputWriters::FileWriter> &outputWriter,
-                                 std::unique_ptr<calculators::Calculator> &calculator) {
+                                 std::unique_ptr<calculators::Calculator> &calculator,
+                                 std::map<boundaries::BoundaryDirection, boundaries::BoundaryType> &boundaryMap) {
         cxxopts::Options options("MolSim", "Molecular Simulation Program");
 
         options.add_options()
@@ -79,7 +80,8 @@ class MolSim {
                 ("delta_t", "Set delta_t", cxxopts::value<double>()->default_value("0.014"))
                 ("end_time", "Set end_time", cxxopts::value<double>()->default_value("1000"))
                 ("output", "Output writer (vtk or xyz)", cxxopts::value<std::string>())
-                ("calculator", "Calculator (sv, lj or dummy)", cxxopts::value<std::string>());
+                ("calculator", "Calculator (sv, lj or dummy)", cxxopts::value<std::string>())
+                ("boundaries", "Boundary conditions (0000 to 2222)", cxxopts::value<std::string>()->default_value("1111"));
 
         auto result = options.parse(argc, argv);
 
@@ -140,7 +142,46 @@ class MolSim {
             SPDLOG_ERROR("Invalid input; please select a calculator.");
             return false;
         }
-
+        if (result.count("boundaries")) {
+            std::string boundariesArg = result["boundaries"].as<std::string>();
+            if (boundariesArg.size() == 4) {
+                std::array<boundaries::BoundaryDirection, 4> directions = {
+                        boundaries::BoundaryDirection::LEFT,
+                        boundaries::BoundaryDirection::BOTTOM,
+                        boundaries::BoundaryDirection::RIGHT,
+                        boundaries::BoundaryDirection::TOP
+                };
+
+                for (int i = 0; i < 4; ++i) {
+                    switch (boundariesArg[i]) {
+                        case '0':
+                            boundaryMap[directions[i]] = boundaries::BoundaryType::OFF;
+                            break;
+                        case '1':
+                            boundaryMap[directions[i]] = boundaries::BoundaryType::REFLECTING;
+                            break;
+                        case '2':
+                            boundaryMap[directions[i]] = boundaries::BoundaryType::OUTFLOW;
+                            break;
+                        default:
+                            SPDLOG_WARN("Invalid boundary value '{}', defaulting all to REFLECTING", boundariesArg[i]);
+                            boundaryMap[directions[i]] = boundaries::BoundaryType::REFLECTING;
+                    }
+                }
+            } else {
+                SPDLOG_WARN("Invalid boundary string length, defaulting all boundaries to REFLECTING");
+                boundaryMap[boundaries::BoundaryDirection::LEFT] = boundaries::BoundaryType::REFLECTING;
+                boundaryMap[boundaries::BoundaryDirection::BOTTOM] = boundaries::BoundaryType::REFLECTING;
+                boundaryMap[boundaries::BoundaryDirection::RIGHT] = boundaries::BoundaryType::REFLECTING;
+                boundaryMap[boundaries::BoundaryDirection::TOP] = boundaries::BoundaryType::REFLECTING;
+            }
+        } else {
+            SPDLOG_WARN("No boundary conditions specified, defaulting all boundaries to REFLECTING");
+            boundaryMap[boundaries::BoundaryDirection::LEFT] = boundaries::BoundaryType::REFLECTING;
+            boundaryMap[boundaries::BoundaryDirection::BOTTOM] = boundaries::BoundaryType::REFLECTING;
+            boundaryMap[boundaries::BoundaryDirection::RIGHT] = boundaries::BoundaryType::REFLECTING;
+            boundaryMap[boundaries::BoundaryDirection::TOP] = boundaries::BoundaryType::REFLECTING;
+        }
         return true;
     }
 
@@ -158,17 +199,13 @@ class MolSim {
     */
     static void performSimulation(ParticleContainer &particleContainer, double &delta_t, double &end_time,
                                   std::unique_ptr<outputWriters::FileWriter> &outputWriter,
-                                  std::unique_ptr<calculators::Calculator> &calculator) {
+                                  std::unique_ptr<calculators::Calculator> &calculator,
+                                  const std::map<boundaries::BoundaryDirection, boundaries::BoundaryType> &boundaryMap) {
         const std::string &filename = "MD";
 
         double current_time = 0.0; // start_time
         int iteration = 0;
 
-        std::map<boundaries::BoundaryDirection, boundaries::BoundaryType> boundaryMap{};
-//        boundaryMap.emplace(boundaries::BoundaryDirection::BOTTOM, boundaries::BoundaryType::OUTFLOW);
-//        boundaryMap.emplace(boundaries::BoundaryDirection::RIGHT, boundaries::BoundaryType::OUTFLOW);
-//        boundaryMap.emplace(boundaries::BoundaryDirection::BOTTOM, boundaries::BoundaryType::OUTFLOW);
-
         std::array<double, 2> domain = {50, 30.0};
 
         const boundaries::BoundaryController controller{boundaryMap, calculator.get(), domain, 1.0};

src/calculators/Thermostat.cpp

@@ -0,0 +1,77 @@
+//
+// Created by kimj2 on 11.06.2024.
+//
+#include "../objects/ParticleContainer.h"
+#include "Thermostat.h"
+#include "../utils/ArrayUtils.h"
+#include "../utils/MaxwellBoltzmannDistribution.h"
+
+const double COMPARISON_TOLERANCE = 1e-9;
+
+void Thermostat::initializeTemp(ParticleContainer &particleContainer) const {
+    if (particleContainer.hasZeroVelocities()) { // Check if all velocities are zero
+        // Initialize velocities by applying Brownian motion
+        double maxwellBoltzmannFactor;
+
+        for (auto &particle : particleContainer) {
+            maxwellBoltzmannFactor = sqrt(start_temp / particle.getM());
+
+            std::array<double, 3> v = particle.getV();
+            const std::array<double, 3> deltaV = maxwellBoltzmannDistributedVelocity(maxwellBoltzmannFactor, 3);
+
+            for (int i = 0; i < dimension; ++i) {
+                v[i] += deltaV[i];
+            }
+
+            particle.setV(v);
+        }
+    }
+
+    // Velocity scaling to initial temperature
+    setTemp(particleContainer, start_temp);
+}
+
+void Thermostat::setTempDirectly(ParticleContainer &particleContainer) const {
+    setTemp(particleContainer, target_temp);
+}
+
+void Thermostat::setTempGradually(ParticleContainer &particleContainer) const {
+    const double currentTemp = calculateCurrentTemp(particleContainer);
+    const double delta_temp = std::min(std::abs(target_temp - currentTemp), this->max_delta_temp);
+    const double updatedTemp = this->target_temp > currentTemp ? currentTemp + delta_temp : currentTemp - delta_temp;
+    setTemp(particleContainer, updatedTemp);
+}
+
+double Thermostat::calculateCurrentTemp(ParticleContainer &particleContainer) const {
+    double kinEnergy = calculateKinEnergy(particleContainer);
+    return 2.0 / (this->dimension * particleContainer.getSize()) * kinEnergy;
+}
+
+//----------------------------------------Helper functions----------------------------------------------------
+
+double Thermostat::calculateKinEnergy(ParticleContainer &particleContainer) const {
+    double kinEnergy = 0.0;
+    for (auto &particle : particleContainer) {
+        kinEnergy += 0.5 * particle.getM() * ArrayUtils::dotProduct(particle.getV(), particle.getV());
+    }
+    return kinEnergy;
+}
+
+void Thermostat::scaleV(double beta, ParticleContainer &particleContainer) const {
+    for (auto &particle : particleContainer) {
+        std::array<double, 3> v = particle.getV();
+        v = beta * v;
+        particle.setV(v);
+    }
+}
+
+void Thermostat::setTemp(ParticleContainer &particleContainer, double newTemp) const {
+    const double currentTemp = calculateCurrentTemp(particleContainer);
+    double beta;
+    if (std::abs(currentTemp) < COMPARISON_TOLERANCE) { // Check if currentTemp is zero
+        beta = sqrt(newTemp / COMPARISON_TOLERANCE); // or some default value
+    } else {
+        beta = sqrt(newTemp / currentTemp);
+    }
+    scaleV(beta, particleContainer);
+}

src/calculators/Thermostat.h

@@ -0,0 +1,57 @@
+//
+// Created by kimj2 on 11.06.2024.
+//
+
+#pragma once
+
+#include "../objects/ParticleContainer.h"
+
+class Thermostat {
+public:
+    Thermostat() = default;
+
+    Thermostat(double start_temp, double target_temp, int applyFrequency, double max_delta_temp, int dimension) :
+            start_temp(start_temp), target_temp(target_temp), applyFrequency(applyFrequency),
+            max_delta_temp(max_delta_temp), dimension(dimension) {}
+
+    void setTempDirectly(ParticleContainer &particleContainer) const;
+
+    void setTempGradually(ParticleContainer &particleContainer) const;
+
+    void initializeTemp(ParticleContainer &particleContainer) const;
+
+    double calculateCurrentTemp(ParticleContainer &particleContainer) const;
+
+    [[nodiscard]] double getStartTemp() const {
+        return start_temp;
+    }
+
+    [[nodiscard]] double getTargetTemp() const {
+        return target_temp;
+    }
+
+    [[nodiscard]] int getApplyFrequency() const {
+        return applyFrequency;
+    }
+
+    [[nodiscard]] double getMaxDeltaTemp() const {
+        return max_delta_temp;
+    }
+
+    [[nodiscard]] int getDimension() const {
+        return dimension;
+    }
+
+private:
+    double start_temp;
+    double target_temp;
+    int applyFrequency;
+    double max_delta_temp;
+    int dimension;
+
+    double calculateKinEnergy(ParticleContainer &particleContainer) const;
+
+    void scaleV(double beta, ParticleContainer &particleContainer) const;
+
+    void setTemp(ParticleContainer &particleContainer, double newTemp) const;
+};

src/objects/Particle.cpp

@@ -10,6 +10,8 @@
 #include <iostream>
 #include "../utils/ArrayUtils.h"
 
+const double COMPARISON_TOLERANCE = 1e-6;
+
 Particle::Particle(int type_arg) {
     type = type_arg;
     f = {0., 0., 0.};
@@ -90,6 +92,12 @@ bool Particle::operator==(const Particle &other) const {
            (volume == other.volume) and (old_f == other.old_f);
 }
 
+bool Particle::hasZeroVelocities() const {
+    const auto &velocities = this->getV();
+    return std::abs(velocities[0]) < COMPARISON_TOLERANCE && std::abs(velocities[1]) < COMPARISON_TOLERANCE &&
+           std::abs(velocities[2]) < COMPARISON_TOLERANCE;
+}
+
 std::ostream &operator<<(std::ostream &stream, Particle &p) {
     stream << p.toString();
     return stream;

src/objects/Particle.h

@@ -91,6 +91,8 @@ class Particle {
     bool operator==(const Particle &other) const;
 
     std::string toString() const;
+
+    bool hasZeroVelocities() const;
 };
 
 std::ostream &operator<<(std::ostream &stream, Particle &p);

src/objects/ParticleContainer.cpp

@@ -88,4 +88,13 @@ size_t ParticleContainer::getSize() const {
 
 std::vector<Particle> ParticleContainer::getParticles() const {
     return particles;
-};
+}
+
+bool ParticleContainer::hasZeroVelocities() const {
+    for (const auto& particle : this->particles) {
+        if (!particle.hasZeroVelocities()) {
+            return false;
+        }
+    }
+    return true;
+}

src/objects/ParticleContainer.h

@@ -87,4 +87,6 @@ class ParticleContainer {
     size_t getSize() const;
 
     std::vector<Particle> getParticles() const;
+
+    bool hasZeroVelocities() const;
 };

src/utils/ArrayUtils.h

@@ -164,6 +164,19 @@ inline Container elementWiseScalarOp(const Scalar &lhs, const Container &rhs,
   SPDLOG_TRACE("Element wise scalar operation: {} {} = {}", lhs, to_string(rhs), to_string(ret));
   return ret;
 }
+/**
+ * Calculates the dot product of two containers.
+ * @tparam Container
+ * @param lhs
+ * @param rhs
+ * @return result of the dot product as double value.
+ */
+template <class Container>
+double dotProduct(const Container &lhs, const Container &rhs) {
+  double result = std::inner_product(lhs.begin(), lhs.end(), rhs.begin(), 0.0);
+  SPDLOG_TRACE("Dot product operation: {} {} = {}", to_string(lhs), to_string(rhs), result);
+  return result;
+}
 
 /**
  * Calculates the L2 norm for a given container.

src/utils/MaxwellBoltzmannDistribution.cpp

@@ -0,0 +1,18 @@
+//
+// Created by kimj2 on 12.06.2024.
+//
+#include "MaxwellBoltzmannDistribution.h"
+#include <random>
+#include <array>
+#include <spdlog/spdlog.h>
+
+std::array<double, 3> maxwellBoltzmannDistributedVelocity(double averageVelocity, size_t dimensions) {
+    static std::default_random_engine randomEngine(42);
+    std::normal_distribution<double> normalDistribution{0, 1};
+    std::array<double, 3> randomVelocity{};
+    for (size_t i = 0; i < dimensions; ++i) {
+        randomVelocity[i] = averageVelocity * normalDistribution(randomEngine);
+    }
+    SPDLOG_TRACE("Generated random velocity: [{}, {}, {}]", randomVelocity[0], randomVelocity[1], randomVelocity[2]);
+    return randomVelocity;
+}