diff --git a/examples/Hybrid/Network1a.JSON b/examples/Hybrid/Network1a.JSON index e4efcd2..a051cb3 100644 --- a/examples/Hybrid/Network1a.JSON +++ b/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, diff --git a/examples/Hybrid/Network2a.JSON b/examples/Hybrid/Network2a.JSON index 356d5bf..e970841 100644 --- a/examples/Hybrid/Network2a.JSON +++ b/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, diff --git a/examples/Hybrid/Network3a.JSON b/examples/Hybrid/Network3a.JSON index 51f77a1..c5e3702 100644 --- a/examples/Hybrid/Network3a.JSON +++ b/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, diff --git a/examples/Hybrid/Network4a.JSON b/examples/Hybrid/Network4a.JSON index 752a649..a54e9c6 100644 --- a/examples/Hybrid/Network4a.JSON +++ b/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, diff --git a/src/baseSimulator.h b/src/baseSimulator.h index 1446e65..2308e38 100644 --- a/src/baseSimulator.h +++ b/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" diff --git a/src/baseSimulator.hh b/src/baseSimulator.hh index 0420699..78153cf 100644 --- a/src/baseSimulator.hh +++ b/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" diff --git a/src/hybridDynamics/CMakeLists.txt b/src/hybridDynamics/CMakeLists.txt index 69e655e..81194a8 100644 --- a/src/hybridDynamics/CMakeLists.txt +++ b/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 ) diff --git a/src/hybridDynamics/DynamicDamping.h b/src/hybridDynamics/DynamicDamping.h new file mode 100644 index 0000000..4145744 --- /dev/null +++ b/src/hybridDynamics/DynamicDamping.h @@ -0,0 +1,42 @@ +/** + * @file Naive.h + */ + +#pragma once + +#include +#include +#include + +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 +class DynamicDampingScheme : public Scheme { + +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 \ No newline at end of file diff --git a/src/hybridDynamics/DynamicDamping.hh b/src/hybridDynamics/DynamicDamping.hh new file mode 100644 index 0000000..2b9f778 --- /dev/null +++ b/src/hybridDynamics/DynamicDamping.hh @@ -0,0 +1,33 @@ +#include "DynamicDamping.h" + +namespace mmft { + +template +DynamicDampingScheme::DynamicDampingScheme(int theta_, int size_) : + Scheme() +{ + /* + this->theta.try_emplace(0,theta_); + for (int i = 0; i < size_; i++) { + this->alpha.try_emplace(i, 0.0); + } + */ +} + +template +void DynamicDampingScheme::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 \ No newline at end of file diff --git a/src/hybridDynamics/LinearDecoupling.h b/src/hybridDynamics/LinearDecoupling.h new file mode 100644 index 0000000..a379711 --- /dev/null +++ b/src/hybridDynamics/LinearDecoupling.h @@ -0,0 +1,49 @@ +/** + * @file Naive.h + */ + +#pragma once + +#include +#include +#include + +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 +class LinearDecouplingScheme : public Scheme { + +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 \ No newline at end of file diff --git a/src/hybridDynamics/LinearDecoupling.hh b/src/hybridDynamics/LinearDecoupling.hh new file mode 100644 index 0000000..b6a7029 --- /dev/null +++ b/src/hybridDynamics/LinearDecoupling.hh @@ -0,0 +1,74 @@ +#include "LinearDecoupling.h" + +namespace mmft { + +template +LinearDecouplingScheme::LinearDecouplingScheme() : + Scheme() +{ + /* + this->theta.try_emplace(0,theta_); + for (int i = 0; i < size_; i++) { + this->alpha.try_emplace(i, 0.0); + } + */ +} + +template +void LinearDecouplingScheme::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 \ No newline at end of file diff --git a/src/hybridDynamics/Naive.h b/src/hybridDynamics/Naive.h index be18f41..1566371 100644 --- a/src/hybridDynamics/Naive.h +++ b/src/hybridDynamics/Naive.h @@ -71,6 +71,11 @@ class NaiveScheme : public Scheme { * @param[in] theta The amount of LBM stream and collide cycles between updates for a module. */ NaiveScheme(std::unordered_map alpha, std::unordered_map beta, std::unordered_map theta); + + /** + * @brief Calculate the update values according to the update scheme. + */ + virtual void compute() override; }; diff --git a/src/hybridDynamics/Naive.hh b/src/hybridDynamics/Naive.hh index 8a5a897..6607cff 100644 --- a/src/hybridDynamics/Naive.hh +++ b/src/hybridDynamics/Naive.hh @@ -33,4 +33,14 @@ template NaiveScheme::NaiveScheme(std::unordered_map alpha_, std::unordered_map beta_, std::unordered_map theta_) : Scheme(alpha_, beta_, theta_) { } +template +void NaiveScheme::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 \ No newline at end of file diff --git a/src/hybridDynamics/Scheme.h b/src/hybridDynamics/Scheme.h index 77f4811..af039a4 100644 --- a/src/hybridDynamics/Scheme.h +++ b/src/hybridDynamics/Scheme.h @@ -8,6 +8,11 @@ #include #include +#include "Eigen/Dense" + +using Eigen::MatrixXd; +using Eigen::VectorXd; + namespace arch { // Forward declared dependencies @@ -28,9 +33,12 @@ class Scheme { protected: - std::unordered_map alpha; // relaxation value for pressure value updates on the Abstract-CFD interface nodes - std::unordered_map beta; // relaxation value for pressure value flow rate on the Abstract-CFD interface nodes - std::unordered_map theta; // Amount of LBM collide and stream iterations between update steps + Eigen::VectorXd alpha; // relaxation value for pressure value updates on the Abstract-CFD interface nodes + Eigen::VectorXd beta; // relaxation value for pressure value flow rate on the Abstract-CFD interface nodes + Eigen::VectorXd theta; // Amount of LBM collide and stream iterations between update steps + + Eigen::VectorXd x_star; + Eigen::VectorXd x; /** * @brief Default constructor of a Scheme. @@ -132,30 +140,46 @@ class Scheme { void setTheta(std::unordered_map theta); /** - * @brief Returns the relaxation value for pressure updates for node with nodeId. - * @param[in] nodeId The node for which alpha is returned. - * @returns alpha. + * @brief Sets the latest x* values. + * @param[in] x_star The vector x*. */ - T getAlpha(int nodeId) const; + void setXstar(const Eigen::VectorXd& x_star); /** - * @brief Returns the relaxation value for pressure updates for all nodes. - * @returns Map of alpha values. + * @brief Calculate the update values according to the update scheme. */ - const std::unordered_map& getAlpha() const; + virtual void compute() = 0; /** - * @brief Returns the relaxation value for flow rate updates for node with nodeId. - * @param[in] nodeId The node for which beta is returned. - * @returns beta. + * @brief Returns the latest x value, computed by the update scheme. */ - T getBeta(int nodeId) const; + const Eigen::VectorXd& getX(); - /** - * @brief Returns the relaxation value for flow rate updates for all nodes. - * @returns Map of beta values. - */ - const std::unordered_map& getBeta() const; + // /** + // * @brief Returns the relaxation value for pressure updates for node with nodeId. + // * @param[in] nodeId The node for which alpha is returned. + // * @returns alpha. + // */ + // T getAlpha(int nodeId) const; + + // /** + // * @brief Returns the relaxation value for pressure updates for all nodes. + // * @returns Map of alpha values. + // */ + // const std::unordered_map& getAlpha() const; + + // /** + // * @brief Returns the relaxation value for flow rate updates for node with nodeId. + // * @param[in] nodeId The node for which beta is returned. + // * @returns beta. + // */ + // T getBeta(int nodeId) const; + + // /** + // * @brief Returns the relaxation value for flow rate updates for all nodes. + // * @returns Map of beta values. + // */ + // const std::unordered_map& getBeta() const; /** * @brief Returns the number of LBM iterations between update steps for module with moduleId. @@ -164,11 +188,11 @@ class Scheme { */ int getTheta(int moduleId) const; - /** - * @brief Returns the number of LBM iterations between update steps for all modules. - * @returns Map of theta values. - */ - const std::unordered_map& getTheta() const; + // /** + // * @brief Returns the number of LBM iterations between update steps for all modules. + // * @returns Map of theta values. + // */ + // const std::unordered_map& getTheta() const; }; diff --git a/src/hybridDynamics/Scheme.hh b/src/hybridDynamics/Scheme.hh index 9f34417..04713c8 100644 --- a/src/hybridDynamics/Scheme.hh +++ b/src/hybridDynamics/Scheme.hh @@ -6,18 +6,19 @@ template Scheme::Scheme() { } template -Scheme::Scheme(std::unordered_map alpha_, std::unordered_map beta_, std::unordered_map theta_) : - alpha(alpha_), beta(beta_), theta(theta_) { } +Scheme::Scheme(std::unordered_map alpha_, std::unordered_map beta_, std::unordered_map theta_) + { } + //alpha(alpha_), beta(beta_), theta(theta_) { } template Scheme::Scheme(const std::unordered_map>>& modules_, T alpha_, T beta_, int theta_) { for (auto& [key, module] : modules_) { for (auto& [nodeId, node] : module->getNodes()) { - alpha.try_emplace(nodeId, alpha_); - beta.try_emplace(nodeId, beta_); + //alpha.try_emplace(nodeId, alpha_); + //beta.try_emplace(nodeId, beta_); } - theta.try_emplace(module->getId(), theta_); + //theta.try_emplace(module->getId(), theta_); } } @@ -25,98 +26,109 @@ template Scheme::Scheme(const std::shared_ptr> module_, T alpha_, T beta_, int theta_) { for (auto& [nodeId, node] : module_->getNodes()) { - alpha.try_emplace(nodeId, alpha_); - beta.try_emplace(nodeId, beta_); + //alpha.try_emplace(nodeId, alpha_); + //beta.try_emplace(nodeId, beta_); } - theta.try_emplace(module_->getId(), theta_); + //theta.try_emplace(module_->getId(), theta_); } template -Scheme::Scheme(const std::shared_ptr> module_, std::unordered_map alpha_, std::unordered_map beta_, int theta_) : - alpha(alpha_), beta(beta_) +Scheme::Scheme(const std::shared_ptr> module_, std::unordered_map alpha_, std::unordered_map beta_, int theta_) + //alpha(alpha_), beta(beta_) { - theta.try_emplace(module_->getId(), theta_); + //theta.try_emplace(module_->getId(), theta_); } template void Scheme::setAlpha(T alpha_) { - for (auto a : alpha) { + /*for (auto a : alpha) { a = alpha_; - } + }*/ } template void Scheme::setAlpha(int nodeId_, T alpha_) { - alpha.at(nodeId_) = alpha_; + //alpha.at(nodeId_) = alpha_; } template void Scheme::setAlpha(std::unordered_map alpha_) { - alpha = alpha_; + //alpha = alpha_; } template void Scheme::setBeta(T beta_) { - for (auto b : beta) { + /*for (auto b : beta) { b = beta_; - } + }*/ } template void Scheme::setBeta(int nodeId_, T beta_) { - beta.at(nodeId_) = beta_; + //beta.at(nodeId_) = beta_; } template void Scheme::setBeta(std::unordered_map beta_) { - beta = beta_; + //beta = beta_; } template void Scheme::setTheta(int theta_) { - for (auto t : theta) { + /*for (auto t : theta) { t = theta_; - } + }*/ } template void Scheme::setTheta(int moduleId_, int theta_) { - theta.at(moduleId_) = theta_; + //theta.at(moduleId_) = theta_; } template void Scheme::setTheta(std::unordered_map theta_) { - theta = theta_; + //theta = theta_; } template -T Scheme::getAlpha(int nodeId_) const { - return alpha.at(nodeId_); +void Scheme::setXstar(const Eigen::VectorXd& x_star_) { + x_star = x_star_; } template -const std::unordered_map& Scheme::getAlpha() const { - return alpha; +const Eigen::VectorXd& Scheme::getX() { + return x; } -template -T Scheme::getBeta(int nodeId_) const { - return beta.at(nodeId_); -} +// template +// T Scheme::getAlpha(int nodeId_) const { +// return alpha.at(nodeId_); +// } -template -const std::unordered_map& Scheme::getBeta() const { - return beta; -} +// template +// const std::unordered_map& Scheme::getAlpha() const { +// return alpha; +// } + +// template +// T Scheme::getBeta(int nodeId_) const { +// return beta.at(nodeId_); +// } + +// template +// const std::unordered_map& Scheme::getBeta() const { +// return beta; +// } template int Scheme::getTheta(int moduleId_) const { - return theta.at(moduleId_); + //return theta.at(moduleId_); + return 1; } -template -const std::unordered_map& Scheme::getTheta() const { - return theta; -} +// template +// const std::unordered_map& Scheme::getTheta() const { +// return theta; +// } } // namespace mmft \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index ca4b9f1..0dce3b0 100755 --- a/src/main.cpp +++ b/src/main.cpp @@ -6,6 +6,8 @@ #include #include +#include + #ifdef USE_ESSLBM #include #endif @@ -18,7 +20,12 @@ int main(int argc, char const* argv []) { MPI_Init(NULL,NULL); #endif - std::string file = argv[1]; + //std::string file = argv[1]; + std::string n = argv[1]; + std::string m = argv[2]; + std::string i = argv[3]; + + std::string file = "../../network_generation/mVLSI-Benchmark/n_"+n+"_m_"+m+"/Network_n"+n+"_m"+m+"_"+i+".json"; // Load and set the network from a JSON file std::cout << "[Main] Create network object..." << std::endl; @@ -30,15 +37,38 @@ int main(int argc, char const* argv []) { std::cout << "[Main] Simulation..." << std::endl; // Perform simulation and store results + auto start_time = std::chrono::high_resolution_clock::now(); testSimulation.simulate(); + auto end_time = std::chrono::high_resolution_clock::now(); + auto time = end_time - start_time; std::cout << "[Main] Results..." << std::endl; + std::cout << time/std::chrono::milliseconds(1) << "," << std::to_string(testSimulation.getNodalAnalysis()->iteration) << std::endl; + // Print the results - testSimulation.getSimulationResults()->printStates(); + //testSimulation.getSimulationResults()->printStates(); //std::cout << "Write diffusive mixtures" << std::endl; //testSimulation.getSimulationResults()->writeMixture(1); + // Write csv + try { + std::string outputFileName = "Time_n"+n+"_m"+m+"_i"+i+".csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + outputFile << time/std::chrono::milliseconds(1) << "," << std::to_string(testSimulation.getNodalAnalysis()->iteration) << std::endl; + + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + #ifdef USE_ESSLBM MPI_Finalize(); #endif diff --git a/src/nodalAnalysis/NodalAnalysis.h b/src/nodalAnalysis/NodalAnalysis.h index 0d736ab..7b0c0e3 100644 --- a/src/nodalAnalysis/NodalAnalysis.h +++ b/src/nodalAnalysis/NodalAnalysis.h @@ -22,6 +22,14 @@ class Network; } // namespace arch +namespace mmft { + +// Forward declared dependencies +template +class Scheme; + +} // namespace scheme + namespace sim { // Forward declared dependencies @@ -43,8 +51,11 @@ class NodalAnalysis { bool pressureConvergence; Eigen::MatrixXd A; // matrix A = [G, B; C, D] + Eigen::MatrixXd Aclean; // matrix A = [G, B; C, D] Eigen::VectorXd z; // vector z = [i; e] Eigen::VectorXd x; // vector x = [v; j] + Eigen::VectorXd x_prev; // vector x^(n-1) + Eigen::VectorXd lambda; // vector of (eigen frequency)^2 const std::vector baseline1 = {0.0, 1000.0, 1000.0, 1000.0, 859.216, 791.962, 859.216, 753.628, 753.628, 422.270, 0.0}; @@ -63,17 +74,25 @@ class NodalAnalysis { 522.332, 454.959, 355.563, 217.833, 0.0}; std::vector> errors; + std::vector> xn; // x_n + std::vector> xn_dt; // x_dot_n (1st "time" derivative) + std::vector> xn_dt2; // x_double dot_n (2nd "time" derivative) + std::vector> lambdan; // x_double dot_n (2nd "time" derivative) + std::vector> Alphan; // x_double dot_n (2nd "time" derivative) + + std::vector Alpha; std::unordered_set conductingNodeIds; std::unordered_map groundNodeIds; - void readConductance(); // loop through channels and build matrix G - void updateReferenceP(); // update the reference pressure for each group - void readPressurePumps(); // loop through pressure pumps and build matrix B, C and vector e - void readFlowRatePumps(); // loop through flowRate pumps and build vector i - void solve(); // solve equation x = A^(-1) * z - void setResults(); // set pressure of nodes to v and flow rate at pressure pumps to j - void initGroundNodes(); // initialize the ground nodes of the groups + void readConductance(); // loop through channels and build matrix G + void updateReferenceP(); // update the reference pressure for each group + void readPressurePumps(); // loop through pressure pumps and build matrix B, C and vector e + void readFlowRatePumps(); // loop through flowRate pumps and build vector i + void solve(); // solve equation x = A^(-1) * z + void solve(std::shared_ptr>& scheme); // solve equation x = updateScheme(A^(-1) * z), for hybrid simulation + void setResults(); // set pressure of nodes to v and flow rate at pressure pumps to j + void initGroundNodes(); // initialize the ground nodes of the groups void clear(); // For hybrid simulations @@ -109,8 +128,9 @@ class NodalAnalysis { * * @param[in] network Pointer to the network on which the nodal analysis is conducted. Results are stored in the network's nodes. * @param[in] cfdSimulators The cfd simulators for a hybrid simulation + * @param[in] scheme The iterative update scheme used for a hybrid simulation */ - bool conductNodalAnalysis(std::unordered_map>>& cfdSimulators); + bool conductNodalAnalysis(std::unordered_map>>& cfdSimulators, std::shared_ptr>& scheme); void writeSystem(); diff --git a/src/nodalAnalysis/NodalAnalysis.hh b/src/nodalAnalysis/NodalAnalysis.hh index 81caecb..4ce1dae 100644 --- a/src/nodalAnalysis/NodalAnalysis.hh +++ b/src/nodalAnalysis/NodalAnalysis.hh @@ -38,8 +38,10 @@ NodalAnalysis::NodalAnalysis(const arch::Network* network_) { int nNodesAndPressurePumps = nNodes + nPressurePumps + network->getModules().size(); A = Eigen::MatrixXd::Zero(nNodesAndPressurePumps, nNodesAndPressurePumps); + Aclean = Eigen::MatrixXd::Zero(12, 12); z = Eigen::VectorXd::Zero(nNodesAndPressurePumps); x = Eigen::VectorXd::Zero(nNodesAndPressurePumps); + x_prev = Eigen::VectorXd::Zero(nNodesAndPressurePumps); } @@ -70,10 +72,12 @@ void NodalAnalysis::clear() { int nNodesAndPressurePumps = nNodes + nPressurePumps + groundNodeIds.size(); A = Eigen::MatrixXd::Zero(nNodesAndPressurePumps, nNodesAndPressurePumps); + Aclean = Eigen::MatrixXd::Zero(12, 12); z = Eigen::VectorXd::Zero(nNodesAndPressurePumps); - x = Eigen::VectorXd::Zero(nNodesAndPressurePumps); + //x = Eigen::VectorXd::Zero(nNodesAndPressurePumps); A.setZero(); + Aclean.setZero(); z.setZero(); } @@ -89,14 +93,14 @@ void NodalAnalysis::conductNodalAnalysis() { } template -bool NodalAnalysis::conductNodalAnalysis(std::unordered_map>>& cfdSimulators) { +bool NodalAnalysis::conductNodalAnalysis(std::unordered_map>>& cfdSimulators, std::shared_ptr>& scheme) { clear(); readConductance(); readCfdSimulators(cfdSimulators); updateReferenceP(); readPressurePumps(); readFlowRatePumps(); - solve(); + solve(scheme); setResults(); writeCfdSimulators(cfdSimulators); initGroundNodes(cfdSimulators); @@ -190,14 +194,36 @@ template void NodalAnalysis::solve() { // solve equation x = A^(-1) * z x = A.colPivHouseholderQr().solve(z); + + iteration++; +} - if (iteration % 1 == 0) { - std::vector tmpError; - for (long unsigned int i = 0; i < baseline3.size(); i++) { - tmpError.push_back(baseline3[i] - x[i]); +template +void NodalAnalysis::solve(std::shared_ptr>& scheme) { + // solve equation x_star = A^(-1) * z + Eigen::VectorXd x_star = A.colPivHouseholderQr().solve(z); + + if (x_star.size() == x.size() && x.size() == x_prev.size()) { + std::vector tmp_xn; + std::vector tmp_xn_dt; + std::vector tmp_xn_dt2; + if (iteration % 1 == 0) { + for (int i = 0; i < x.size(); i++) { + tmp_xn.push_back(x[i]); + tmp_xn_dt.push_back(x_star[i] - x[i]); + tmp_xn_dt2.push_back(x_star[i] - 2*x[i] + x_prev[i]); + } + xn.push_back(tmp_xn); + xn_dt.push_back(tmp_xn_dt); + xn_dt2.push_back(tmp_xn_dt2); } - errors.push_back(tmpError); } + + x_prev = x; + + scheme->setXstar(x_star); + scheme->compute(); + x = scheme->getX(); iteration++; } @@ -379,15 +405,8 @@ void NodalAnalysis::writeCfdSimulators(std::unordered_mapgetPressure(); - T set_pressure = 0.0; - if (old_pressure > 0 ) { - set_pressure = old_pressure + cfdSimulator.second->getAlpha(key) * ( new_pressure - old_pressure ); - } else { - set_pressure = new_pressure; - } - pressures_.at(key) = set_pressure; - - if (abs(old_pressure - new_pressure) > 1e-2) { + pressures_.at(key) = new_pressure; + if (abs(old_pressure - new_pressure) > 1e-1) { pressureConvergence = false; } } @@ -395,15 +414,8 @@ void NodalAnalysis::writeCfdSimulators(std::unordered_mapgetOpenings().at(key).width; - T set_flowRate = 0.0; - if (old_flowRate > 0 ) { - set_flowRate = old_flowRate + cfdSimulator.second->getBeta(key) * ( new_flowRate - old_flowRate ); - } else { - set_flowRate = new_flowRate; - } - flowRates_.at(key) = set_flowRate; - - if (abs(old_flowRate - new_flowRate) > 1e-2) { + flowRates_.at(key) = new_flowRate; + if (abs(old_flowRate - new_flowRate) > 1e-1) { pressureConvergence = false; } } @@ -445,10 +457,12 @@ void NodalAnalysis::printSystem() { template void NodalAnalysis::writeSystem() { + std::vector indices = {1,2,3,4,5,6,7,8,9,11,12}; //std::vector indices = {3,5,6,7,8}; // #2 - std::vector indices = {3, 5, 8, 9, 10, 11, 12, 13, 14}; // #3 + //std::vector indices = {3, 5, 8, 9, 10, 11, 12, 13, 14}; // #3 //std::vector indices = {2,8,9,4,10,11,5,12,13,7,14}; // #4 + // Write errors try { std::string outputFileName = "Errors.csv"; std::cout << "Generating CSV file: " << outputFileName << std::endl; @@ -463,7 +477,142 @@ void NodalAnalysis::writeSystem() { for (const auto& index : indices) { outputFile << vec[index] << ","; } - outputFile << vec[15]; + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write xn + try { + std::string outputFileName = "x.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : xn) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write dxn/dt + try { + std::string outputFileName = "dxn_dt.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : xn_dt) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write d2xn/dt2 + try { + std::string outputFileName = "d2xn_dt2.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : xn_dt2) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write lambda + try { + std::string outputFileName = "lambda.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : lambdan) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write alpha + try { + std::string outputFileName = "alpha.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : Alphan) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; outputFile << "\n"; it++; } diff --git a/src/nodalAnalysis/NodalAnalysis_backup.h b/src/nodalAnalysis/NodalAnalysis_backup.h new file mode 100644 index 0000000..e2996ae --- /dev/null +++ b/src/nodalAnalysis/NodalAnalysis_backup.h @@ -0,0 +1,131 @@ +/** + * @file NodalAnalysis.h + */ + +#pragma once + +#include +#include +#include +#include + +#include "Eigen/Dense" + +using Eigen::MatrixXd; +using Eigen::VectorXd; + +namespace arch { + +// Forward declared dependencies +template +class Network; + +} // namespace arch + +namespace sim { + +// Forward declared dependencies +template +class CFDSimulator; + +} // namespace sim + +namespace nodal { + +template +class NodalAnalysis { +private: + const arch::Network* network; + + int nNodes; // Number of nodes + int nPressurePumps; // Number of pressurePumps + + bool pressureConvergence; + + Eigen::MatrixXd A; // matrix A = [G, B; C, D] + Eigen::MatrixXd Aclean; // matrix A = [G, B; C, D] + Eigen::VectorXd z; // vector z = [i; e] + Eigen::VectorXd x; // vector x = [v; j] + Eigen::VectorXd x_prev; // vector x^(n-1) + Eigen::VectorXd lambda; // vector of (eigen frequency)^2 + + const std::vector baseline1 = {0.0, 1000.0, 1000.0, 1000.0, 859.216, + 791.962, 859.216, 753.628, 753.628, 422.270, 0.0}; + + const std::vector baseline2 = {0.0, 1000.0, 577.014, 729.934, 491.525, + 469.468, 641.131, 607.210, 427.408, 270.394, 0.0}; + + const std::vector baseline3 = {0.0, 1000.0, 687.204, 801.008, 625.223, + 601.422, 428.435, 428.007, 331.775, 733.690, + 703.349, 578.736, 525.094, 529.993, 326.021, + 198.163, 0.0}; + + const std::vector baseline4 = {0.0, 1000.0, 781.127, 616.813, + 604.130, 535.202, 447.967, 371.026, + 698.425, 695.905, 552.151, 540.371, + 522.332, 454.959, 355.563, 217.833, 0.0}; + + std::vector> errors; + std::vector> xn; // x_n + std::vector> xn_dt; // x_dot_n (1st "time" derivative) + std::vector> xn_dt2; // x_double dot_n (2nd "time" derivative) + std::vector> lambdan; // x_double dot_n (2nd "time" derivative) + std::vector> Alphan; // x_double dot_n (2nd "time" derivative) + + std::vector Alpha; + + std::unordered_set conductingNodeIds; + std::unordered_map groundNodeIds; + + void readConductance(); // loop through channels and build matrix G + void updateReferenceP(); // update the reference pressure for each group + void readPressurePumps(); // loop through pressure pumps and build matrix B, C and vector e + void readFlowRatePumps(); // loop through flowRate pumps and build vector i + void solve(); // solve equation x = A^(-1) * z + void setResults(); // set pressure of nodes to v and flow rate at pressure pumps to j + void initGroundNodes(); // initialize the ground nodes of the groups + void clear(); + + // For hybrid simulations + void readCfdSimulators(std::unordered_map>>& cfdSimulators); + void writeCfdSimulators(std::unordered_map>>& cfdSimulators); + void initGroundNodes(std::unordered_map>>& cfdSimulators); + + // Helper functions + bool contains( const std::unordered_set& set, int key); + bool contains( const std::unordered_map& map, int key); + void printSystem(); + +public: + + int iteration = 0; + + /** + * @brief Creates a NodalAnalysis object + */ + NodalAnalysis(const arch::Network* network); + + /** + * @brief Conducts the Modifed Nodal Analysis (e.g., http://qucs.sourceforge.net/tech/node14.html) and computes the pressure levels for each node. + * Hence, the passed nodes contain the final pressure levels when the function is finished. + * + * @param[in] network Pointer to the network on which the nodal analysis is conducted. Results are stored in the network's nodes. + */ + void conductNodalAnalysis(); + + /** + * @brief Conducts the Modifed Nodal Analysis (e.g., http://qucs.sourceforge.net/tech/node14.html) and computes the pressure levels for each node. + * Hence, the passed nodes contain the final pressure levels when the function is finished. + * + * @param[in] network Pointer to the network on which the nodal analysis is conducted. Results are stored in the network's nodes. + * @param[in] cfdSimulators The cfd simulators for a hybrid simulation + */ + bool conductNodalAnalysis(std::unordered_map>>& cfdSimulators); + + + void writeSystem(); + +}; + + +} // namespace nodal \ No newline at end of file diff --git a/src/nodalAnalysis/NodalAnalysis_backup.hh b/src/nodalAnalysis/NodalAnalysis_backup.hh new file mode 100644 index 0000000..bd38cfe --- /dev/null +++ b/src/nodalAnalysis/NodalAnalysis_backup.hh @@ -0,0 +1,688 @@ +#include "NodalAnalysis.h" + +namespace nodal { + +template +NodalAnalysis::NodalAnalysis(const arch::Network* network_) { + network = network_; + nNodes = network->getNodes().size() + network->getVirtualNodes(); + + // loop through modules + for (const auto& [key, module] : network->getModules()) { + // For now only LBM modules implemented. + #ifndef USE_ESSLBM + assert(module->getModuleType() == arch::ModuleType::LBM); + #elif USE_ESSLBM + assert(module->getModuleType() == arch::ModuleType::ESS_LBM); + #endif + } + + // Sort nodes into conducting nodes and ground nodes. + // First loop, all nodes with id > 0 are conducting nodes. + int iPump = nNodes; + for (const auto& [key, group] : network->getGroups()) { + for (const auto& nodeId : group->nodeIds) { + // The node is a conducting node + if(!network->getNodes().at(nodeId)->getGround() && nodeId != group->groundNodeId) { + conductingNodeIds.emplace(nodeId); + } + // The node is an overall ground node, or counts as ground to a group + else if (!network->getNodes().at(nodeId)->getGround() && nodeId == group->groundNodeId) { + groundNodeIds.emplace(nodeId, iPump); + iPump++; + } + } + } + + nPressurePumps = network->getPressurePumps().size() + groundNodeIds.size(); + int nNodesAndPressurePumps = nNodes + nPressurePumps + network->getModules().size(); + + A = Eigen::MatrixXd::Zero(nNodesAndPressurePumps, nNodesAndPressurePumps); + Aclean = Eigen::MatrixXd::Zero(12, 12); + z = Eigen::VectorXd::Zero(nNodesAndPressurePumps); + x = Eigen::VectorXd::Zero(nNodesAndPressurePumps); + x_prev = Eigen::VectorXd::Zero(nNodesAndPressurePumps); + +} + +template +void NodalAnalysis::clear() { + + pressureConvergence = true; + + conductingNodeIds.clear(); + groundNodeIds.clear(); + + int iPump = network->getNodes().size() + network->getVirtualNodes() + network->getPressurePumps().size(); + + for (const auto& [key, group] : network->getGroups()) { + for (const auto& nodeId : group->nodeIds) { + // The node is a conducting node + if(!network->getNodes().at(nodeId)->getGround() && nodeId != group->groundNodeId) { + conductingNodeIds.emplace(nodeId); + } + // The node is an overall ground node, or counts as ground to a group + else if (!network->getNodes().at(nodeId)->getGround() && nodeId == group->groundNodeId) { + groundNodeIds.emplace(nodeId, iPump); + iPump++; + } + } + } + + int nNodesAndPressurePumps = nNodes + nPressurePumps + groundNodeIds.size(); + + A = Eigen::MatrixXd::Zero(nNodesAndPressurePumps, nNodesAndPressurePumps); + Aclean = Eigen::MatrixXd::Zero(12, 12); + z = Eigen::VectorXd::Zero(nNodesAndPressurePumps); + //x = Eigen::VectorXd::Zero(nNodesAndPressurePumps); + + A.setZero(); + Aclean.setZero(); + z.setZero(); +} + +template +void NodalAnalysis::conductNodalAnalysis() { + clear(); + readConductance(); + readPressurePumps(); + readFlowRatePumps(); + solve(); + setResults(); + initGroundNodes(); +} + +template +bool NodalAnalysis::conductNodalAnalysis(std::unordered_map>>& cfdSimulators) { + clear(); + readConductance(); + readCfdSimulators(cfdSimulators); + updateReferenceP(); + readPressurePumps(); + readFlowRatePumps(); + solve(); + setResults(); + writeCfdSimulators(cfdSimulators); + initGroundNodes(cfdSimulators); + return pressureConvergence; +} + +template +void NodalAnalysis::readConductance() { + // loop through channels and build matrix G + for (const auto& channel : network->getChannels()) { + auto nodeAMatrixId = channel.second->getNodeA(); + auto nodeBMatrixId = channel.second->getNodeB(); + const T conductance = 1. / channel.second->getResistance(); + + // main diagonal elements of G + if (!network->getNodes().at(nodeAMatrixId)->getGround()) { + A(nodeAMatrixId, nodeAMatrixId) += conductance; + } + + if (!network->getNodes().at(nodeBMatrixId)->getGround()) { + A(nodeBMatrixId, nodeBMatrixId) += conductance; + } + + // minor diagonal elements of G (if no ground node was present) + if (!network->getNodes().at(nodeAMatrixId)->getGround() && !network->getNodes().at(nodeBMatrixId)->getGround()) { + A(nodeAMatrixId, nodeBMatrixId) -= conductance; + A(nodeBMatrixId, nodeAMatrixId) -= conductance; + } + } +} + +template +void NodalAnalysis::updateReferenceP() { + // Update the reference pressure for each group + for (const auto& [key, group] : network->getGroups()) { + if (group->initialized) { + auto& node = network->getNodes().at(group->groundNodeId); + group->pRef = node->getPressure(); + int pumpId = groundNodeIds.at(group->groundNodeId); + + A(group->groundNodeId, pumpId) = 1; // matrix B + A(pumpId, group->groundNodeId) = 1; // matrix C + + z(pumpId) = node->getPressure(); + } + } +} + +template +void NodalAnalysis::readPressurePumps() { + int iPump = nNodes; + // loop through pressurePumps and build matrix B, C and vector e + for (const auto& pressurePump : network->getPressurePumps()) { + auto nodeAMatrixId = pressurePump.second->getNodeA(); + auto nodeBMatrixId = pressurePump.second->getNodeB(); + + if (contains(conductingNodeIds, nodeAMatrixId)) { + A(nodeAMatrixId, iPump) = -1; // matrix B + A(iPump, nodeAMatrixId) = -1; // matrix C + } + + if (contains(conductingNodeIds, nodeBMatrixId)) { + A(nodeBMatrixId, iPump) = 1; // matrix B + A(iPump, nodeBMatrixId) = 1; // matrix C + } + + z(iPump) = pressurePump.second->getPressure(); + + iPump++; + } +} + +template +void NodalAnalysis::readFlowRatePumps() { + // loop through flowRatePumps and build vector i + for (const auto& flowRatePump : network->getFlowRatePumps()) { + auto nodeAMatrixId = flowRatePump.second->getNodeA(); + auto nodeBMatrixId = flowRatePump.second->getNodeB(); + const T flowRate = flowRatePump.second->getFlowRate(); + + if (contains(conductingNodeIds, nodeAMatrixId)){ + z(nodeAMatrixId) = -flowRate; + } + if (contains(conductingNodeIds, nodeBMatrixId)){ + z(nodeBMatrixId) = flowRate; + } + } +} + +template +void NodalAnalysis::solve() { + // solve equation x = A^(-1) * z + Eigen::VectorXd x_next = A.colPivHouseholderQr().solve(z); + + std::vector tmp_Alpha (x_next.size(), 0.0); + + if (x_next.size() == x.size() && x.size() == x_prev.size()) { + lambda = (x_next.array() - 2*x.array() + x_prev.array())/x.array(); + std::vector tmp_xn; + std::vector tmp_xn_dt; + std::vector tmp_xn_dt2; + std::vector tmp_lambdan; + if (iteration % 1 == 0) { + for (int i = 0; i < x.size(); i++) { + tmp_xn.push_back(x[i]); + tmp_xn_dt.push_back(x_next[i] - x[i]); + tmp_xn_dt2.push_back(x_next[i] - 2*x[i] + x_prev[i]); + tmp_lambdan.push_back(lambda[i]); + tmp_Alpha[i] = (1.0/(1.0 + std::abs(x_next[i] - 2*x[i] + x_prev[i]))); + } + xn.push_back(tmp_xn); + xn_dt.push_back(tmp_xn_dt); + xn_dt2.push_back(tmp_xn_dt2); + lambdan.push_back(tmp_lambdan); + Alphan.push_back(tmp_Alpha); + } + } + + Alpha = tmp_Alpha; + + if (x_next.size() == x.size() && x.size() == x_prev.size()) { + Eigen::VectorXd eigenAlpha = (1.0/(1.0 + 10*((x_next.array() - 2*x.array() + x_prev.array()).abs()))); + //std::cout << "2nd derivative" << std::endl; + //std::cout << (x_next.array() - 2*x.array() + x_prev.array()).abs() << std::endl; + x_prev = x; + x = (1 - eigenAlpha.array()) * x.array() + (eigenAlpha.array()) * x_next.array(); + } else { + x_prev = x; + x = x_next; + } + + + +/* + if (x_next.size() == x.size() && x.size == x_prev.size()) { + for (long unsigned int i = 0; i < x.size(); i++) { + lamda = ()/ + } + } +*/ +/* + std::vector matrixIndices = {1,2,3,4,5,6,7,8,9,11,12,13}; + + std::cout << "\nA:\n" << A << std::endl; + if (A.size() > 145) { + for (int i=0; i<12; i++) { + for (int j=0; j<12; j++) { + Aclean(i,j) = A(matrixIndices[i], matrixIndices[j]); + } + } + } + + std::cout << "\nAclean:\n" << Aclean << std::endl; + std::cout << "\ndet(A):\n" << Aclean.determinant() << std::endl; + std::cout << "\nA-1:\n" << Aclean.inverse() << std::endl; + std::cout << "\nx:\n" << x << std::endl; + std::cout << "\nz:\n" << z << "\n" << std::endl; +*/ + if (iteration % 1 == 0) { + std::vector tmpError; + for (long unsigned int i = 0; i < baseline2.size(); i++) { + tmpError.push_back(baseline2[i] - x[i]); + } + errors.push_back(tmpError); + } + + iteration++; +} + +template +void NodalAnalysis::setResults() { + // set pressure of nodes to result value + for (const auto& [key, group] : network->getGroups()) { + for (auto nodeMatrixId : group->nodeIds) { + auto& node = network->getNodes().at(nodeMatrixId); + if (contains(conductingNodeIds, nodeMatrixId)) { + node->setPressure(x(nodeMatrixId)); + } else if (node->getGround()) { + node->setPressure(0.0); + } + } + } + + for (auto& channel : network->getChannels()) { + auto& nodeA = network->getNodes().at(channel.second->getNodeA()); + auto& nodeB = network->getNodes().at(channel.second->getNodeB()); + channel.second->setPressure(nodeA->getPressure() - nodeB->getPressure()); + } + + // set flow rate at pressure pumps + int iPump = nNodes; + for (auto& pressurePump : network->getPressurePumps()){ + pressurePump.second->setFlowRate(x(iPump)); + iPump++; + } +} + +template +void NodalAnalysis::initGroundNodes() { + // Initialize the ground nodes of the groups + for (const auto& [key, group] : network->getGroups()) { + if (!group->initialized && !group->grounded) { + T pMin = -1.0; + for (auto nodeId : group->nodeIds) { + if (pMin < 0.0) { + pMin = x(nodeId); + group->groundNodeId = nodeId; + } + if (x(nodeId) < pMin) { + pMin = x(nodeId); + group->groundNodeId = nodeId; + } + } + for (auto channelId : group->channelIds) { + if (network->getChannels().at(channelId)->getNodeA() == group->groundNodeId) { + group->groundChannelId = channelId; + } + else if (network->getChannels().at(channelId)->getNodeB() == group->groundNodeId) { + group->groundChannelId = channelId; + } + } + groundNodeIds.emplace(group->groundNodeId, 0); + conductingNodeIds.erase(group->groundNodeId); + group->initialized = true; + } + } +} + +template +void NodalAnalysis::initGroundNodes(std::unordered_map>>& cfdSimulators) { + // Initialize the ground nodes of the groups + for (const auto& [key, group] : network->getGroups()) { + if (!group->initialized && !group->grounded) { + T pMin = -1.0; + for (auto nodeId : group->nodeIds) { + if (pMin < 0.0) { + pMin = x(nodeId); + group->groundNodeId = nodeId; + } + if (x(nodeId) < pMin) { + pMin = x(nodeId); + group->groundNodeId = nodeId; + } + } + for (auto channelId : group->channelIds) { + if (network->getChannels().at(channelId)->getNodeA() == group->groundNodeId) { + group->groundChannelId = channelId; + } + else if (network->getChannels().at(channelId)->getNodeB() == group->groundNodeId) { + group->groundChannelId = channelId; + } + } + groundNodeIds.emplace(group->groundNodeId, 0); + conductingNodeIds.erase(group->groundNodeId); + group->initialized = true; + } + } + + // Loop over modules to set the ground nodes + for (const auto& [key, cfdSimulator] : cfdSimulators) { + if ( ! cfdSimulator->getInitialized() ) { + std::unordered_map flowRates_ = cfdSimulator->getFlowRates(); + std::unordered_map groundNodes; + for (const auto& [nodeId, node] : cfdSimulator->getModule()->getNodes()) { + T flowRate = 0.0; + if (contains(groundNodeIds, nodeId)) { + groundNodes.try_emplace(nodeId, true); + for (auto& [key, group] : network->getGroups()) { + if (nodeId == group->groundNodeId) { + T height = network->getChannels().at(group->groundChannelId)->getHeight(); + flowRate = network->getChannels().at(group->groundChannelId)->getFlowRate()/height; + } + } + } else { + groundNodes.try_emplace(nodeId, false); + } + flowRates_.at(nodeId) = flowRate; + } + cfdSimulator->setGroundNodes(groundNodes); + cfdSimulator->storeFlowRates(flowRates_); + cfdSimulator->setInitialized(true); + } + } +} + +template +void NodalAnalysis::readCfdSimulators(std::unordered_map>>& cfdSimulators) { + for (const auto& [key, cfdSimulator] : cfdSimulators) { + // If module is not initialized (1st loop), loop over channels of fully connected graph + if ( ! cfdSimulator->getInitialized() ) { + for (const auto& [key, channel] : cfdSimulator->getNetwork()->getChannels()) { + auto nodeAMatrixId = channel->getNodeA(); + auto nodeBMatrixId = channel->getNodeB(); + const T conductance = 1. / channel->getResistance(); + + // main diagonal elements of G + if (contains(conductingNodeIds, nodeAMatrixId)) { + A(nodeAMatrixId, nodeAMatrixId) += conductance; + } + + if (contains(conductingNodeIds, nodeBMatrixId)) { + A(nodeBMatrixId, nodeBMatrixId) += conductance; + } + + // minor diagonal elements of G (if no ground node was present) + if (contains(conductingNodeIds, nodeAMatrixId) && contains(conductingNodeIds, nodeBMatrixId)) { + A(nodeAMatrixId, nodeBMatrixId) -= conductance; + A(nodeBMatrixId, nodeAMatrixId) -= conductance; + } + } + } + /* If module is initialized, loop over boundary nodes and: + * - read pressure for ground nodes + * - read flow rate for conducting nodes + */ + else if ( cfdSimulator->getInitialized() ) { + for (const auto& [key, node] : cfdSimulator->getModule()->getNodes()) { + // Write the module's flowrates into vector i if the node is not a group's ground node + if (contains(conductingNodeIds, key)) { + T flowRate = cfdSimulator->getFlowRates().at(key) * cfdSimulator->getOpenings().at(key).height; + z(key) = -flowRate; + } + // Write module's pressure into matrix B, C and vector e + else if (contains(groundNodeIds, key)) { + T pressure = cfdSimulator->getPressures().at(key); + node->setPressure(pressure); + } + } + } + } +} + +template +void NodalAnalysis::writeCfdSimulators(std::unordered_map>>& cfdSimulators) { + + // Set the pressures and flow rates on the boundary nodes of the modules + for (auto& cfdSimulator : cfdSimulators) { + std::unordered_map old_pressures = cfdSimulator.second->getPressures(); + std::unordered_map old_flowrates = cfdSimulator.second->getFlowRates(); + std::unordered_map pressures_ = cfdSimulator.second->getPressures(); + std::unordered_map flowRates_ = cfdSimulator.second->getFlowRates(); + for (auto& [key, node] : cfdSimulator.second->getModule()->getNodes()){ + // Communicate pressure to the module + if (contains(conductingNodeIds, key)) { + T old_pressure = old_pressures.at(key); + T new_pressure = node->getPressure(); + T set_pressure = 0.0; + if (old_pressure > 0 ) { + //set_pressure = (1.0 - Alpha[key]) * old_pressure + Alpha[key] * new_pressure; + set_pressure = new_pressure; + } else { + set_pressure = new_pressure; + } + pressures_.at(key) = set_pressure; + + if (abs(old_pressure - new_pressure) > 1e-2) { + pressureConvergence = false; + } + } + // Communicate the flow rate to the module + else if (contains(groundNodeIds, key)) { + T old_flowRate = old_flowrates.at(key) ; + T new_flowRate = x(groundNodeIds.at(key)) / cfdSimulator.second->getOpenings().at(key).width; + T set_flowRate = 0.0; + if (old_flowRate > 0 ) { + //set_flowRate = (1.0 - Alpha[key]) * old_flowRate + Alpha[key] * new_flowRate; + set_flowRate = new_flowRate; + } else { + set_flowRate = new_flowRate; + } + flowRates_.at(key) = set_flowRate; + + if (abs(old_flowRate - new_flowRate) > 1e-2) { + pressureConvergence = false; + } + } + } + cfdSimulator.second->storePressures(pressures_); + cfdSimulator.second->storeFlowRates(flowRates_); + } +} + +template +bool NodalAnalysis::contains( const std::unordered_set& set, int key) { + bool contain = false; + for (auto& nodeId : set) { + if (key == nodeId) { + contain = true; + } + } + return contain; +} + +template +bool NodalAnalysis::contains( const std::unordered_map& map, int key) { + bool contain = false; + for (auto& [nodeId, pumpId] : map) { + if (key == nodeId) { + contain = true; + } + } + return contain; +} + +template +void NodalAnalysis::printSystem() { + std::cout << "Matrix A:\n" << A << "\n\n" << std::endl; + std::cout << "Vector z:\n" << z << "\n\n" << std::endl; + std::cout << "Vector x:\n" << x << "\n\n" << std::endl; +} + +template +void NodalAnalysis::writeSystem() { + + std::vector indices = {1,2,3,4,5,6,7,8,9,11,12}; + //std::vector indices = {3,5,6,7,8}; // #2 + //std::vector indices = {3, 5, 8, 9, 10, 11, 12, 13, 14}; // #3 + //std::vector indices = {2,8,9,4,10,11,5,12,13,7,14}; // #4 + + // Write errors + try { + std::string outputFileName = "Errors.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : errors) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write xn + try { + std::string outputFileName = "x.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : xn) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write dxn/dt + try { + std::string outputFileName = "dxn_dt.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : xn_dt) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write d2xn/dt2 + try { + std::string outputFileName = "d2xn_dt2.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : xn_dt2) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write lambda + try { + std::string outputFileName = "lambda.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : lambdan) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + // Write alpha + try { + std::string outputFileName = "alpha.csv"; + std::cout << "Generating CSV file: " << outputFileName << std::endl; + + std::ofstream outputFile; + outputFile.open(outputFileName); + + int it = 1; + + for (const auto& vec : Alphan) { + outputFile << std::setprecision(4) << it << ","; + for (const auto& index : indices) { + outputFile << vec[index] << ","; + } + outputFile << vec[13]; + outputFile << "\n"; + it++; + } + // Close the file + outputFile.close(); + + } + catch(...) { + throw std::invalid_argument("Couldn't write csv file"); + } + + std::cout << "CSV file has been generated " << std::endl; +} + +} // namespace nodal \ No newline at end of file diff --git a/src/porting/jsonPorter.hh b/src/porting/jsonPorter.hh index 96d3cc9..7402c4f 100644 --- a/src/porting/jsonPorter.hh +++ b/src/porting/jsonPorter.hh @@ -107,12 +107,17 @@ void simulationFromJSON(json jsonString, arch::Network* network_, sim::Simula // Read a Hybrid simulation definition else if (simType == sim::Type::Hybrid) { + std::cout << "[readContinuousPase]" << std::endl; + readContinuousPhase(jsonString, simulation, activeFixture); readResistanceModel(jsonString, simulation); if (platform == sim::Platform::Continuous) { + std::cout << "[readSimulators]" << std::endl; readSimulators(jsonString, simulation, network_); + std::cout << "[readUpdateScheme]" << std::endl; readUpdateScheme(jsonString, simulation); + std::cout << "[sortGroups]" << std::endl; network_->sortGroups(); } else if (platform == sim::Platform::Mixing) { readMixingModel(jsonString, simulation); @@ -136,8 +141,9 @@ void simulationFromJSON(json jsonString, arch::Network* network_, sim::Simula } else { throw std::invalid_argument("Invalid platform for Hybrid simulation. Please select one of the following:\n\tcontinuous"); } - + std::cout << "[readPumps]" << std::endl; readPumps(jsonString, network_); + std::cout << "[isValid]" << std::endl; network_->isNetworkValid(); } diff --git a/src/porting/jsonReaders.hh b/src/porting/jsonReaders.hh index 29cf8b8..47316cb 100644 --- a/src/porting/jsonReaders.hh +++ b/src/porting/jsonReaders.hh @@ -334,13 +334,12 @@ void readUpdateScheme(json jsonString, sim::Simulation& simulation) { if (simulator.contains("alpha") && simulator.contains("beta") && simulator.contains("theta")) { std::cout << "am here" << std::endl; if (simulator["alpha"].is_number() && simulator["beta"].is_number()) { - std::cout << "doing this 1" << std::endl; T alpha = simulator["alpha"]; T beta = simulator["beta"]; int theta = simulator["theta"]; - simulation.setNaiveHybridScheme(moduleCounter, alpha, beta, theta); + //simulation.setNaiveHybridScheme(moduleCounter, alpha, beta, theta); + simulation.setNaiveHybridScheme(alpha, beta, theta); } else if (simulator["alpha"].is_array() && simulator["beta"].is_array()) { - std::cout << "doing this 2" << std::endl; int nodeCounter = 0; std::unordered_map alpha; std::unordered_map beta; @@ -362,6 +361,14 @@ void readUpdateScheme(json jsonString, sim::Simulation& simulation) { } } } + } else if (jsonString["simulation"]["updateScheme"]["scheme"] == "DynamicDamping") { + //int theta = jsonString["simulation"]["updateScheme"]["theta"]; + int theta = 1; + simulation.setDynamicHybridScheme(theta); + return; + } else if (jsonString["simulation"]["updateScheme"]["scheme"] == "LinearDecoupling") { + simulation.setLinearDecouplingScheme(); + return; } else { throw std::invalid_argument("Update scheme was not set."); } @@ -389,9 +396,14 @@ void readPumps(json jsonString, arch::Network* network) { if (!jsonString["simulation"].contains("pumps") || jsonString["simulation"]["pumps"].empty()) { throw std::invalid_argument("No pumps are defined. Please define at least 1 pump."); } + for (auto& [key,channel] : network->getChannels()) { + std::cout << "key:\t" << key << "\t\tchannelId:\t" << channel->getId() << std::endl; + } + std::cout << "channels size: " << network->getChannels().size() << std::endl; for (auto& pump : jsonString["simulation"]["pumps"]) { if (pump.contains("channel") && pump.contains("type")) { int channelId = pump["channel"]; + std::cout << channelId << std::endl; if (pump["type"] == "PumpPressure") { if (pump.contains("deltaP")) { T pressure = pump["deltaP"]; diff --git a/src/simulation/Simulation.h b/src/simulation/Simulation.h index 2319a84..51fded8 100644 --- a/src/simulation/Simulation.h +++ b/src/simulation/Simulation.h @@ -33,6 +33,12 @@ class Opening; namespace mmft { // Forward declared dependencies +template +class DynamicDampingScheme; + +template +class LinearDecouplingScheme; + template class Scheme; @@ -143,7 +149,7 @@ class Simulation { ResistanceModel* resistanceModel; ///< The resistance model used for the simulation. MembraneModel* membraneModel; ///< The membrane model used for an OoC simulation. MixingModel* mixingModel; ///< The mixing model used for a mixing simulation. - std::unordered_map>> updateSchemes; ///< The update scheme for Abstract-CFD coupling + std::shared_ptr> updateScheme; ///< The update scheme for Abstract-CFD coupling int continuousPhase = 0; ///< Fluid of the continuous phase. int iteration = 0; int maxIterations = 1e5; @@ -305,16 +311,20 @@ class Simulation { * @param[in] theta The amount of LBM stream and collide cycles between updates for the modules. * @returns A shared_ptr to the created naive update scheme. */ - std::shared_ptr> setNaiveHybridScheme(int moduleId, T alpha, T beta, int theta); + std::shared_ptr> setNaiveHybridScheme(int moduleId, std::unordered_map alpha, std::unordered_map beta, int theta); /** - * @brief Define and set the naive update scheme for a hybrid simulation. - * @param[in] alpha The relaxation value for the pressure value update for all nodes of the module. - * @param[in] beta The relaxation value for the flow rate value update for all nodes of the module. - * @param[in] theta The amount of LBM stream and collide cycles between updates for the modules. - * @returns A shared_ptr to the created naive update scheme. + * @brief Define and set the dynamic damping update scheme for a hybrid simulation. + * @param[in] theta The amount of LBM stream and collide cycles between updates for all modules. + * @returns A shared_ptr to the created dynamic damping update scheme. */ - std::shared_ptr> setNaiveHybridScheme(int moduleId, std::unordered_map alpha, std::unordered_map beta, int theta); + std::shared_ptr> setDynamicHybridScheme(int theta); + + /** + * @brief Define and set the dynamic damping update scheme for a hybrid simulation. + * @returns A shared_ptr to the created dynamic damping update scheme. + */ + std::shared_ptr> setLinearDecouplingScheme(); /** * @brief Create injection. @@ -629,6 +639,10 @@ class Simulation { * @brief Print the results as pressure at the nodes and flow rates at the channels */ void printResults(); + + std::shared_ptr> getNodalAnalysis() { + return nodalAnalysis; + } }; } // namespace sim diff --git a/src/simulation/Simulation.hh b/src/simulation/Simulation.hh index e6b6b82..33dd1b9 100644 --- a/src/simulation/Simulation.hh +++ b/src/simulation/Simulation.hh @@ -201,27 +201,29 @@ namespace sim { template std::shared_ptr> Simulation::setNaiveHybridScheme(T alpha, T beta, int theta) { auto naiveScheme = std::make_shared>(network->getModules(), alpha, beta, theta); - for (auto& [key, simulator] : cfdSimulators) { - updateSchemes.try_emplace(simulator->getId(), naiveScheme); - simulator->setUpdateScheme(updateSchemes.at(simulator->getId())); - } + updateScheme = naiveScheme; return naiveScheme; } template - std::shared_ptr> Simulation::setNaiveHybridScheme(int moduleId, T alpha, T beta, int theta) { + std::shared_ptr> Simulation::setNaiveHybridScheme(int moduleId, std::unordered_map alpha, std::unordered_map beta, int theta) { auto naiveScheme = std::make_shared>(network->getModule(moduleId), alpha, beta, theta); - updateSchemes.try_emplace(moduleId, naiveScheme); - cfdSimulators.at(moduleId)->setUpdateScheme(updateSchemes.at(moduleId)); + updateScheme = naiveScheme; return naiveScheme; } template - std::shared_ptr> Simulation::setNaiveHybridScheme(int moduleId, std::unordered_map alpha, std::unordered_map beta, int theta) { - auto naiveScheme = std::make_shared>(network->getModule(moduleId), alpha, beta, theta); - updateSchemes.try_emplace(moduleId, naiveScheme); - cfdSimulators.at(moduleId)->setUpdateScheme(updateSchemes.at(moduleId)); - return naiveScheme; + std::shared_ptr> Simulation::setDynamicHybridScheme(int theta) { + auto dynamicScheme = std::make_shared>(theta, 10); + updateScheme = dynamicScheme; + return dynamicScheme; + } + + template + std::shared_ptr> Simulation::setLinearDecouplingScheme() { + auto linearDecouplingScheme = std::make_shared>(); + updateScheme = linearDecouplingScheme; + return linearDecouplingScheme; } template @@ -589,8 +591,7 @@ namespace sim { // conduct CFD simulations allConverged = conductCFDSimulation(cfdSimulators); // compute nodal analysis again - pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators); - + pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators, updateScheme); } #ifdef VERBOSE @@ -626,7 +627,7 @@ namespace sim { // conduct CFD simulations allConverged = conductCFDSimulation(cfdSimulators); // compute nodal analysis again - pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators); + pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators,updateScheme); } #ifdef VERBOSE @@ -666,7 +667,7 @@ namespace sim { // conduct CFD simulations allConverged = conductCFDSimulation(cfdSimulators); // compute nodal analysis again - pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators); + pressureConverged = nodalAnalysis->conductNodalAnalysis(cfdSimulators, updateScheme); } #ifdef VERBOSE @@ -890,7 +891,7 @@ namespace sim { #ifdef VERBOSE std::cout << "[Simulation] Conduct initial nodal analysis..." << std::endl; #endif - nodalAnalysis->conductNodalAnalysis(cfdSimulators); + nodalAnalysis->conductNodalAnalysis(cfdSimulators, updateScheme); // Prepare CFD geometry and lattice #ifdef VERBOSE @@ -919,7 +920,7 @@ namespace sim { #ifdef VERBOSE std::cout << "[Simulation] Conduct initial nodal analysis..." << std::endl; #endif - nodalAnalysis->conductNodalAnalysis(cfdSimulators); + nodalAnalysis->conductNodalAnalysis(cfdSimulators, updateScheme); // Prepare CFD geometry and lattice #ifdef VERBOSE @@ -948,7 +949,7 @@ namespace sim { #ifdef VERBOSE std::cout << "[Simulation] Conduct initial nodal analysis..." << std::endl; #endif - nodalAnalysis->conductNodalAnalysis(cfdSimulators); + nodalAnalysis->conductNodalAnalysis(cfdSimulators, updateScheme); // Prepare CFD geometry and lattice #ifdef VERBOSE diff --git a/src/simulation/simulators/cfdSimulator.h b/src/simulation/simulators/cfdSimulator.h index fadb0ce..f7aeb7c 100755 --- a/src/simulation/simulators/cfdSimulator.h +++ b/src/simulation/simulators/cfdSimulator.h @@ -48,7 +48,6 @@ class CFDSimulator { std::shared_ptr> moduleNetwork; ///< Fully connected graph as network for the initial approximation. std::unordered_map> moduleOpenings; ///< Map of openings. std::unordered_map groundNodes; ///< Map of nodes that communicate the pressure to the 1D solver. - std::shared_ptr> updateScheme; ///< The update scheme for Abstract-CFD coupling /** * @brief Constructor of a CFDSimulator, which acts as a base definition for other simulators. diff --git a/src/simulation/simulators/cfdSimulator.hh b/src/simulation/simulators/cfdSimulator.hh index 552436d..8309f5f 100755 --- a/src/simulation/simulators/cfdSimulator.hh +++ b/src/simulation/simulators/cfdSimulator.hh @@ -27,7 +27,7 @@ template CFDSimulator::CFDSimulator (int id_, std::string name_, std::string stlFile_, std::shared_ptr> cfdModule_, std::unordered_map> openings_, std::shared_ptr> updateScheme_, ResistanceModel* resistanceModel_) : CFDSimulator (id_, name_, stlFile_, cfdModule_, openings_, resistanceModel_) { - updateScheme = updateScheme_; + //updateScheme = updateScheme_; } template @@ -58,7 +58,7 @@ void CFDSimulator::setInitialized(bool initialization_) { template void CFDSimulator::setUpdateScheme(const std::shared_ptr>& updateScheme_) { - this->updateScheme = updateScheme_; + //this->updateScheme = updateScheme_; } template @@ -73,12 +73,12 @@ std::string CFDSimulator::getVtkFile() { template T CFDSimulator::getAlpha(int nodeId_) { - return updateScheme->getAlpha(nodeId_); + //return updateScheme->getAlpha(nodeId_); } template T CFDSimulator::getBeta(int nodeId_) { - return updateScheme->getBeta(nodeId_); + //return updateScheme->getBeta(nodeId_); } template diff --git a/src/simulation/simulators/olbContinuous.hh b/src/simulation/simulators/olbContinuous.hh index 46b0370..b44087c 100644 --- a/src/simulation/simulators/olbContinuous.hh +++ b/src/simulation/simulators/olbContinuous.hh @@ -99,7 +99,7 @@ void lbmSimulator::writeVTK (int iT) { this->vtkFile = olb::singleton::directories().getVtkOutDir() + olb::createFileName( this->name ) + ".pvd"; } - if (iT % 1000 == 0) { + if (iT % 10000 == 0) { olb::SuperLatticePhysVelocity2D velocity(getLattice(), getConverter()); olb::SuperLatticePhysPressure2D pressure(getLattice(), getConverter()); @@ -116,13 +116,14 @@ void lbmSimulator::writeVTK (int iT) { if (iT %1000 == 0) { #ifdef VERBOSE std::cout << "[writeVTK] " << this->name << " currently at timestep " << iT << std::endl; + /* for (auto& [key, Opening] : this->moduleOpenings) { if (this->groundNodes.at(key)) { meanPressures.at(key)->print(); } else { fluxes.at(key)->print(); } - } + } */ #endif } @@ -138,7 +139,8 @@ void lbmSimulator::writeVTK (int iT) { template void lbmSimulator::solve() { - int theta = this->updateScheme->getTheta(this->cfdModule->getId()); + //int theta = this->updateScheme->getTheta(this->cfdModule->getId()); + int theta = 1; this->setBoundaryValues(step); for (int iT = 0; iT < theta; ++iT){ writeVTK(step);