WIP: mVLSI paper

examples/Hybrid/Network1a.JSON

@@ -172,7 +172,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [1.0, 1.0, 1.0, 1.0],
                 "beta": [1.0, 1.0, 1.0, 1.0],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-1,
                 "tau": 0.55,

examples/Hybrid/Network2a.JSON

@@ -143,7 +143,8 @@
     ],
     "activeFixture": 0,
     "updateScheme": {
-        "scheme": "Naive"
+        "scheme": "LinearDecoupling",
+        "theta": 1
     },
     "settings": {
         "simulators": [
@@ -154,8 +155,8 @@
                 "charPhysLength": 1e-4,
                 "charPhysVelocity": 1e-1,
                 "alpha": [0.5, 0.005, 0.005],
-                "beta": [1.0, 0.025, 0.025],
-                "theta": 10,
+                "beta": [0.0, 0.0, 0.0],
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-1,
                 "tau": 0.55,
@@ -199,9 +200,9 @@
                 "stlFile": "../examples/STL/T-shape.stl",
                 "charPhysLength": 1e-4,
                 "charPhysVelocity": 1e-1,
-                "alpha": [0.005, 0.005, 0.05],
-                "beta": [0.025, 0.025, 1.0],
-                "theta": 10,
+                "alpha": [0.0, 0.0, 0.05],
+                "beta": [0.025, 0.025, 0.0],
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-1,
                 "tau": 0.55,

examples/Hybrid/Network3a.JSON

@@ -214,7 +214,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [1.0, 0.03, 0.02],
                 "beta": [0.0, 0.0, 0.0],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-1,
                 "tau": 0.55,
@@ -261,7 +261,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [0.0, 0.0, 0.005, 0.005],
                 "beta": [1.0, 1.0, 0.0, 0.0],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-1,
                 "tau": 0.55,
@@ -318,7 +318,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [0.0, 0.0, 0.1],
                 "beta": [0.05, 0.05, 0.0],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-1,
                 "tau": 0.55,

examples/Hybrid/Network4a.JSON

@@ -213,7 +213,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [0.003, 0.001, 0.001],
                 "beta": [0.015, 0.01, 0.01],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-2,
                 "tau": 0.55,
@@ -259,7 +259,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [0.003, 0.003, 0.003],
                 "beta": [0.015, 0.015, 0.015],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-2,
                 "tau": 0.55,
@@ -306,7 +306,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [0.003, 0.003, 0.003],
                 "beta": [0.015, 0.015, 0.015],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-2,
                 "tau": 0.55,
@@ -353,7 +353,7 @@
                 "charPhysVelocity": 1e-1,
                 "alpha": [0.003, 0.003, 0.003],
                 "beta": [0.015, 0.015, 0.015],
-                "theta": 10,
+                "theta": 1,
                 "resolution": 20,
                 "epsilon": 1e-2,
                 "tau": 0.55,

src/baseSimulator.h

@@ -26,8 +26,10 @@
 
 #include "nodalAnalysis/NodalAnalysis.h"
 
+#include "hybridDynamics/DynamicDamping.h"
 #include "hybridDynamics/Scheme.h"
 #include "hybridDynamics/Naive.h"
+#include "hybridDynamics/LinearDecoupling.h"
 
 #include "architecture/Channel.h"
 #include "architecture/ChannelPosition.h"

src/baseSimulator.hh

@@ -20,8 +20,10 @@
 
 #include "nodalAnalysis/NodalAnalysis.hh"
 
+#include "hybridDynamics/DynamicDamping.hh"
 #include "hybridDynamics/Scheme.hh"
 #include "hybridDynamics/Naive.hh"
+#include "hybridDynamics/LinearDecoupling.hh"
 
 #include "architecture/Channel.hh"
 #include "architecture/ChannelPosition.hh"

src/hybridDynamics/CMakeLists.txt

@@ -1,9 +1,13 @@
 set(SOURCE_LIST
+    DynamicDamping.hh
+    LinearDecoupling.hh
     Scheme.hh
     Naive.hh
 )
 
 set(HEADER_LIST
+    DynamicDamping.h
+    LinearDecoupling.h
     Scheme.h
     Naive.h
 )

src/hybridDynamics/DynamicDamping.h

@@ -0,0 +1,42 @@
+/**
+ * @file Naive.h
+ */
+
+#pragma once
+
+#include <memory>
+#include <unordered_map>
+#include <vector>
+
+namespace mmft{
+
+/**
+ * @brief The Naive Scheme is an update scheme that sets the relaxation factor of an iterative update scheme to a 
+ * given constant for all nodes.
+ */
+template<typename T>
+class DynamicDampingScheme : public Scheme<T> {
+
+private:
+
+    Eigen::VectorXd x_prev;
+    Eigen::VectorXd factor;
+
+public:
+    /**
+     * @brief Constructor of the Naive Scheme with provided constants.
+     * @param[in] modules The map of modules with boundary nodes upon which this scheme acts.
+     * @param[in] alpha The relaxation value for the pressure value update.
+     * @param[in] beta The relaxation value of the flow rate value update.
+     * @param[in] theta The amount of LBM stream and collide cycles between updates for a module.
+     */
+    DynamicDampingScheme(int theta, int size);
+
+    /**
+     * @brief Compute the update value alpha according to alpha = 1/(1 + 10*d2xdt2).
+     */
+    void compute() override;
+
+};
+
+}   // namespace mmft

src/hybridDynamics/DynamicDamping.hh

@@ -0,0 +1,33 @@
+#include "DynamicDamping.h"
+
+namespace mmft {
+
+template<typename T>
+DynamicDampingScheme<T>::DynamicDampingScheme(int theta_, int size_) :
+    Scheme<T>() 
+{ 
+    /*
+    this->theta.try_emplace(0,theta_);
+    for (int i = 0; i < size_; i++) {
+        this->alpha.try_emplace(i, 0.0);
+    }
+    */
+}
+
+template<typename T>
+void DynamicDampingScheme<T>::compute() {
+
+    if (this->x_star.size() == this->x.size() && this->x.size() == x_prev.size()) {
+        // Central difference scheme for second order derivative
+        Eigen::VectorXd d2x_dt2 = this->x_star.array() - 2*this->x.array() + x_prev.array();
+        factor = (1.0/(1.0 + 10*(d2x_dt2.array().abs())));
+        x_prev = this->x;
+        this->x = (1 - factor.array()) * this->x.array() + factor.array() * this->x_star.array();
+    } else {
+        x_prev = this->x;
+        this->x = this->x_star;
+    }
+
+}
+
+}   // namespace mmft

src/hybridDynamics/LinearDecoupling.h

@@ -0,0 +1,49 @@
+/**
+ * @file Naive.h
+ */
+
+#pragma once
+
+#include <memory>
+#include <unordered_map>
+#include <vector>
+
+namespace mmft{
+
+/**
+ * @brief The Naive Scheme is an update scheme that sets the relaxation factor of an iterative update scheme to a 
+ * given constant for all nodes.
+ */
+template<typename T>
+class LinearDecouplingScheme : public Scheme<T> {
+
+private:
+
+    Eigen::VectorXd x_prev;
+    Eigen::VectorXd factor;
+    Eigen::VectorXd lambda;
+    Eigen::VectorXd omega;
+
+    T zeta = 10;
+    Eigen::VectorXd b = Eigen::VectorXd::Zero(2);
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(2,2);
+    Eigen::VectorXd ab = Eigen::VectorXd::Zero(2);
+
+public:
+    /**
+     * @brief Constructor of the Naive Scheme with provided constants.
+     * @param[in] modules The map of modules with boundary nodes upon which this scheme acts.
+     * @param[in] alpha The relaxation value for the pressure value update.
+     * @param[in] beta The relaxation value of the flow rate value update.
+     * @param[in] theta The amount of LBM stream and collide cycles between updates for a module.
+     */
+    LinearDecouplingScheme();
+
+    /**
+     * @brief Compute the update value alpha according to alpha = 1/(1 + 10*d2xdt2).
+     */
+    void compute() override;
+
+};
+
+}   // namespace mmft

src/hybridDynamics/LinearDecoupling.hh

@@ -0,0 +1,74 @@
+#include "LinearDecoupling.h"
+
+namespace mmft {
+
+template<typename T>
+LinearDecouplingScheme<T>::LinearDecouplingScheme() :
+    Scheme<T>() 
+{ 
+    /*
+    this->theta.try_emplace(0,theta_);
+    for (int i = 0; i < size_; i++) {
+        this->alpha.try_emplace(i, 0.0);
+    }
+    */
+}
+
+template<typename T>
+void LinearDecouplingScheme<T>::compute() {
+
+    if (this->x_star.size() == this->x.size() && this->x.size() == x_prev.size()) {
+        // Central difference scheme for second order derivative
+        Eigen::VectorXd d2x_dt2 = this->x_star.array() - 2*this->x.array() + x_prev.array();
+
+        lambda = d2x_dt2.array()/this->x.array();
+
+        int n = 0;
+        T l_sum = 0.0;
+        T w_sum = 0.0;
+        T inv_l_sum = 0.0;
+        T inv_w_sum = 0.0;
+
+        for (int i = 0; i < lambda.size(); i++) {
+            if (lambda[i] > 1e-12) {
+                T w = sqrt(lambda[i]);
+                l_sum += lambda[i];
+                w_sum += w;
+                inv_l_sum += 1.0/lambda[i];
+                inv_w_sum += 1.0/w;
+                n++;
+            }
+        }
+
+
+        A(0,0) = 0.25*inv_l_sum;
+        A(1,0) = 0.25*n;
+        A(0,1) = 0.25*n;
+        A(1,1) = 0.25*l_sum;
+
+        b(0) = 0.5*inv_w_sum*zeta;
+        b(1) = 0.5*w_sum*zeta;
+
+        ab = A.colPivHouseholderQr().solve(b);
+
+        T factor_a = (1.0)/(1.0 + ab[0]);
+        T factor_b = (ab[1] - 1.0)/(1.0 + ab[0]);
+        T factor_c = 1.0 + (-1.0)/(1.0 + ab[0]);
+
+        std::cout << "alpha: " << ab[0] << "\tbeta: " << ab[1] << std::endl;
+        std::cout << "a: " << factor_a << "\tb: " << factor_b << "\tc: " << factor_c << std::endl;
+
+        factor = (1.0/(1.0 + 10*(d2x_dt2.array().abs())));
+
+        x_prev = this->x;
+
+        this->x = factor_a*this->x_star.array() + factor_b*lambda.array()*this->x_star.array() + factor_c*this->x.array();
+        //this->x = (1 - factor.array()) * this->x.array() + factor.array() * this->x_star.array();
+    } else {
+        x_prev = this->x;
+        this->x = this->x_star;
+    }
+
+}
+
+}   // namespace mmft

src/hybridDynamics/Naive.h

@@ -71,6 +71,11 @@ class NaiveScheme : public Scheme<T> {
      * @param[in] theta The amount of LBM stream and collide cycles between updates for a module.
      */
     NaiveScheme(std::unordered_map<int, T> alpha, std::unordered_map<int, T> beta, std::unordered_map<int, int> theta);
+
+    /**
+     * @brief Calculate the update values according to the update scheme.
+     */
+    virtual void compute() override;
 
 };
 

src/hybridDynamics/Naive.hh

@@ -33,4 +33,14 @@ template<typename T>
 NaiveScheme<T>::NaiveScheme(std::unordered_map<int, T> alpha_, std::unordered_map<int, T> beta_, std::unordered_map<int, int> theta_) :
     Scheme<T>(alpha_, beta_, theta_) { }
 
+template<typename T>
+void NaiveScheme<T>::compute() { 
+    if (this->x.size() == this->x_star.size()) {
+        this->x = (1 - this->alpha.array()) * this->x.array() + this->alpha.array() * this->x_star.array();
+    } else {
+        // The first iteration the vector lengths might differ
+        this->x = this->x_star;
+    }
+}
+
 }   // namespace mmft