Merge branch 'ws5' of https://github.com/JoHaHu/MolSim into ws5

JoHaHu · Jul 18, 2024 · c54166a · c54166a
2 parents d75d8f3 + ced209d
commit c54166a
Show file tree

Hide file tree

Showing 25 changed files with 6,138 additions and 3,490 deletions.
diff --git a/input/eingabe-membrane-ws5-step1.xml b/input/eingabe-membrane-ws5-step1.xml
@@ -0,0 +1,44 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<scenario xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+          xs:noNamespaceSchemaLocation="../src/lib/simulator/io/xml_reader/SimulationInputSchema.xsd">
+    <header base_name="Membrane"
+            output_file="../../output/membrane"
+            output_frequency="0"
+            t_end="150" seed="1337"
+            delta_t="0.01"
+            dimensions="3"
+            vectorized="false"
+            parallelized="true"
+    />
+    <checkpoints path="../../output/checkpoint.xml"/>
+    <membrane fzup="0.8" k="300" rzero="2.2"/>
+    <container>
+        <linked_cells cutoff_radius="4">
+            <domain_size x="148" y="148" z="148"/>
+            <boundary_conditions>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+            </boundary_conditions>
+        </linked_cells>
+        <!--        <vector/>-->
+    </container>
+    <forces>
+        <lennard_jones>
+            <gravity>-0.001</gravity>
+            <particleTypes>
+                <particleType id="0" sigma="1" epsilon="1" mass="1"/>
+            </particleTypes>
+            <particles>
+                <membrane spacing="2.2" particleTypeId="0">
+                    <coordinate x="15" y="15" z="15"/>
+                    <dimensions x="50" y="20" z="1"/>
+                    <velocity x="0" y="0" z="0"/>
+                </membrane>
+            </particles>
+        </lennard_jones>
+    </forces>
+</scenario>
diff --git a/input/eingabe-membrane-ws5-step2.xml b/input/eingabe-membrane-ws5-step2.xml
@@ -0,0 +1,39 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<scenario xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+          xs:noNamespaceSchemaLocation="../src/lib/simulator/io/xml_reader/SimulationInputSchema.xsd">
+    <header base_name="Membrane"
+            output_file="../../output/membrane2"
+            output_frequency="0"
+            t_end="350" seed="1337"
+            delta_t="0.01"
+            dimensions="3"
+            vectorized="false"
+            parallelized="true"
+    />
+    <checkpoints>
+        <checkpoint>../../output/checkpoint.xml</checkpoint>
+    </checkpoints>
+    <membrane fzup="0" k="300" rzero="2.2"/>
+    <container>
+        <linked_cells cutoff_radius="4">
+            <domain_size x="148" y="148" z="148"/>
+            <boundary_conditions>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+                <boundary_condition type="reflecting"/>
+            </boundary_conditions>
+        </linked_cells>
+        <!--        <vector/>-->
+    </container>
+    <forces>
+        <lennard_jones>
+            <gravity>-0.001</gravity>
+            <particleTypes>
+                <particleType id="0" sigma="1" epsilon="1" mass="1"/>
+            </particleTypes>
+        </lennard_jones>
+    </forces>
+</scenario>
diff --git a/src/bin/MolSim.cpp b/src/bin/MolSim.cpp
@@ -14,7 +14,7 @@
 #include "utils/LoggerManager.h"
 
 template<const size_t DIMENSIONS>
-auto setup(std::shared_ptr<config::Config> config) -> auto{
+auto setup(std::shared_ptr<config::Config> config) -> auto {
   auto particle_loader = simulator::io::ParticleGenerator<DIMENSIONS>(config);
 
   auto particles_vector = particle_loader.load_particles();
@@ -70,12 +70,18 @@ auto setup(std::shared_ptr<config::Config> config) -> auto{
   }
 
   pc->refresh();
+  auto membrane_force = std::optional<simulator::physics::MembraneForce>();
+
+  if (config->k != 0) {
+    membrane_force = std::optional(simulator::physics::MembraneForce(config->k, config->rzero));
+  }
 
   auto simulator = simulator::Simulator<DIMENSIONS>(
       std::move(pc),
       std::move(plotter),
       config,
-      checkpointer);
+      checkpointer,
+      std::move(membrane_force));
   return simulator;
 }
 

diff --git a/src/lib/Particle.h b/src/lib/Particle.h
@@ -1,10 +1,12 @@
 #pragma once
 
 #include "range/v3/algorithm.hpp"
+#include "range/v3/view/iota.hpp"
 #include "range/v3/view/zip.hpp"
 #include <array>
 #include <execution>
 #include <experimental/simd>
+#include <range/v3/range/conversion.hpp>
 #include <vector>
 
 /**
@@ -50,14 +52,16 @@ class Particle {
      */
   uint8_t fixed{};
 
+  uint8_t is_membrane{};
+
   Particle(const std::array<double, DIMENSIONS> &position,
            const std::array<double, DIMENSIONS> &velocity,
            const std::array<double, DIMENSIONS> &force,
            const std::array<double, DIMENSIONS> &old_force,
            double mass,
            int type,
-           int fixed) : position(position), velocity(velocity), force(force), old_force(old_force), mass(mass),
-                        type(type), fixed(fixed) {}
+           int fixed, int is_membrane) : position(position), velocity(velocity), force(force), old_force(old_force), mass(mass),
+                                         type(type), fixed(fixed), is_membrane(is_membrane) {}
 
   Particle(
       // for visualization, we need always 3 coordinates
@@ -71,6 +75,13 @@ class Particle {
   }
 };
 
