🚧 Add Abstract and CFD coupling for mixtures

examples/Hybrid/Mixing1a.JSON

@@ -174,7 +174,7 @@
             "mixtureInjections": [
                 {
                     "mixture": 0,
-                    "channel": 1,
+                    "channel": 4,
                     "t0": 0.0
                 }
             ]

src/simulation/MixingModels.h

@@ -212,6 +212,11 @@ class InstantaneousMixingModel : public MixingModel<T> {
     */
     void clean(arch::Network<T>* network);
 
+    /**
+     * @brief Propagate all the species through a network for a steady-state simulation
+     */
+    void propagateSpecies(arch::Network<T>* network, Simulation<T>* sim) override;
+
     /**
      * @brief From the mixtureInjections and CFD simulators, generate temporary mxtures that 
      * flow into the network at correspondingnode entry points.
@@ -220,21 +225,16 @@ class InstantaneousMixingModel : public MixingModel<T> {
     void initNodeOutflow(Simulation<T>* sim, std::vector<Mixture<T>>& tmpMixtures);
 
     /**
-     * @brief Propagate all the species through a network for a steady-state simulation
+     * @brief Propagate the mixtures through the corresponding channel entirely, without considering time steps
      */
-    void propagateSpecies(arch::Network<T>* network, Simulation<T>* sim) override;
+    void channelPropagation(arch::Network<T>* network);
 
     /**
      * @brief From the node's inflows, generate the node outflow
      */
-    void updateNodeInflow(arch::Network<T>* network);
-
-    /**
-     * @brief Propagate the mixtures through the corresponding channel entirely, without considering time steps
-     */
-    bool channelPropagation(arch::Network<T>* network);
+    bool updateNodeOutflow(Simulation<T>* sim, std::vector<Mixture<T>>& tmpMixtures);
 
-    void storeConcentrations(Simulation<T>* sim);
+    void storeConcentrations(Simulation<T>* sim, const std::vector<Mixture<T>>& tmpMixtures);
 
     /**
      * @brief Print all mixtures and their positions in the network.

src/simulation/MixingModels.hh

@@ -68,28 +68,50 @@ void InstantaneousMixingModel<T>::propagateSpecies(arch::Network<T>* network, Si
 
     std::vector<Mixture<T>> tmpMixtures;
 
+    std::cout << "[propagateSpecies] Define inflow volume" << std::endl;
+
+    // Define total inflow volume at nodes
+    for (auto& [nodeId, node] : network->getNodes()) {
+        for (auto& channel : network->getChannelsAtNode(nodeId) ) {
+            // Check if the channel flows into the node
+            if ((channel->getFlowRate() > 0.0 && channel->getNodeB() == nodeId) || (channel->getFlowRate() < 0.0 && channel->getNodeA() == nodeId)) {
+                T inflowVolume = std::abs(channel->getFlowRate());
+                auto [iterator, inserted] = totalInflowVolumeAtNode.try_emplace(nodeId, inflowVolume);
+                if (!inserted) {
+                    iterator->second += inflowVolume;
+                }
+            }
+        }
+    }
+
+    std::cout << "[propagateSpecies] Initialize the node inflow" << std::endl;
     // Initial node outflow from mixtureInjections and CFD simulators, stored in mixtureOutflowAtNode
     initNodeOutflow(sim, tmpMixtures);
 
-    // Propagate the mixtures through the entire channel, without considering time steps
-    bool networkCleared = channelPropagation(network);
+    std::cout << "[propagateSpecies] Channel propagation" << std::endl;
+    // Propagate the mixtures through the entire channel, without considering time steps    
+    channelPropagation(network);
 
-    while (!networkCleared) {
+    bool inflowUpdated = true;
+    while (inflowUpdated) {
+        std::cout << "[propagateSpecies] Update node outflow" << std::endl;
         // From node inflow, generate the node's outflow
-        updateNodeInflow(network);
+        inflowUpdated = updateNodeOutflow(sim, tmpMixtures);
         // Propagate the mixtures through the entire channel
-        networkCleared = channelPropagation(network);
+        std::cout << "[propagateSpecies] Channel propagation" << std::endl;
+        channelPropagation(network);
     }
 
     // Store the concentrations of the final state in the concentration buffer of olbMixingSolver.
-    storeConcentrations(sim);
+    storeConcentrations(sim, tmpMixtures);
 
     clean(network);
 
 }
 
 template<typename T>
 void InstantaneousMixingModel<T>::initNodeOutflow(Simulation<T>* sim, std::vector<Mixture<T>>& tmpMixtures) {
+    std::cout << "[initNodeOutflow] Add mixture injections" << std::endl;
     // Add mixture injections
     for (auto& [key, mixtureInjection] : sim->getMixtureInjections()) {
         int tmpMixtureIndex = tmpMixtures.size();
@@ -98,7 +120,7 @@ void InstantaneousMixingModel<T>::initNodeOutflow(Simulation<T>* sim, std::vecto
         tmpMixtures.push_back(Mixture<T>(*sim->getMixture(mixtureInjection->getMixtureId())));
         mixtureOutflowAtNode.try_emplace(nodeId, tmpMixtureIndex);
     }
-
+    std::cout << "[initNodeOutflow] Add CFD outflows" << std::endl;
     // Add CFD Simulator outflows
     for (auto& [key, cfdSimulator] : sim->getCFDSimulators()) {
         for (auto& [nodeId, opening] : cfdSimulator->getOpenings()) {
@@ -107,11 +129,11 @@ void InstantaneousMixingModel<T>::initNodeOutflow(Simulation<T>* sim, std::vecto
                 int tmpMixtureIndex = tmpMixtures.size();
                 int tmpMixtureId = std::numeric_limits<int>::max();
                 std::unordered_map<int, Specie<T>*> species;
-                std::unordered_map<int, T> specieConcentrations;
-                /** TODO:
-                 * for specie in species, add species and specieConcentration to tmpMixture
-                 */
-                Mixture<T> tmpMixture = Mixture<T>(tmpMixtureId, species, specieConcentrations, sim->getContinuousPhase());
+                std::unordered_map<int, T> speciesConcentrations(cfdSimulator->getConcentrations().at(nodeId));
+                for (auto& [speciesId, concentration] : cfdSimulator->getConcentrations().at(nodeId)) {
+                    species.try_emplace(speciesId, sim->getSpecie(speciesId));
+                }
+                Mixture<T> tmpMixture = Mixture<T>(tmpMixtureId, species, speciesConcentrations, sim->getContinuousPhase());
                 tmpMixtures.push_back(tmpMixture);
                 mixtureOutflowAtNode.try_emplace(nodeId, tmpMixtureIndex);
             }
