✨ Add Abstract and CFD coupling for mixtures

src/simulation/CFDSim.hh

@@ -62,7 +62,7 @@ namespace sim {
             #endif
             cfdSimulator.second->adSolve();
 
-            if (!cfdSimulator.second->hasConverged()) {
+            if (!cfdSimulator.second->hasAdConverged()) {
                 allConverge = false;
             }
 

src/simulation/MixingModels.hh

@@ -68,8 +68,6 @@ 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) ) {
@@ -84,21 +82,17 @@ void InstantaneousMixingModel<T>::propagateSpecies(arch::Network<T>* network, Si
         }
     }
 
-    std::cout << "[propagateSpecies] Initialize the node inflow" << std::endl;
     // Initial node outflow from mixtureInjections and CFD simulators, stored in mixtureOutflowAtNode
     initNodeOutflow(sim, tmpMixtures);
 
-    std::cout << "[propagateSpecies] Channel propagation" << std::endl;
     // Propagate the mixtures through the entire channel, without considering time steps    
     channelPropagation(network);
 
     bool inflowUpdated = true;
     while (inflowUpdated) {
-        std::cout << "[propagateSpecies] Update node outflow" << std::endl;
         // From node inflow, generate the node's outflow
         inflowUpdated = updateNodeOutflow(sim, tmpMixtures);
         // Propagate the mixtures through the entire channel
-        std::cout << "[propagateSpecies] Channel propagation" << std::endl;
         channelPropagation(network);
     }
 
@@ -111,7 +105,6 @@ void InstantaneousMixingModel<T>::propagateSpecies(arch::Network<T>* network, Si
 
 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();
@@ -120,7 +113,6 @@ 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()) {
@@ -216,25 +208,18 @@ bool InstantaneousMixingModel<T>::updateNodeOutflow(Simulation<T>* sim, std::vec
 
 template<typename T>
 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;
+                if (mixtureOutflowAtNode.count(nodeId)) {
+                    for (auto& [specieId, specieConcentration] : tmpMixtures[mixtureOutflowAtNode.at(nodeId)].getSpecieConcentrations()) {
+                        concentrations.at(nodeId).at(specieId) = specieConcentration;
+                    }
                 }
-                */
             }
         }
-        std::cout << std::endl;
         cfdSimulator->storeConcentrations(concentrations);
     }
 }

src/simulation/Simulation.hh

@@ -604,9 +604,7 @@ 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

@@ -158,6 +158,8 @@ class CFDSimulator {
 
     virtual void storeCfdResults (int iT) {}
 
+    virtual bool hasAdConverged() const { return false; }
+
     friend void coupleNsAdResults<T>(const std::unordered_map<int, std::unique_ptr<CFDSimulator<T>>>& cfdSimulators);
 
 };

src/simulation/simulators/olbMixing.h

@@ -50,6 +50,9 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
 
     std::unordered_map<int, Specie<T>*> species;
 
+    std::unordered_map<int, T> averageDensities;
+    std::unordered_map<int, bool> custConverges;
+
     std::unordered_map<int, std::shared_ptr<olb::SuperLattice<T, ADDESCRIPTOR>>> adLattices;      ///< The LBM lattice on the geometry.
     std::unordered_map<int, std::unique_ptr<olb::util::ValueTracer<T>>> adConverges;            ///< Value tracer to track convergence.
     std::unordered_map<int, std::shared_ptr<const olb::AdeUnitConverterFromResolutionAndRelaxationTime<T, ADDESCRIPTOR>>> adConverters;      ///< Object that stores conversion factors from phyical to lattice parameters.
@@ -174,9 +177,7 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
      * @brief Returns whether the module has converged or not.
      * @returns Boolean for module convergence.
     */
-    bool hasAdConverged(int key) const {
-        return adConverges.at(key)->hasConverged();
-    };
+    bool hasAdConverged() const override;
 
 };
 

src/simulation/simulators/olbMixing.hh

@@ -153,7 +153,15 @@ void lbmMixingSimulator<T>::writeVTK (int iT) {
 
         // write vtk to file system
         this->vtkFile = olb::singleton::directories().getVtkOutDir() + "data/" + olb::createFileName( this->name, iT ) + ".vtm";
-        this->converge->takeValue(this->getLattice().getStatistics().getAverageEnergy(), print);
+        this->converge->takeValue(this->getLattice().getStatistics().getAverageEnergy(), !print);
+        for (auto& [key, adConverge] : adConverges) {
+            //adConverge->takeValue(getAdLattice(key).getStatistics().getAverageRho(), print);
+            T newRho = getAdLattice(key).getStatistics().getAverageRho();
+            if (std::abs(averageDensities.at(key) - newRho) < 1e-5) {
+                custConverges.at(key) = true;
+            }
+            averageDensities.at(key) = newRho;
+        }
     }
     if (iT %1000 == 0) {
         #ifdef VERBOSE
@@ -257,7 +265,9 @@ template<typename T>
 void lbmMixingSimulator<T>::initAdConvergenceTracker () {
     // Initialize a convergence tracker for concentrations
     for (auto& [speciesId, specie] : species) {
-        this->adConverges.try_emplace(speciesId, std::make_unique<olb::util::ValueTracer<T>> (this->stepIter, this->epsilon));
+        this->adConverges.try_emplace(speciesId, std::make_unique<olb::util::ValueTracer<T>> (1000, 1e-1));
+        this->averageDensities.try_emplace(speciesId, T(0.0));
+        this->custConverges.try_emplace(speciesId, false);
     }
 }
 
@@ -379,4 +389,15 @@ std::unordered_map<int, std::unordered_map<int, T>> lbmMixingSimulator<T>::getCo
     return this->concentrations;
 }
 
+template<typename T>
+bool lbmMixingSimulator<T>::hasAdConverged() const {
+    bool c = true;
+    for (auto& [key, converge] : custConverges) {
+        if (!converge) {
+            c = false;
+        }
+    }
+    return c;
+};
+
 }   // namespace arch