WIP: Define Diffusive Mixing Coupling

cda-tum · Nov 26, 2024 · d815d98 · d815d98
1 parent 9ed5daf
commit d815d98
Show file tree

Hide file tree

Showing 14 changed files with 184 additions and 37 deletions.
diff --git a/src/olbProcessors/AdeConcBoundaries2D.h b/src/olbProcessors/AdeConcBoundaries2D.h
@@ -0,0 +1,57 @@
+#pragma once
+
+#include <olb2D.h>
+#include <olb2D.hh>
+
+namespace arch {
+
+// Forward declared dependencies
+template<typename T>
+class Opening;
+
+template <typename T>
+class Node;
+
+}
+
+namespace olb {
+
+template<typename T>
+class AdeConcBC : public AnalyticalF2D<T,T> {
+
+private:
+  T width;                ///< Physical width of the opening.
+  T xc;                   ///< The physical x coordinate, xc, of the opening node (center).
+  T yc;                   ///< The physical y coordinate, yc, of the opening node (center).
+  T x0;                   ///< The physical x coordinate, x0, at which the opening "starts" (bottom for right facing opening).
+  T y0;                   ///< The physical y coordinate, y0, at which the opening "starts" (bottom for right facing opening).
+  std::vector<T> array;   ///< Vector that contains the array of concentration values along the opening.
+
+public:
+
+  AdeConcBC(const arch::Opening<T>& opening, std::vector<T> concentrationField);
+
+  bool operator()(T output[], const T input[]) override;
+
+};
+
+template<typename T>
+class AdeConc1D : public SuperPlaneIntegralF2D<T,T> {
+
+private:
+  T width;                ///< Physical width of the opening.
+  T xc;                   ///< The physical x coordinate, xc, of the opening node (center).
+  T yc;                   ///< The physical y coordinate, yc, of the opening node (center).
+  T x0;                   ///< The physical x coordinate, x0, at which the opening "starts" (bottom for right facing opening).
+  T y0;                   ///< The physical y coordinate, y0, at which the opening "starts" (bottom for right facing opening).
+  std::vector<T> array;   ///< Vector that contains the array of concentration values along the opening.
+
+public:
+
+  AdeConc1D(const arch::Opening<T>& opening, std::vector<T> concentrationField);
+
+  bool operator()(T output[], const T input[]) override;
+
+};
+
+}
diff --git a/src/olbProcessors/AdeConcBoundaries2D.hh b/src/olbProcessors/AdeConcBoundaries2D.hh
@@ -0,0 +1,37 @@
+#include <AdeConcBoundary2D.h>
+
+namespace olb {
+
+template<typename T>
+AdeConcBoundary2D<T>::AdeConcBoundary2D(const arch::Opening<T>& opening, std::vector<T> concentrationField) :
+  width(opening.width),
+  xc(opening.node->getPosition()[0]),
+  yc(opening.node->getPosition()[0]),
+  x0(xc - opening.tangent[0]*0.5*width),
+  y0(yc - opening.tangent[1]*0.5*width),
+  array(concentrationField) { }
+
+template<typename T>
+bool AdeConcBoundary2D<T>::operator()(T output[], const T input[]) {
+  T s   = std::sqrt((input[0]-x0)*(input[0]-x0) + (input[1]-y0)*(input[1]-y0));
+  T d   = std::sqrt((input[0]-xc)*(input[0]-xc) + (input[1]-yc)*(input[1]-yc));
+  T res = width / array.size();
+  if (s < width) {
+    output[0] = array[std::floor(s/res)];
+  }
+  if (0.5*width - d < 0.0) {
+    output[0] = T();
+  }
+  return true;
+}
+
+template<typename T>
+AdeConc1D<T>::AdeConc1D(const arch::Opening<T>& opening, std::vector<T> concentrationField) 
+  : SuperPlaneIntegralF2D<T>()
+
+template<typename T>
+bool AdeConc1D<T>::operator()(T output[], const T input[]) {
+
+}
+
+} // namespace olb
diff --git a/src/simulation/CFDSim.hh b/src/simulation/CFDSim.hh
@@ -47,7 +47,7 @@ namespace sim {
     }
 
     template<typename T>
