diff --git a/src/architecture/Network.h b/src/architecture/Network.h
index 66063da..f3c3cfb 100755
--- a/src/architecture/Network.h
+++ b/src/architecture/Network.h
@@ -303,6 +303,13 @@ class Network {
      */
     bool isGround(int nodeId) const;
 
+    /**
+     * @brief Checks and returns if a node is a boundary node (connects to a CFD module).
+     * @param[in] nodeId Id of the node that should be checked.
+     * @return If the node with the specified id is a boundary node.
+     */
+    bool isBoundary(int nodeId) const;
+
     /**
      * @brief Checks and returns if an edge is a channel
     */
diff --git a/src/architecture/Network.hh b/src/architecture/Network.hh
index 16aa0a0..fc98117 100755
--- a/src/architecture/Network.hh
+++ b/src/architecture/Network.hh
@@ -427,6 +427,16 @@ bool Network<T>::isGround(int nodeId_) const {
     return nodes.at(nodeId_)->getGround();
 }
 
+template<typename T>
+bool Network<T>::isBoundary(int nodeId_) const {
+    if (modularReach.count(nodeId_)) {
+        if (modularReach.at(nodeId_)->getModuleType() == ModuleType::LBM) {
+            return true;
+        }
+    }
+    return false;
+}
+
 template<typename T>
 bool Network<T>::isChannel(int edgeId_) const {
     return channels.count(edgeId_);
diff --git a/src/simulation/MixingModels.h b/src/simulation/MixingModels.h
index 7fbd0c3..4c3ff70 100644
--- a/src/simulation/MixingModels.h
+++ b/src/simulation/MixingModels.h
@@ -34,6 +34,9 @@ class Simulation;
 template<typename T>
 class Specie;
 
+template<typename T>
+class lbmSimulator;
+
 // Structure to define the mixture inflow into a node
 template<typename T>
 struct MixtureInFlow {
@@ -251,13 +254,16 @@ template<typename T>
 class DiffusionMixingModel : public MixingModel<T> {
 
 private:
-    int resolution = 10;
+    int noFourierTerms = 100;
     std::set<int> mixingNodes;
     std::vector<std::vector<RadialPosition<T>>> concatenatedFlows;
     std::unordered_map<int, std::vector<FlowSection<T>>> outflowDistributions;
     std::unordered_map<int, int> filledEdges;                                   ///< Which edges are currently filled and what mixture is at the front <EdgeID, MixtureID>
+    std::unordered_map<int, std::unordered_map<int, std::vector<T>>>  concentrationFieldsOut; ///< Defines which concentration fields are defined at nodes at the interface between 1D into CFD <channelId, <specieId, concentrationField>>
     void generateInflows();
 
+    std::shared_ptr<lbmSimulator<T>> simulator; // Pointer to the lbm simulator
+
 public:
 
     DiffusionMixingModel();
@@ -290,13 +296,38 @@ class DiffusionMixingModel : public MixingModel<T> {
     
     void printTopology();
 
-    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalSolutionConstant(T channelLength, T channelWidth, int resolution, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters);
+    /**
+     * @brief Calculate the species concentration distribution across the channel width (at the end of the channel) for constant concentration flow sections. To get the analytical solution, a fourier series is employed. 
+     */
+    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalFunction(T channelLength, T channelWidth, int noFourierTerms, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters);
 
-    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalSolutionFunction(T channelLength, T channelWidth, int resolution, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters, std::function<T(T)> fConstant);
+    /**
+     * @brief Calculate the species concentration distribution across the channel width (at the end of the channel) for non-constant concentration inputs flow sections. To get the analytical solution, a fourier series is employed. 
+     */
+    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalFunction(T channelLength, T channelWidth, int noFourierTerms, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters, std::function<T(T)> fConstant);
 
-    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalSolutionTotal(T channelLength, T currChannelFlowRate, T channelWidth, int resolution, int speciesId, T pecletNr, 
+    /**
+     * @brief Calculate the species concentration distribution in a channel by calculating the get analytical function for constant and non-constant flow sections.
+     */
+    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalSolution(T channelLength, T currChannelFlowRate, T channelWidth, int noFourierTerms, int specieId, T pecletNr, 
         const std::vector<FlowSection<T>>& flowSections, std::unordered_map<int, std::unique_ptr<Mixture<T>>>& diffusiveMixtures);
 
+    /**
+     * @brief Calculate the species concentration across the channel width (at the end of the channel) for concentation fields flowing out of a CFD module. For this the constant concentration values at each lattice point of the CFD module are translated into linear segments.
+     * This way a more accurate connection, while adhering to the law of mass conservation is achieved. For the concentration distribution function at the end of a channel the analytical solution is solved by employing a fourier series. 
+     */
+    std::tuple<std::function<T(T)>,std::vector<T>,T> getAnalyticalSolutionHybrid(T channelLength, T currChannelFlowRate, T channelWidth, int resolution, T pecletNr, std::vector<T> concentrationField, T dx);
+
+    /**
+     * @brief Use piecewise linerar interpolation to translate the analytical solution of the species concentration across the channel width (at the end of the channel) into the concentration values in the lattice of the CFD module. This considers the conservation of mass. 
+     */
+    std::vector<T> defineConcentrationNodeFieldForCfdInput(int resolutionIntoCfd, int specieId, int channelId, T channelWidth, std::unordered_map<int, std::unique_ptr<Mixture<T>>>& Mixtures, int noFourierTerms);
+    
+    /**
+     * @brief Use a fifth order polynomial function to translate the analytical solution of the species concentration across the channel width (at the end of the channel) into the concentration values in the lattice of the CFD module.
+     */
+    std::tuple<std::function<T(T)>, std::vector<T>, T> getAnalyticalSolutionHybridInput(T channelLength, T channelWidth, int noFourierTerms, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters);
+
     void clean(arch::Network<T>* network);
 
     void printMixturesInNetwork();
@@ -309,6 +340,22 @@ class DiffusionMixingModel : public MixingModel<T> {
 
     bool isDiffusive() override { return 1; };
 
+    /**
+     * @brief Reset the number of Fourier terms used in the analytical solution. The higher the number the higher the precision, although this is negligible at some point.
+     * @param[in] newnoFourierTerms The new number of Fourier terms.
+     */
+    void resetnoFourierTerms(int newnoFourierTerms) {
+        noFourierTerms = newnoFourierTerms;
+    }
+
+    /**
+     * @brief Return the number of Fourier terms used in the analytical solution.
+     * @param[out] noFourierTerms The current number of Fourier terms.
+     */
+    int getnoFourierTerms() const {
+        return noFourierTerms;
+    }
+
 };
 
 }   // namespace sim
\ No newline at end of file
diff --git a/src/simulation/MixingModels.hh b/src/simulation/MixingModels.hh
index 04574e9..01fbb44 100644
--- a/src/simulation/MixingModels.hh
+++ b/src/simulation/MixingModels.hh
@@ -4,6 +4,7 @@
 #include <deque>
 #include <iostream>
 #include <cmath>
+#include <vector>
 
 #define M_PI 3.14159265358979323846
 
@@ -446,8 +447,8 @@ void DiffusionMixingModel<T>::generateInflows(T timeStep, arch::Network<T>* netw
                 std::unordered_map<int, std::tuple<std::function<T(T)>,std::vector<T>,T>> newDistributions;
                 for (auto& specieId : presentSpecies) {
                     T pecletNr = (std::abs(channel->getFlowRate()) / channel->getHeight()) / (sim->getSpecie(specieId))->getDiffusivity();
-                    std::tuple<std::function<T(T)>, std::vector<T>, T> analyticalResult = getAnalyticalSolutionTotal(
-                        channel->getLength(), std::abs(channel->getFlowRate()), channel->getWidth(), 100, specieId, 
+                    std::tuple<std::function<T(T)>, std::vector<T>, T> analyticalResult = getAnalyticalSolution(
+                        channel->getLength(), std::abs(channel->getFlowRate()), channel->getWidth(), noFourierTerms, specieId, 
                         pecletNr, outflowDistributions.at(channel->getId()), mixtures);
                     newDistributions.try_emplace(specieId, analyticalResult);
                 }
@@ -458,6 +459,42 @@ void DiffusionMixingModel<T>::generateInflows(T timeStep, arch::Network<T>* netw
             }
         }
     }
+    // Define the outflow of the CFD simulators for diffusive mixing
+    // These nodes are not defined as mixing nodes but they might still have a complex concentration distribution
+    for (auto& [key, cfdSimulator] : sim->getCFDSimulators()) {
+        std::unordered_map<int, std::unordered_map<int, std::vector<T>>> concentrationFieldsIn = cfdSimulator->getNodeConcentrationFields();
+            
+        for (auto& [nodeId, opening] : cfdSimulator->getOpenings()) {
+            // If the node is an inflow
+            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;
+                    for (const auto& [specieId, concentrationField] : concentrationFieldsIn.at(nodeId)) {
+                        std::cout << "Getting here 2" << std::endl;
+                        T dx = channel->getWidth() / concentrationField.size();
+                        T pecletNr = (std::abs(channel->getFlowRate()) / channel->getHeight()) / (sim->getSpecie(specieId))->getDiffusivity();
+                        std::tuple<std::function<T(T)>, std::vector<T>, T> analyticalResult = getAnalyticalSolutionHybrid(
+                            channel->getLength(), channel->getFlowRate(), channel->getWidth(), noFourierTerms, pecletNr, concentrationField, dx);
+                        newDistributions.try_emplace(specieId, analyticalResult);
+                    }
+                    //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)) {
+                    int resolutionIntoCfd = cfdSimulator->getResolution(nodeId);
+                    concentrationFieldsOut.try_emplace(nodeId, std::unordered_map<int, std::vector<T>>());
+                    for (const auto& [specieId, speciePtr] : sim->getSpecies()) {
+                        std::vector<T> concentrationFieldForCfdInput = defineConcentrationNodeFieldForCfdInput(resolutionIntoCfd, specieId, channelId, channel->getWidth(), mixtures, noFourierTerms);
+                        concentrationFieldsOut.at(nodeId).try_emplace(specieId, concentrationFieldForCfdInput);
+                    }                    
+                }
+            }
+        }
+        cfdSimulator->storeNodeConcentrationFields(concentrationFieldsOut);
+    }
 } 
 
 template<typename T>
@@ -476,7 +513,7 @@ void DiffusionMixingModel<T>::topologyAnalysis( arch::Network<T>* network, int n
     arch::Node<T>* node = network->getNode(nodeId).get();
     concatenatedFlows.clear();
     outflowDistributions.clear();
-    if (!node->getGround()){
+    if (!node->getGround() && !network->isBoundary(nodeId)){
 
         // 1. List all connecting channels with angle and inflow
         std::vector<RadialPosition<T>> channelOrder;
@@ -657,25 +694,48 @@ void DiffusionMixingModel<T>::topologyAnalysis( arch::Network<T>* network, int n
                 outflowDistributions.at(opposed.at(0).channelId).push_back(flow4);
             }
         } else {
-            throw std::invalid_argument("Too many channels at node " + std::to_string(nodeId) + " for DiffusionMixingModel.");
+            throw std::invalid_argument("Too many channels (" + std::to_string(concatenatedFlows.size()) + ") at node " + std::to_string(nodeId) + " for DiffusionMixingModel.");
         }
     }
 }
 
 template<typename T>
 void DiffusionMixingModel<T>::propagateSpecies(arch::Network<T>* network, Simulation<T>* sim) {
-    // TODO
+
+    // Initialization
+    T timeStep = 0.0;
+    bool inflowUpdated = true;
+
+    // As long as the remaining timestep does not return to 0.0, we conduct the following loop
+    while (inflowUpdated) {
+        // Propagate all mixtures and check if a mixture reaches a channel end
+        std::cout << "[propagateSpecies][updateNodeInflow]" << std::endl;
+        updateNodeInflow(timeStep, network);
+        // Generate a new inflow in case a mixture has reached channel end.
+        std::cout << "[propagateSpecies][generateInflows]" << std::endl;
+        generateInflows(timeStep, network, sim, sim->getMixtures());
+        // Clear the buffers for the diffusion mixing model
+        std::cout << "[propagateSpecies][clean]" << std::endl;
+        clean(network);
+        // Define and store the new minimal timestep for the next iteration
+        this->updateMinimalTimeStep(network);
+        timeStep = this->getMinimalTimeStep();
+        if (timeStep < 1e-12) {
+            inflowUpdated = false;
+        }
+    }
+
 }
 
 template<typename T>
-std::tuple<std::function<T(T)>, std::vector<T>, T> DiffusionMixingModel<T>::getAnalyticalSolutionConstant(T channelLength, T channelWidth, int resolution, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters) { 
+std::tuple<std::function<T(T)>, std::vector<T>, T> DiffusionMixingModel<T>::getAnalyticalFunction(T channelLength, T channelWidth, int noFourierTerms, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters) { 
     T a_0 = 0.0;
     std::vector<T> segmentedResult;
 
     for (const auto& parameter : parameters) {
-        for (int n = 1; n < resolution; n++) {
+        for (int n = 1; n < noFourierTerms; n++) {
             T a_n = (2/(n * M_PI))  * (parameter.concentrationAtChannelEnd) * (std::sin(n * M_PI * parameter.endWidth) - std::sin(n * M_PI * parameter.startWidth)); 
-                segmentedResult.push_back(a_n * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth)));
+            segmentedResult.push_back(a_n * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth)));
         }
     }
 
@@ -683,11 +743,11 @@ std::tuple<std::function<T(T)>, std::vector<T>, T> DiffusionMixingModel<T>::getA
             a_0 += 2 * parameter.concentrationAtChannelEnd  * (parameter.endWidth - parameter.startWidth);
         }
 