+struct MembranePair {
+ public:
+  size_t first;
+  size_t second;
+  bool diagonal;
+};
+
 /***
  * A structure of Arrays for the particles
  *
@@ -94,6 +105,12 @@ class Particles {
   std::vector<size_t> cell{};
   std::vector<size_t> block{};
   std::vector<size_t> color{};
+  std::vector<uint8_t> membrane{};
+  std::vector<MembranePair> membrane_pairs{};
+  /**
+   * Used to detect swapped particles
+   * */
+  std::vector<size_t> swap_index{};
 
 #if DEBUG
   /**
@@ -124,22 +141,29 @@ class Particles {
     cell.emplace_back(0);
     block.emplace_back(0);
     color.emplace_back(0);
+    membrane.emplace_back(p.is_membrane);
+    swap_index.emplace_back(0);
 #if DEBUG
     ids.emplace_back(size);
 #endif
     size++;
   }
 
+  auto insert_membrane_pair(MembranePair p) -> void {
+    membrane_pairs.emplace_back(p);
+  }
+
   size_t size{};
 
   template<typename Callable>
   void sort(Callable c) {
-    // recaulcates the cell
+    // recalculates the cell
     for (size_t i = 0; i < size; ++i) {
       auto [_color, _block, _cell] = c(i);
       color[i] = _color;
       block[i] = _block;
       cell[i] = _cell;
+      swap_index[i] = i;
     }
     // This is ugly but i haven't found a better way
     std::apply([&](auto &...ps) {
@@ -154,7 +178,7 @@ class Particles {
                 vs...,
                 fs...,
                 ofs...,
-                mass, type, active, fixed
+                mass, type, active, fixed, membrane, swap_index
 #if DEBUG
                 ,
                 ids
@@ -176,6 +200,11 @@ class Particles {
                         velocities);
     },
                positions);
+
+    for (auto &pair : membrane_pairs) {
+      pair.first = swap_index[pair.first];
+      pair.second = swap_index[pair.second];
+    }
   }
 
   void
@@ -189,4 +218,20 @@ class Particles {
       }
     }
   }
+
+  void apply_membrane(std::function<void(size_t, size_t, bool)> const &f) {
+    for (auto &pair : membrane_pairs) {
+      f(pair.first, pair.second, pair.diagonal);
+    }
+  }
+
+ private:
+  template<typename T>
+  void reorder_vec(std::vector<T> &input, const std::vector<size_t> &permutation) {
+    auto result = std::vector<T>(permutation.size());
+    for (size_t i = 0; i < permutation.size(); ++i) {
+      result[permutation[i]] = input[i];
+    }
+    std::swap(input, result);
+  }
 };
diff --git a/src/lib/config/Config.h b/src/lib/config/Config.h
@@ -166,6 +166,10 @@ class Config {
    */
   std::vector<Cuboid> cuboids;
 
+  /**
+   * a vector that can store multiple cuboids for simulation defined in the Cuboid class
+   */
+  std::vector<Cuboid> membrane;
   /**
    * a vector that can store multiple discs for simulation defined in the Disc class
    */
@@ -212,6 +216,10 @@ class Config {
    */
   std::string input_filename;
 
+  double k = 0;
+  double rzero = 0;
+  double fzup = 0;
+
   /**
    * prints the help message
    */

diff --git a/src/lib/config/Cuboid.h b/src/lib/config/Cuboid.h
@@ -28,15 +28,15 @@ class Cuboid {
          std::vector<double> &velocity,
          double mass, double spacing, int type)
       : coordinates(coordinates), particles(particles), velocity(velocity), mass(mass), spacing(spacing),
-        type(type), fixed(0){};
+        type(type), fixed(0) {};
 
   // Constructor for fixed cuboids (additional parameter)
   Cuboid(std::vector<double> &coordinates,
          std::vector<double> &particles,
          std::vector<double> &velocity,
          double mass, double spacing, int type, int fixed)
       : coordinates(coordinates), particles(particles), velocity(velocity), mass(mass), spacing(spacing),
-        type(type), fixed(fixed){};
+        type(type), fixed(fixed) {};
 
   std::string toString() const {
     std::string output;

diff --git a/src/lib/container/Container.h b/src/lib/container/Container.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include "Particle.h"
+#include "simulator/physics/MembraneForce.h"
 namespace container {
 
 template<const size_t DIMENSIONS>
@@ -13,6 +14,13 @@ class Container {
   * @param f Function to apply to each particle.
   */
   virtual void linear(std::function<void(Particles<DIMENSIONS> &, size_t)>) = 0;
+
+  /**
+  * @brief Applies a function to each particle in the container.
+  * @param f Function to apply to each particle.
+  */
+  virtual void membrane(simulator::physics::MembraneForce &) = 0;
+
   /**
   * @brief Applies a force to each pair of particles in the container.
   *
@@ -47,6 +55,13 @@ class Container {
   * @param p Particle to insert.
   */
   virtual void insert(Particle<DIMENSIONS> p) = 0;
+
+  /**
+  * @brief Inserts a membrane pair into the container.
+  *
+  * @param p Particle to insert.
+  */
+  virtual void insert_membrane_pair(MembranePair p) = 0;
 };
 
 }// namespace container