@@ -120,18 +142,101 @@ void InstantaneousMixingModel<T>::initNodeOutflow(Simulation<T>* sim, std::vecto
 }
 
 template<typename T>
-void InstantaneousMixingModel<T>::updateNodeInflow(arch::Network<T>* network) {
-    // TODO
+void InstantaneousMixingModel<T>::channelPropagation(arch::Network<T>* network) {
+    for (auto& [nodeId, mixtureId] : mixtureOutflowAtNode) {
+        for (auto& channel : network->getChannelsAtNode(nodeId)) {
+            // Find the nodeId that is across the channel
+            int oppositeNode;
+            if (channel->getFlowRate() > 0.0 && channel->getNodeA() == nodeId) {
+                oppositeNode = channel->getNodeB();
+            } else if (channel->getFlowRate() < 0.0 && channel->getNodeB() == nodeId) {
+                oppositeNode = channel->getNodeA();
+            } else {
+                continue;
+            }
+            // Update the mixture inflow at the node across the channel
+            MixtureInFlow<T> mixtureInflow = {mixtureId, std::abs(channel->getFlowRate())};
+            auto [iterator, inserted] = mixtureInflowAtNode.try_emplace(oppositeNode, std::vector<MixtureInFlow<T>>(1, mixtureInflow));
+            if (!inserted) {
+                mixtureInflowAtNode.at(oppositeNode).push_back(mixtureInflow);
+            }
+        }
+    }
 }
 
 template<typename T>
-bool InstantaneousMixingModel<T>::channelPropagation(arch::Network<T>* network) {
-    // TODO
+bool InstantaneousMixingModel<T>::updateNodeOutflow(Simulation<T>* sim, std::vector<Mixture<T>>& tmpMixtures) {
+    bool updated = false;
+    for (auto& [nodeId, mixtureInflowList] : mixtureInflowAtNode) {
+        bool createMixture = false;
+        std::unordered_map<int, Specie<T>*> speciePtrs;
+        std::unordered_map<int, T> newConcentrations;
+        for (auto& mixtureInflow : mixtureInflowList) {
+            for (auto& [specieId, oldConcentration] : tmpMixtures[mixtureInflow.mixtureId].getSpecieConcentrations()) {
+                speciePtrs.try_emplace(specieId, sim->getSpecie(specieId));
+                T newConcentration = oldConcentration * mixtureInflow.inflowVolume / totalInflowVolumeAtNode.at(nodeId);
+                auto [iterator, inserted] = newConcentrations.try_emplace(specieId, newConcentration);
+                if (!inserted) {
+                    iterator->second = iterator->second + newConcentration;
+                }
+            }
+            if (mixtureInflow.mixtureId != mixtureInflowList[0].mixtureId) {
+                createMixture = true;
+            }
+        }
+        int outflowMixtureId;
+        Mixture<T> newMixture (tmpMixtures.size(), speciePtrs, newConcentrations, sim->getContinuousPhase());
+        if (createMixture) {
+            outflowMixtureId = tmpMixtures.size();
+        } else {
+            outflowMixtureId = mixtureInflowList[0].mixtureId;
+        }
+        if (!mixtureOutflowAtNode.count(nodeId)) {
+            if (createMixture) {
+                tmpMixtures.push_back(newMixture);
+            }
+            mixtureOutflowAtNode.try_emplace(nodeId, outflowMixtureId);
+            updated = true;
+
+        } else {
+            // Check if the ouflow gets updated or not
+            if (tmpMixtures[mixtureOutflowAtNode.at(nodeId)] == tmpMixtures[outflowMixtureId]) {
+                continue;
+            } else {
+                if (createMixture) {
+                    tmpMixtures.push_back(newMixture);
+                }
+                mixtureOutflowAtNode.at(nodeId) = outflowMixtureId;
+                updated = true;
+            }
+        }
+    }
+    return updated;
 }
 
 template<typename T>
-void InstantaneousMixingModel<T>::storeConcentrations(Simulation<T>* sim) {
-    // TODO
+void InstantaneousMixingModel<T>::storeConcentrations(Simulation<T>* sim, const std::vector<Mixture<T>>& tmpMixtures) {
+    for (auto& [nodeId, mixtureId] : mixtureOutflowAtNode) {
+        std::cout << "nodeId: " << nodeId << "\tmixtureId: " << mixtureId << std::endl;
+    }
+    std:: cout << "There are " << tmpMixtures.size() << " mixtures in tmpMixtures" <<std::endl;
+    for (auto& [key, cfdSimulator] : sim->getCFDSimulators()) {
+        std::unordered_map<int, std::unordered_map<int, T>> concentrations = cfdSimulator->getConcentrations();
+        std::cout << "The inflow nodes are: ";
+        for (auto& [nodeId, opening] : cfdSimulator->getOpenings()) {
+            // If the node is an inflow
+            if (cfdSimulator->getFlowRates().at(nodeId) > 0.0) {
+                std::cout << nodeId << " ";
+                /*
+                for (auto& [specieId, specieConcentration] : tmpMixtures[mixtureOutflowAtNode.at(nodeId)].getSpecieConcentrations()) {
+                    concentrations.at(nodeId).at(specieId) = specieConcentration;
+                }
+                */
+            }
+        }
+        std::cout << std::endl;
+        cfdSimulator->storeConcentrations(concentrations);
+    }
 }
 
 template<typename T>

src/simulation/Mixture.h

@@ -90,6 +90,8 @@ class Mixture {
      */
     Mixture(int id, std::unordered_map<int, Specie<T>*> species, std::unordered_map<int, T> specieConcentrations, Fluid<T>* carrierFluid);
 
+    bool operator== (const Mixture<T> &t);
+
     /**
      * @brief Get the id of this mixture
      * @return Unique identifier of the mixture.

src/simulation/Mixture.hh

@@ -21,6 +21,17 @@ Mixture<T>::Mixture(int id, std::unordered_map<int, Specie<T>*> species, std::un
                     viscosity(carrierFluid->getViscosity()), density(carrierFluid->getDensity()), 
                     largestMolecularSize(0.0) { }
 
+template<typename T>
+bool Mixture<T>::operator== (const Mixture<T> &t) {
+    if (species == t.species && 
+        specieConcentrations == t.specieConcentrations && 
+        viscosity == t.viscosity &&
+        density == t.density) { 
+        return true; 
+    }
+    return false;
+}
+
 template<typename T>
 int Mixture<T>::getId() const {
     return id;

src/simulation/Simulation.hh

@@ -604,7 +604,9 @@ namespace sim {
             // Obtain overal steady-state concentration results
             bool concentrationConverged = false;
             while (!concentrationConverged) {
+                std::cout << "[Simulation] Conducting AD solve" << std::endl;
                 concentrationConverged = conductADSimulation(cfdSimulators);
+                std::cout << "[Simulation] Propagating Species" << std::endl;
                 this->mixingModel->propagateSpecies(network, this);
             }
         }

src/simulation/simulators/cfdSimulator.h

@@ -132,9 +132,9 @@ class CFDSimulator {
 
     virtual std::unordered_map<int, T> getFlowRates() const = 0;
 
-    virtual void storeConcentrations(std::unordered_map<int, std::vector<T>> concentrations) { }
+    virtual void storeConcentrations(std::unordered_map<int, std::unordered_map<int, T>> concentrations) { }
 
-    virtual std::unordered_map<int, std::vector<T>> getConcentrations() const { return std::unordered_map<int, std::vector<T>>(); }
+    virtual std::unordered_map<int, std::unordered_map<int, T>> getConcentrations() const { return std::unordered_map<int, std::unordered_map<int, T>>(); }
 
     virtual void setBoundaryValues(int iT) = 0;
 

src/simulation/simulators/olbMixing.h

@@ -162,13 +162,13 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
      * @brief Store the abstract concentrations at the nodes on the module boundary in the simulator.
      * @param[in] concentrations Map of concentrations and node ids.
      */
-    void storeConcentrations(std::unordered_map<int, std::vector<T>> concentrations) override;
+    void storeConcentrations(std::unordered_map<int, std::unordered_map<int, T>> concentrations) override;
 
     /**
      * @brief Get the concentrations at the boundary nodes.
      * @returns Concentrations
      */
-    std::unordered_map<int, std::vector<T>> getConcentrations() const override;
+    std::unordered_map<int, std::unordered_map<int, T>> getConcentrations() const override;
 
     /**
      * @brief Returns whether the module has converged or not.

src/simulation/simulators/olbMixing.hh

@@ -190,6 +190,7 @@ void lbmMixingSimulator<T>::solve() {
 template<typename T>
 void lbmMixingSimulator<T>::executeCoupling() {
     this->lattice->executeCoupling();
+    std::cout << "[lbmSimulator] Execute NS-AD coupling " << this->name << "... OK" << std::endl;
 }
 
 
@@ -225,7 +226,11 @@ void lbmMixingSimulator<T>::initValueContainers () {
     for (auto& [key, node] : this->moduleOpenings) {
         this->pressures.try_emplace(key, (T) 0.0);
         this->flowRates.try_emplace(key, (T) 0.0);
-        this->concentrations.try_emplace(key, std::vector<T>(species.size(), 0.0));
+        std::unordered_map<int, T> tmpConcentrations;
+        for (auto& [speciesId, speciesPtr] : species) {
+            tmpConcentrations.try_emplace(speciesId, 0.0);
+        }
+        this->concentrations.try_emplace(key, tmpConcentrations);
     }
 }
 
@@ -365,12 +370,12 @@ void lbmMixingSimulator<T>::setConcentration2D (int key) {
 }
 
 template<typename T>
-void lbmMixingSimulator<T>::storeConcentrations(std::unordered_map<int, std::vector<T>> concentrations_) {
+void lbmMixingSimulator<T>::storeConcentrations(std::unordered_map<int, std::unordered_map<int, T>> concentrations_) {
     this->concentrations = concentrations_;
 }
 
 template<typename T>
-std::unordered_map<int, std::vector<T>> lbmMixingSimulator<T>::getConcentrations() const {
+std::unordered_map<int, std::unordered_map<int, T>> lbmMixingSimulator<T>::getConcentrations() const {
     return this->concentrations;
 }