-    bool conductADSimulation(const std::unordered_map<int, std::unique_ptr<CFDSimulator<T>>>& cfdSimulators) {
+    bool conductADSimulation(const std::unordered_map<int, std::unique_ptr<CFDSimulator<T>>>& cfdSimulators, bool fieldValues) {
 
         bool allConverge = true;
 
@@ -61,7 +61,7 @@ namespace sim {
             assert(cfdSimulator.second->getModule()->getModuleType() == arch::ModuleType::ESS_LBM);
             throw std::runtime_error("Simulation of Advection Diffusion not defined for ESS LBM.");
             #endif
-            cfdSimulator.second->adSolve();
+            cfdSimulator.second->adSolve(fieldValues);
 
             if (!cfdSimulator.second->hasAdConverged()) {
                 allConverge = false;

diff --git a/src/simulation/MixingModels.hh b/src/simulation/MixingModels.hh
@@ -469,24 +469,30 @@ void DiffusionMixingModel<T>::generateInflows(T timeStep, arch::Network<T>* netw
             for (auto& [channelId, channel] : network->getChannels()) { // There is only one channel connected to a node at a CFD hybrid node connection
                 if (cfdSimulator->getFlowRates().at(nodeId) > 0.0 && (channel->getNodeA() == nodeId || channel->getNodeB() == nodeId)) {
                     std::unordered_map<int, std::tuple<std::function<T(T)>,std::vector<T>,T>> newDistributions;
-                    std::cout << "Getting here, nodeId is " + std::to_string(nodeId) << std::endl;
+                    std::cout << "Getting here 1, nodeId is " + std::to_string(nodeId) << std::endl;
                     for (const auto& [specieId, concentrationField] : concentrationFieldsIn.at(nodeId)) {
-                        std::cout << "Getting here 2" << std::endl;
+                        std::cout << "Getting here 1.1" << std::endl;
                         T dx = channel->getWidth() / concentrationField.size();
+                        std::cout << "Getting here 1.2" << std::endl;
                         T pecletNr = (std::abs(channel->getFlowRate()) / channel->getHeight()) / (sim->getSpecie(specieId))->getDiffusivity();
+                        std::cout << "Getting here 1.3" << std::endl;
                         std::tuple<std::function<T(T)>, std::vector<T>, T> analyticalResult = getAnalyticalSolutionHybrid(
                             channel->getLength(), channel->getFlowRate(), channel->getWidth(), noFourierTerms, pecletNr, concentrationField, dx);
+                        std::cout << "Getting here 1.4" << std::endl;
                         newDistributions.try_emplace(specieId, analyticalResult);
+                        std::cout << "Getting here 1.5" << std::endl;
                     }
                     //Create new DiffusiveMixture
                     DiffusiveMixture<T>* newMixture = dynamic_cast<DiffusiveMixture<T>*>(sim->addDiffusiveMixture(newDistributions));
                     newMixture->setNonConstant();
                     this->injectMixtureInEdge(newMixture->getId(), channelId);
                 }
                 else if (cfdSimulator->getFlowRates().at(nodeId) < 0.0 && (channel->getNodeA() == nodeId || channel->getNodeB() == nodeId)) {
+                    std::cout << "Getting here 2, nodeId is " + std::to_string(nodeId) << std::endl;
                     int resolutionIntoCfd = cfdSimulator->getResolution(nodeId);
                     concentrationFieldsOut.try_emplace(nodeId, std::unordered_map<int, std::vector<T>>());
                     for (const auto& [specieId, speciePtr] : sim->getSpecies()) {
+                        std::cout << "Getting here 2.1" << std::endl;
                         std::vector<T> concentrationFieldForCfdInput = defineConcentrationNodeFieldForCfdInput(resolutionIntoCfd, specieId, channelId, channel->getWidth(), mixtures, noFourierTerms);
                         concentrationFieldsOut.at(nodeId).try_emplace(specieId, concentrationFieldForCfdInput);
                     }                    

diff --git a/src/simulation/Simulation.hh b/src/simulation/Simulation.hh
@@ -639,7 +639,7 @@ namespace sim {
             // Obtain overal steady-state concentration results
             bool concentrationConverged = false;
             while (!concentrationConverged) {
-                concentrationConverged = conductADSimulation(cfdSimulators);
+                concentrationConverged = conductADSimulation(cfdSimulators, mixingModel->isDiffusive());
                 this->mixingModel->propagateSpecies(network, this);
             }
         }
@@ -682,7 +682,7 @@ namespace sim {
             // Obtain overal steady-state concentration results
             bool concentrationConverged = false;
             while (!concentrationConverged) {
-                concentrationConverged = conductADSimulation(cfdSimulators);
+                concentrationConverged = conductADSimulation(cfdSimulators, mixingModel->isDiffusive());
                 this->mixingModel->propagateSpecies(network, this);
             }
         }

diff --git a/src/simulation/simulators/cfdSimulator.h b/src/simulation/simulators/cfdSimulator.h
@@ -256,7 +256,7 @@ class CFDSimulator {
      * @brief Set the boundary values on the lattice at the module nodes.
      * @param[in] iT Iteration step.
     */
-    virtual void setBoundaryValues(int iT) = 0;
+    virtual void setBoundaryValues(int iT, bool fieldValues=false) = 0;
 
     /**
      * @brief Prepare the LBM geometry of this simulator.
@@ -285,7 +285,7 @@ class CFDSimulator {
     /**
      * @brief Conducts the collide and stream operations of the AD lattice(s).
     */
-    virtual void adSolve() 
+    virtual void adSolve(bool fieldValues) 
     {
         throw std::runtime_error("The function prepareLattice is undefined for this CFD simulator.");
     }

diff --git a/src/simulation/simulators/essContinuous.h b/src/simulation/simulators/essContinuous.h
@@ -83,7 +83,7 @@ namespace sim {
              * @brief Set the boundary values on the lattice at the module nodes.
              * @param[in] iT Iteration step.
             */
-            void setBoundaryValues(int iT);
+            void setBoundaryValues(int iT, bool fieldValues=false);
 
             /**
              * @brief Conducts the collide and stream operations of the lattice.

diff --git a/src/simulation/simulators/essContinuous.hh b/src/simulation/simulators/essContinuous.hh
@@ -17,7 +17,7 @@ namespace sim{
     }
 
     template<typename T>
-    void essLbmSimulator<T>::setBoundaryValues(int iT)
+    void essLbmSimulator<T>::setBoundaryValues(int iT, bool fieldValues)
     {
         setFlowRates(flowRates);
         setPressures(pressures);

diff --git a/src/simulation/simulators/essMixing.h b/src/simulation/simulators/essMixing.h
@@ -83,7 +83,7 @@ namespace sim {
              * @brief Set the boundary values on the lattice at the module nodes.
              * @param[in] iT Iteration step.
             */
-            void setBoundaryValues(int iT);
+            void setBoundaryValues(int iT, bool fieldValues=false);
 
             /**
              * @brief Conducts the collide and stream operations of the lattice.

diff --git a/src/simulation/simulators/essMixing.hh b/src/simulation/simulators/essMixing.hh
@@ -18,7 +18,7 @@ namespace sim{
     }
 
     template<typename T>
-    void essLbmSimulator<T>::setBoundaryValues(int iT)
+    void essLbmSimulator<T>::setBoundaryValues(int iT, bool fieldValues)
     {
         setFlowRates(flowRates);
         setPressures(pressures);

diff --git a/src/simulation/simulators/olbContinuous.h b/src/simulation/simulators/olbContinuous.h
@@ -174,7 +174,7 @@ using BounceBack = olb::BounceBack<T,DESCRIPTOR>;
      * @brief Set the boundary values on the lattice at the module nodes.
      * @param[in] iT Iteration step.
     */
-    void setBoundaryValues(int iT) override;
+    void setBoundaryValues(int iT, bool fieldValues=false) override;
 
     /**
      * @brief Conducts the collide and stream operations of the lattice.

diff --git a/src/simulation/simulators/olbContinuous.hh b/src/simulation/simulators/olbContinuous.hh
@@ -74,7 +74,7 @@ void lbmSimulator<T>::prepareLattice () {
 }
 
 template<typename T>
-void lbmSimulator<T>::setBoundaryValues (int iT) {
+void lbmSimulator<T>::setBoundaryValues (int iT, bool fieldValues) {
 
     for (auto& [key, Opening] : this->moduleOpenings) {
         if (this->groundNodes.at(key)) {
@@ -389,15 +389,15 @@ template<typename T>
 void lbmSimulator<T>::setFlowProfile2D (int openingKey, T openingWidth)  {
     T maxVelocity = (3./2.)*(flowRates[openingKey]/(openingWidth));
     T distance2Wall = getConverter().getConversionFactorLength()/2.;
-    this->flowProfiles.at(openingKey) = std::make_shared<olb::Poiseuille2D<T>>(getGeometry(), openingKey+3, getConverter().getLatticeVelocity(maxVelocity), distance2Wall);
-    getLattice().defineU(getGeometry(), openingKey+3, *this->flowProfiles.at(openingKey));
+    flowProfiles.at(openingKey) = std::make_shared<olb::Poiseuille2D<T>>(getGeometry(), openingKey+3, getConverter().getLatticeVelocity(maxVelocity), distance2Wall);
+    getLattice().defineU(getGeometry(), openingKey+3, *flowProfiles.at(openingKey));
 }
 
 template<typename T>
 void lbmSimulator<T>::setPressure2D (int openingKey)  {
     T rhoV = getConverter().getLatticeDensityFromPhysPressure((pressures[openingKey]));
-    this->densities.at(openingKey) = std::make_shared<olb::AnalyticalConst2D<T,T>>(rhoV);
-    getLattice().defineRho(getGeometry(), openingKey+3, *this->densities.at(openingKey));
+    densities.at(openingKey) = std::make_shared<olb::AnalyticalConst2D<T,T>>(rhoV);
+    getLattice().defineRho(getGeometry(), openingKey+3, *densities.at(openingKey));
 }
 
 template<typename T>

diff --git a/src/simulation/simulators/olbMixing.h b/src/simulation/simulators/olbMixing.h
@@ -44,8 +44,8 @@ using ADDynamics = olb::AdvectionDiffusionBGKdynamics<T,ADDESCRIPTOR>;
 using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
 
 protected:
-    std::unordered_map<int, std::unordered_map<int, T>> concentrations;   ///< Vector of concentration values at module nodes. <nodeId, <speciId, conc>>
-    std::unordered_map<int, std::unordered_map<int, std::vector<T>>> nodeConcentrationFields;   ///< Vector of concentration values at module nodes. <nodeId, <speciId, <conc>>>
+    std::unordered_map<int, std::unordered_map<int, T>> concentrations;   ///< Vector of concentration values at module nodes. <nodeId, <specieId, conc>>
+    std::unordered_map<int, std::unordered_map<int, std::vector<T>>> nodeConcentrationFields;   ///< Vector of concentration values at module nodes. <nodeId, <specieId, <conc>>>
 
     std::unordered_map<int, Specie<T>*> species;
 
@@ -61,7 +61,7 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
     std::unordered_map<int, T*> fluxWall;
     T zeroFlux = 0.0;
 
-    std::unordered_map<int, std::unordered_map<int, std::shared_ptr<olb::AnalyticalConst2D<T,T>>>> concentrationProfiles;
+    std::unordered_map<int, std::unordered_map<int, std::shared_ptr<olb::AdeConcBoundary2D<T,T>>>> concentrationProfiles;   ///< Map of concentration field in olb format at module nodes. <nodeId, <specieId, <conc>>>
     std::unordered_map<int, std::unordered_map<int, std::shared_ptr<olb::SuperPlaneIntegralFluxPressure2D<T>>>> meanConcentrations;       ///< Map of mean pressure values at module nodes.
 
     auto& getAdConverter(int key) {
@@ -94,13 +94,15 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
     */
     void executeCoupling() override;
 
-    void setConcentration2D(int key);
+    void setConcentrationConst2D(int key);
+
+    void setConcentrationField2D(int key);
 
     /**
      * @brief Update the values at the module nodes based on the simulation result after stepIter iterations.
      * @param[in] iT Iteration step.
     */
-    void storeCfdResults(int iT);
+    void storeCfdResults(int iT, bool fieldValues);
 
 public:
     /**
@@ -159,7 +161,7 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
      * @brief Set the boundary values on the lattice at the module nodes.
      * @param[in] iT Iteration step.
     */
-    void setBoundaryValues(int iT) override;
+    void setBoundaryValues(int iT, bool fieldValues) override;
 
     /**
      * @brief Conducts the collide and stream operations of the lattice.
@@ -174,7 +176,7 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
     /**
      * @brief Conducts the collide and stream operations of the AD lattice(s).
     */
-    void adSolve() override;
+    void adSolve(bool fieldValues) override;
 
     /**
      * @brief Write the vtk file with results of the CFD simulation to file system.