-    auto f = [a_0, channelLength, channelWidth, resolution, pecletNr, parameters](T w) { // This returns C(w, l_1)
+    auto f = [a_0, channelLength, channelWidth, noFourierTerms, pecletNr, parameters](T w) { // This returns C(w, l_1)
         T f_sum = 0.0;
 
         for (const auto& parameter : parameters) {
-            for (int n = 1; n < resolution; n++) {
+            for (int n = 1; n < noFourierTerms; n++) {
                 T a_n = (2/(n * M_PI))  * (parameter.concentrationAtChannelEnd) * (std::sin(n * M_PI * parameter.endWidth) - std::sin(n * M_PI * parameter.startWidth));
                 f_sum += a_n * std::cos(n * M_PI * (w)) * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth));
             }
@@ -698,7 +758,7 @@ std::tuple<std::function<T(T)>, std::vector<T>, T> DiffusionMixingModel<T>::getA
 }
 
 template<typename T> 
-std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAnalyticalSolutionFunction(T channelLength, T channelWidth, int resolution, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters, std::function<T(T)> fConstant) { 
+std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAnalyticalFunction(T channelLength, T channelWidth, int noFourierTerms, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters, std::function<T(T)> fConstant) { 
     // From Channel Start to Channel End for complex input
     std::vector<T> segmentedResult;
     T a_0 = 0.0;
@@ -708,13 +768,13 @@ std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAna
     for (const auto& parameter : parameters) {
         T a_0_old = parameter.a_0_old;
         T scaleFactor = parameter.scaleFactor;
-        for (int n = 1; n < resolution; n++) {
+        for (int n = 1; n < noFourierTerms; n++) {
             a_n = a_0_old / (M_PI * n) * (std::sin(n * M_PI * parameter.endWidth) - std::sin(n * M_PI * parameter.startWidth));
             
             for (long unsigned int i = 0; i < parameter.segmentedResult.size(); i++) {
                 T translateFactor = parameter.translateFactor;
 
-                int oldN = (i % (resolution - 1)) + 1;
+                int oldN = (i % (noFourierTerms - 1)) + 1;
 
                 if (abs(oldN/scaleFactor - n) < 1e-8) { 
                     a_n += 2 * ((0.5 * parameter.endWidth - 0.5 * parameter.startWidth) * std::cos(oldN * M_PI * translateFactor) 
@@ -743,14 +803,14 @@ std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAna
         for (long unsigned int i = 0; i < parameter.segmentedResult.size(); i++) { 
             T scaleFactor = parameter.scaleFactor;
             T translateFactor = parameter.translateFactor;
-            int oldN = (i % (resolution - 1)) + 1;
+            int oldN = (i % (noFourierTerms - 1)) + 1;
             a_0 += 2 * parameter.segmentedResult[i] * scaleFactor/(oldN * M_PI) * 
                 (std::sin(oldN * M_PI * parameter.endWidth / scaleFactor + oldN * M_PI * translateFactor) 
                 - std::sin(oldN * M_PI * parameter.startWidth / scaleFactor + oldN * M_PI * translateFactor));
         }
     }     
 
-    auto f = [a_0, channelLength, channelWidth, resolution, pecletNr, parameters, fConstant](T w) { // This returns C(w, l_1)
+    auto f = [a_0, channelLength, channelWidth, noFourierTerms, pecletNr, parameters, fConstant](T w) { // This returns C(w, l_1)
         T f_sum = 0.0;
         T a_n = 0.0;
 
@@ -758,12 +818,12 @@ std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAna
         for (const auto& parameter : parameters) {
             T a_0_old = parameter.a_0_old;
             T scaleFactor = parameter.scaleFactor;
-            for (int n = 1; n < resolution; n++) {
+            for (int n = 1; n < noFourierTerms; n++) {
                 a_n = a_0_old / (M_PI * n) * (std::sin(n * M_PI * parameter.endWidth) - std::sin(n * M_PI * parameter.startWidth));
                 
                 for (long unsigned int i = 0; i < parameter.segmentedResult.size(); i++) {
                     T translateFactor = parameter.translateFactor;
-                    int oldN = (i % (resolution - 1)) + 1;
+                    int oldN = (i % (noFourierTerms - 1)) + 1;
 
                     if (abs(oldN/scaleFactor - n) < 1e-12) { 
                         a_n += 2 * ((0.5 * parameter.endWidth - 0.5 * parameter.startWidth) * std::cos(oldN * M_PI * translateFactor) 
@@ -792,8 +852,8 @@ std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAna
 
 // From Channel Start to Channel End for complex input
 template<typename T>
-std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAnalyticalSolutionTotal(
-    T channelLength, T currChannelFlowRate, T channelWidth, int resolution, int speciesId, 
+std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAnalyticalSolution(
+    T channelLength, T currChannelFlowRate, T channelWidth, int noFourierTerms, int specieId, 
     T pecletNr, const std::vector<FlowSection<T>>& flowSections, std::unordered_map<int, std::unique_ptr<Mixture<T>>>& Mixtures) { 
     
     T prevEndWidth = 0.0;
@@ -819,18 +879,18 @@ std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAna
         } else {
             T mixtureId = this->filledEdges.at(flowSection.channelId); // get the diffusive mixture at a specific channelId
             DiffusiveMixture<T>* diffusiveMixture = dynamic_cast<DiffusiveMixture<T>*>(Mixtures.at(mixtureId).get()); // Assuming diffusiveMixtures is passed
-            if (diffusiveMixture->getSpecieConcentrations().find(speciesId) == diffusiveMixture->getSpecieConcentrations().end()) { // the species is not in the mixture
+            if (diffusiveMixture->getSpecieConcentrations().find(specieId) == diffusiveMixture->getSpecieConcentrations().end()) { // the species is not in the mixture
                 T concentration = 0.0;
                 std::function<T(T)> zeroFunction = [](T) -> T { return 0.0; };
                 std::vector<T> zeroSegmentedResult = {0};
                 constantFlowSections.push_back({startWidth, endWidth, 1.0, translateFactor, concentration, zeroFunction, zeroSegmentedResult, T(0.0)});
             } else if (diffusiveMixture->getIsConstant()) { 
-                T concentration = diffusiveMixture->getConcentrationOfSpecie(speciesId); 
+                T concentration = diffusiveMixture->getConcentrationOfSpecie(specieId); 
                 std::function<T(T)> zeroFunction = [](T) -> T { return 0.0; };
                 std::vector<T> zeroSegmentedResult = {0};
                 constantFlowSections.push_back({startWidth, endWidth, 1.0, translateFactor, concentration, zeroFunction, zeroSegmentedResult, T(0.0)}); 
             } else {
-                auto specieDistribution = diffusiveMixture->getSpecieDistributions().at(speciesId);
+                auto specieDistribution = diffusiveMixture->getSpecieDistributions().at(specieId);
                 // std::function<T(T)> concentrationFunction = specieDistribution.first;
                 // std::vector<T> segmentedResults = specieDistribution.second;
                 std::function<T(T)> concentrationFunction = std::get<0>(specieDistribution);
@@ -843,8 +903,8 @@ std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAna
         }
     }
 
-    auto [fConstant, segmentedResultConstant, a_0_Constant] = getAnalyticalSolutionConstant(channelLength, channelWidth, resolution, pecletNr, constantFlowSections);
-    auto [fFunction, segmentedResultFunction, a_0_Function] = getAnalyticalSolutionFunction(channelLength, channelWidth, resolution, pecletNr, functionFlowSections, fConstant);
+    auto [fConstant, segmentedResultConstant, a_0_Constant] = getAnalyticalFunction(channelLength, channelWidth, noFourierTerms, pecletNr, constantFlowSections);
+    auto [fFunction, segmentedResultFunction, a_0_Function] = getAnalyticalFunction(channelLength, channelWidth, noFourierTerms, pecletNr, functionFlowSections, fConstant);
 
     segmentedResultFunction.insert(segmentedResultFunction.end(), segmentedResultConstant.begin(), segmentedResultConstant.end());
 
@@ -854,7 +914,165 @@ std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAna
 
 }
 
+template<typename T>
+std::tuple<std::function<T(T)>,std::vector<T>,T> DiffusionMixingModel<T>::getAnalyticalSolutionHybrid(
+    T channelLength, T currChannelFlowRate, T channelWidth, int noFourierTerms, T pecletNr, std::vector<T> concentrationField, T dx) { 
+
+    T distance2Wall = dx / 2.0;
+    std::vector<T> segmentedResult;
+    std::vector<T> a_n_values(noFourierTerms, 0.0);
+    T a_0 = 0.0;
+
+    for (size_t i = 0; i < concentrationField.size() - 1; ++i) {
+        T width1 = distance2Wall + i * dx;
+        T concentration1 = concentrationField[i];
+        T width2 = distance2Wall + (i + 1) * dx;
+        T concentration2 = concentrationField[i + 1];
+
+        T linearSegmentFactor = (concentration2 - ((concentration1 * width2) / width1)) / (1 - (width2 / width1));
+        T linearSegmentSlope = (concentration1 - linearSegmentFactor) / width1;
+
+        a_0 += 2 * linearSegmentFactor * (width2 - width1) + linearSegmentSlope * (pow(width2, 2) - pow(width1, 2));
+
+        for (int n = 1; n < noFourierTerms; n++) {
+            T a_n = 2 / (n * M_PI) * std::sin(n * M_PI * width2) * (linearSegmentFactor + linearSegmentSlope * width2) + 2 / pow((n * M_PI), 2) * linearSegmentSlope * std::cos(n * M_PI * width2)
+                - 2 / (n * M_PI) * std::sin(n * M_PI * width1) * (linearSegmentFactor + linearSegmentSlope * width1) - 2 / pow((n * M_PI), 2) * linearSegmentSlope * std::cos(n * M_PI * width1);
+            segmentedResult.push_back(a_n * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth)));
+            a_n_values[n] += a_n;
+        }
+    }
+    // add the connection to the sides
+    // start segment
+    T concentrationStart = concentrationField[0];
+    T startWidth = 0.0;
+    a_0 += 2 * (concentrationStart * (distance2Wall - startWidth));
+    for (int n = 1; n < noFourierTerms; n++) {
+        T a_n_start = (2/(n * M_PI)) * concentrationStart * (std::sin(n * M_PI * distance2Wall) - std::sin(n * M_PI * startWidth));
+        a_n_values[n] += a_n_start;
+        segmentedResult.insert(segmentedResult.begin(), a_n_start * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth)));
+    }
+    // end segment
+    T concentrationEnd = concentrationField[concentrationField.size() - 1];
+    T endWidth = dx * concentrationField.size();
+    a_0 += 2 * (concentrationEnd * (endWidth - (endWidth - distance2Wall)));
+    for (int n = 1; n < noFourierTerms; n++) {
+        T a_n_end = (2/(n * M_PI)) * concentrationEnd * (std::sin(n * M_PI * endWidth) - std::sin(n * M_PI * (endWidth - distance2Wall)));
+        a_n_values[n] += a_n_end;
+        segmentedResult.push_back(a_n_end * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth)));
+    }
+
+    auto functionHybrid = [a_0, a_n_values, channelLength, channelWidth, noFourierTerms, pecletNr](T w) { // This returns C(w, l_1)
+        T f_sum = 0.0;
+        for (int n = 1; n < noFourierTerms; n++) {
+            f_sum += a_n_values[n] * std::cos(n * M_PI * w) * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength / channelWidth));
+        }
+        return 0.5 * a_0 + f_sum;
+    };
 
+    return {functionHybrid, segmentedResult, a_0};
+}
+
+template<typename T>
+std::vector<T> DiffusionMixingModel<T>::defineConcentrationNodeFieldForCfdInput(int resolutionIntoCfd, int specieId, int channelId, T channelWidth, std::unordered_map<int, std::unique_ptr<Mixture<T>>>& Mixtures, int noFourierTerms) {
+    // To adhere to the conservation of mass, the sum of concentration passed must be equal to that represented by the concentration function across the width.
+    // For this piecewise linerar interpolation is employed.
+    std::vector<T> concentrationField;
+    int mixtureId = 0;
+    if (this->filledEdges.count(channelId)) {
+        mixtureId = this->filledEdges.at(channelId); // get the diffusive mixture at a specific channelId
+    } else {
+        concentrationField = std::vector<T> (resolutionIntoCfd, 0.0);
+        return concentrationField;
+    }
+
+    DiffusiveMixture<T>* diffusiveMixture = dynamic_cast<DiffusiveMixture<T>*>(Mixtures.at(mixtureId).get()); // Assuming diffusiveMixtures is passed
+    if (diffusiveMixture->getSpecieConcentrations().find(specieId) == diffusiveMixture->getSpecieConcentrations().end()) { // the species is not in the mixture
+        // here we could distinguish between constant and non-constant concentration
+        concentrationField = std::vector<T>(resolutionIntoCfd, 0.0);
+    } else if (diffusiveMixture->getIsConstant()) { 
+        T concentration = diffusiveMixture->getConcentrationOfSpecie(specieId);
+        concentrationField = std::vector<T>(resolutionIntoCfd, concentration);
+    } else {
+        auto specieDistribution = diffusiveMixture->getSpecieDistributions().at(specieId);
+        std::vector<T> segmentedResult = std::get<1>(specieDistribution);
+        T a_0_old = std::get<2>(specieDistribution);
+        for (int i = 0; i < resolutionIntoCfd - 1; i++) {
+            T segmentWidth = channelWidth / resolutionIntoCfd;
+            T start = i * segmentWidth;
+            T end = (i + 1) * segmentWidth; 
+            // now we to integrate the function over the width of each dx
+            T concentration = 0.5 * M_PI * a_0_old * (end - start); 
+            for (long unsigned int n = 0; n < segmentedResult.size(); n++) { 
+                int oldN = (n % (noFourierTerms - 1)) + 1;
+                concentration += segmentedResult[n] / oldN * (std::sin(oldN * M_PI * end) - std::sin(oldN * M_PI * start));
+            }
+            concentrationField.push_back(concentration);
+        }
+    }
+    return concentrationField;
+}
+
+template<typename T>
+std::tuple<std::function<T(T)>, std::vector<T>, T> DiffusionMixingModel<T>::getAnalyticalSolutionHybridInput(
+    T channelLength, T channelWidth, int noFourierTerms, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters) { 
+    T a_0 = 0.0;
+    std::vector<T> segmentedResult;
+
+    // TODO - get the parameters a to f from the polynomial fit of the hybrid concentration output
+    for (const auto& parameter : parameters) {
+        for (int n = 1; n < noFourierTerms; n++) {
+            T a_n = 1 / pow(n * M_PI, 6) * (n * M_PI * std::sin(n * M_PI * parameter.endWidth) * (parameter.a * pow(n * M_PI, 4) * pow(parameter.endWidth, 5) - 20 * parameter.a * pow(n * M_PI, 2) * pow(parameter.endWidth, 3) 
+                + 120 * parameter.a * parameter.endWidth + parameter.b * (pow(n * M_PI, 4) * pow(parameter.endWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.endWidth, 2) + 24)
+                + parameter.c * pow(n * M_PI, 4) * pow(parameter.endWidth, 3) - 6 * parameter.c * pow(n * M_PI, 2) * parameter.endWidth + parameter.d * pow(n * M_PI, 4) * pow(parameter.endWidth, 2)
+                - 2 * parameter.d * pow(n * M_PI, 2) + parameter.e * pow(n * M_PI, 4) * parameter.endWidth + parameter.f * pow(n * M_PI, 4))
+                + std::cos(n * M_PI * parameter.endWidth) * (5 * parameter.a * (pow(n * M_PI, 4) * pow(parameter.endWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.endWidth, 2) + 24)
+                + pow(n * M_PI, 2) * (4 * parameter.b * pow(n * M_PI, 2) * pow(parameter.endWidth, 3) - 24 * parameter.b * parameter.endWidth + 3 * parameter.c * pow(n * M_PI, 2) * pow(parameter.endWidth, 2) - 6 * parameter.c 
+                + 2 * parameter.d * pow(n * M_PI, 2) * parameter.endWidth + parameter.e * pow(n * M_PI, 2))))
+                -
+                (1 / pow(n * M_PI, 6) * (n * M_PI * std::sin(n * M_PI * parameter.startWidth) * (parameter.a * pow(n * M_PI, 4) * pow(parameter.startWidth, 5) - 20 * parameter.a * pow(n * M_PI, 2) * pow(parameter.startWidth, 3) 
+                + 120 * parameter.a * parameter.startWidth + parameter.b * (pow(n * M_PI, 4) * pow(parameter.startWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.startWidth, 2) + 24)
+                + parameter.c * pow(n * M_PI, 4) * pow(parameter.startWidth, 3) - 6 * parameter.c * pow(n * M_PI, 2) * parameter.startWidth + parameter.d * pow(n * M_PI, 4) * pow(parameter.startWidth, 2)
+                - 2 * parameter.d * pow(n * M_PI, 2) + parameter.e * pow(n * M_PI, 4) * parameter.startWidth + parameter.f * pow(n * M_PI, 4))
+                + std::cos(n * M_PI * parameter.startWidth) * (5 * parameter.a * (pow(n * M_PI, 4) * pow(parameter.startWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.startWidth, 2) + 24)
+                + pow(n * M_PI, 2) * (4 * parameter.b * pow(n * M_PI, 2) * pow(parameter.startWidth, 3) - 24 * parameter.b * parameter.endWidth + 3 * parameter.c * pow(n * M_PI, 2) * pow(parameter.startWidth, 2) - 6 * parameter.c 
+                + 2 * parameter.d * pow(n * M_PI, 2) * parameter.endWidth + parameter.e * pow(n * M_PI, 2))))); 
+            segmentedResult.push_back(a_n * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth)));
+        }
+    }
+
+    for (const auto& parameter : parameters) { // iterates through all channels that flow into the current channel
+            a_0 += parameter.a / 6 * pow(parameter.endWidth, 6) + parameter.b / 5 * pow(parameter.endWidth, 5) + parameter.c / 4 * pow(parameter.endWidth, 4) + parameter.d / 3 * pow(parameter.endWidth, 3) + parameter.e / 2 * pow(parameter.endWidth, 2) + parameter.f * parameter.endWidth
+                - (parameter.a / 6 * pow(parameter.startWidth, 6) + parameter.b / 5 * pow(parameter.startWidth, 5) + parameter.c / 4 * pow(parameter.startWidth, 4) + parameter.d / 3 * pow(parameter.startWidth, 3) + parameter.e / 2 * pow(parameter.startWidth, 2) + parameter.f * parameter.startWidth);
+        }
+
+    auto f = [a_0, channelLength, channelWidth, noFourierTerms, pecletNr, parameters](T w) { // This returns C(w, l_1)
+        T f_sum = 0.0;
+
+        for (const auto& parameter : parameters) {
+            for (int n = 1; n < noFourierTerms; n++) {
+                T a_n = 1 / pow(n * M_PI, 6) * (n * M_PI * std::sin(n * M_PI * parameter.endWidth) * (parameter.a * pow(n * M_PI, 4) * pow(parameter.endWidth, 5) - 20 * parameter.a * pow(n * M_PI, 2) * pow(parameter.endWidth, 3) 
+                + 120 * parameter.a * parameter.endWidth + parameter.b * (pow(n * M_PI, 4) * pow(parameter.endWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.endWidth, 2) + 24)
+                + parameter.c * pow(n * M_PI, 4) * pow(parameter.endWidth, 3) - 6 * parameter.c * pow(n * M_PI, 2) * parameter.endWidth + parameter.d * pow(n * M_PI, 4) * pow(parameter.endWidth, 2)
+                - 2 * parameter.d * pow(n * M_PI, 2) + parameter.e * pow(n * M_PI, 4) * parameter.endWidth + parameter.f * pow(n * M_PI, 4))
+                + std::cos(n * M_PI * parameter.endWidth) * (5 * parameter.a * (pow(n * M_PI, 4) * pow(parameter.endWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.endWidth, 2) + 24)
+                + pow(n * M_PI, 2) * (4 * parameter.b * pow(n * M_PI, 2) * pow(parameter.endWidth, 3) - 24 * parameter.b * parameter.endWidth + 3 * parameter.c * pow(n * M_PI, 2) * pow(parameter.endWidth, 2) - 6 * parameter.c 
+                + 2 * parameter.d * pow(n * M_PI, 2) * parameter.endWidth + parameter.e * pow(n * M_PI, 2))))
+                -
+                (1 / pow(n * M_PI, 6) * (n * M_PI * std::sin(n * M_PI * parameter.startWidth) * (parameter.a * pow(n * M_PI, 4) * pow(parameter.startWidth, 5) - 20 * parameter.a * pow(n * M_PI, 2) * pow(parameter.startWidth, 3) 
+                + 120 * parameter.a * parameter.startWidth + parameter.b * (pow(n * M_PI, 4) * pow(parameter.startWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.startWidth, 2) + 24)
+                + parameter.c * pow(n * M_PI, 4) * pow(parameter.startWidth, 3) - 6 * parameter.c * pow(n * M_PI, 2) * parameter.startWidth + parameter.d * pow(n * M_PI, 4) * pow(parameter.startWidth, 2)
+                - 2 * parameter.d * pow(n * M_PI, 2) + parameter.e * pow(n * M_PI, 4) * parameter.startWidth + parameter.f * pow(n * M_PI, 4))
+                + std::cos(n * M_PI * parameter.startWidth) * (5 * parameter.a * (pow(n * M_PI, 4) * pow(parameter.startWidth, 4) - 12 * pow(n * M_PI, 2) * pow(parameter.startWidth, 2) + 24)
+                + pow(n * M_PI, 2) * (4 * parameter.b * pow(n * M_PI, 2) * pow(parameter.startWidth, 3) - 24 * parameter.b * parameter.startWidth + 3 * parameter.c * pow(n * M_PI, 2) * pow(parameter.startWidth, 2) - 6 * parameter.c 
+                + 2 * parameter.d * pow(n * M_PI, 2) * parameter.startWidth + parameter.e * pow(n * M_PI, 2)))));
+
+                f_sum += a_n * std::cos(n * M_PI * (w)) * std::exp(-pow(n, 2) * pow(M_PI, 2) * (1 / pecletNr) * (channelLength/channelWidth));
+            }
+        }
+        return 0.5 * a_0 + f_sum;
+    };
+    return {f, segmentedResult, a_0};
+}
 
 template<typename T>
 void DiffusionMixingModel<T>::printTopology() {
@@ -890,6 +1108,7 @@ void DiffusionMixingModel<T>::clean(arch::Network<T>* network) {
             }
         }
     }
+    concentrationFieldsOut.clear();
 }
 
 template<typename T>
diff --git a/src/simulation/simulators/cfdSimulator.h b/src/simulation/simulators/cfdSimulator.h
index fadb0ce..47ca506 100755
--- a/src/simulation/simulators/cfdSimulator.h
+++ b/src/simulation/simulators/cfdSimulator.h
@@ -203,6 +203,17 @@ class CFDSimulator {
      */
     virtual std::unordered_map<int, T> getFlowRates() const = 0;
 
+    /**
+     * @brief Get the lattice resolution at the boundary nodes.
+     * @param[in] nodeId of the boundary node.
+     * @returns Resolution.
+     */
+    virtual int getResolution(int nodeId) const 
+    { 
+        throw std::runtime_error("The function getResolution is undefined for this CFD simulator.");
+        return 0;
+    }
+
     /**
      * @brief Store the abstract concentrations at the nodes on the module boundary in the simulator.
      * @param[in] concentrations Map of concentrations and node ids.
@@ -212,16 +223,35 @@ class CFDSimulator {
         throw std::runtime_error("The function storeConcentrations is undefined for this CFD simulator.");
     }
 
+    /**
+     * @brief Store the abstract concentration fields at the nodes on the module boundary in the simulator.
+     * @param[in] concentrationFieldsOut Map consisting of species id and a vector of concentrations at each grid point and node ids.
+     */
+    virtual void storeNodeConcentrationFields(std::unordered_map<int, std::unordered_map<int, std::vector<T>>> concentrationFieldsOut) 
+    {
+        throw std::runtime_error("The function storeNodeConcentrationFields is undefined for this CFD simulator.");
+    }
+
     /**
      * @brief Get the concentrations at the boundary nodes.
      * @returns Concentrations
      */
     virtual std::unordered_map<int, std::unordered_map<int, T>> getConcentrations() const 
     { 
-        throw std::runtime_error("The function storeConcentrations is undefined for this CFD simulator.");
+        throw std::runtime_error("The function getConcentrations is undefined for this CFD simulator.");
         return std::unordered_map<int, std::unordered_map<int, T>>(); 
     }
 
+    /**
+     * @brief Get the concentrations at the boundary nodes.
+     * @returns Concentrations
+     */
+    virtual std::unordered_map<int, std::unordered_map<int, std::vector<T>>> getNodeConcentrationFields() const 
+    { 
+        throw std::runtime_error("The function getNodeConcentrationFields is undefined for this CFD simulator.");
+        return std::unordered_map<int, std::unordered_map<int, std::vector<T>>>(); 
+    }
+
     /**
      * @brief Set the boundary values on the lattice at the module nodes.
      * @param[in] iT Iteration step.
diff --git a/src/simulation/simulators/olbContinuous.h b/src/simulation/simulators/olbContinuous.h
index 03d355e..9ff6572 100644
--- a/src/simulation/simulators/olbContinuous.h
+++ b/src/simulation/simulators/olbContinuous.h
@@ -74,6 +74,7 @@ using BounceBack = olb::BounceBack<T,DESCRIPTOR>;
     std::shared_ptr<const olb::UnitConverterFromResolutionAndRelaxationTime<T, DESCRIPTOR>> converter;      ///< Object that stores conversion factors from phyical to lattice parameters.
     std::unordered_map<int, std::shared_ptr<olb::SuperPlaneIntegralFluxVelocity2D<T>>> fluxes;              ///< Map of fluxes at module nodes. 
     std::unordered_map<int, std::shared_ptr<olb::SuperPlaneIntegralFluxPressure2D<T>>> meanPressures;       ///< Map of mean pressure values at module nodes.
+    std::unordered_map<int, int> resolutions;                       ///< Vector of Resolutions at module nodes <nodeId, resolution>
 
     auto& getConverter() {
         return *converter;
@@ -262,6 +263,15 @@ using BounceBack = olb::BounceBack<T,DESCRIPTOR>;
         return resolution; 
     };
 
+    /**
+     * @brief Get the resolution at a boundary node.
+     * @param[in] nodeId id of the node.
+     * @returns resolution
+     */
+    int getResolution(int nodeId) const override {
+        return resolutions.at(nodeId);
+    }
+
     /**
      * @brief Get the convergence criterion.
      * @returns epsilon.
diff --git a/src/simulation/simulators/olbContinuous.hh b/src/simulation/simulators/olbContinuous.hh
index 46b0370..5d49fcb 100644
--- a/src/simulation/simulators/olbContinuous.hh
+++ b/src/simulation/simulators/olbContinuous.hh
@@ -12,6 +12,9 @@ lbmSimulator<T>::lbmSimulator (
         epsilon(epsilon_), relaxationTime(relaxationTime_)
 { 
     this->cfdModule->setModuleTypeLbm();
+    for (auto& [key, opening] : openings_) {
+        resolutions.try_emplace(key, resolution_);
+    }
 } 
 
 template<typename T>
diff --git a/src/simulation/simulators/olbMixing.h b/src/simulation/simulators/olbMixing.h
index 368c4e1..82d9400 100644
--- a/src/simulation/simulators/olbMixing.h
+++ b/src/simulation/simulators/olbMixing.h
@@ -45,6 +45,7 @@ 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, Specie<T>*> species;
 
@@ -187,12 +188,24 @@ using NoADDynamics = olb::NoDynamics<T,ADDESCRIPTOR>;
      */
     void storeConcentrations(std::unordered_map<int, std::unordered_map<int, T>> concentrations) override;
 
+    /**
+     * @brief Store the abstract concentration fields at the nodes on the module boundary in the simulator.
+     * @param[in] concentrations Map of concentration vectors and node ids for each speciesId.
+     */
+    void storeNodeConcentrationFields(std::unordered_map<int, std::unordered_map<int, std::vector<T>>> nodeConcentrationFields) override;
+
     /**
      * @brief Get the concentrations at the boundary nodes.
      * @returns Concentrations
      */
     std::unordered_map<int, std::unordered_map<int, T>> getConcentrations() const override;
 
+    /**
+     * @brief Get the concentration field at the boundary nodes, i.e., the concentration at each grid end point in the lattice.
+     * @returns Concentration vectors
+     */
+    std::unordered_map<int, std::unordered_map<int, std::vector<T>>> getNodeConcentrationFields() const override;
+
     /**
      * @brief Returns whether the module has converged or not.
      * @returns Boolean for module convergence.
diff --git a/src/simulation/simulators/olbMixing.hh b/src/simulation/simulators/olbMixing.hh
index 188ddcb..ccecf00 100644
--- a/src/simulation/simulators/olbMixing.hh
+++ b/src/simulation/simulators/olbMixing.hh
@@ -246,10 +246,13 @@ void lbmMixingSimulator<T>::initValueContainers () {
         this->pressures.try_emplace(key, (T) 0.0);
         this->flowRates.try_emplace(key, (T) 0.0);
         std::unordered_map<int, T> tmpConcentrations;
+        std::unordered_map<int, std::vector<T>> tmpFieldConcentrations;
         for (auto& [speciesId, speciesPtr] : species) {
             tmpConcentrations.try_emplace(speciesId, 0.0);
+            tmpFieldConcentrations.try_emplace(speciesId, std::vector<T>(this->getResolution(key), 0.0));
         }
         this->concentrations.try_emplace(key, tmpConcentrations);
+        this->nodeConcentrationFields.try_emplace(key, tmpFieldConcentrations);
     }
 }
 
@@ -395,11 +398,21 @@ void lbmMixingSimulator<T>::storeConcentrations(std::unordered_map<int, std::uno
     this->concentrations = concentrations_;
 }
 
+template<typename T>
+void lbmMixingSimulator<T>::storeNodeConcentrationFields(std::unordered_map<int, std::unordered_map<int, std::vector<T>>> concentrationFieldsOut_) {
+    this->nodeConcentrationFields = concentrationFieldsOut_;
+}
+
 template<typename T>
 std::unordered_map<int, std::unordered_map<int, T>> lbmMixingSimulator<T>::getConcentrations() const {
     return this->concentrations;
 }
 
+template<typename T>
+std::unordered_map<int, std::unordered_map<int, std::vector<T>>> lbmMixingSimulator<T>::getNodeConcentrationFields() const {
+    return this->nodeConcentrationFields;
+}
+
 template<typename T>
 bool lbmMixingSimulator<T>::hasAdConverged() const {
     bool c = true;