Diffusive Mixing (1D <--> CFD)

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 {
@@ -59,6 +62,12 @@ struct FlowSection {
     T width;
 };
 
+template<typename T> // TODO this is redundant -> remove
+struct HybridFlowSection : public FlowSection<T> {
+    T hybridResolution;
+    std::vector<T> concentrations; // Alternatively we could see this as one flow section that contains all concnetrations across the connection (maybe as a vector)
+};
+
 template<typename T>
 struct FlowSectionInput {
     T startWidth;
@@ -251,17 +260,28 @@ template<typename T>
 class DiffusionMixingModel : public MixingModel<T> {
 
 private:
-    int resolution = 10;
+    int numFourierTerms = 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();
 
+    lbmSimulator<T>* simulator; // Pointer to the lbmSimulator
+    // std::shared_ptr<lbmSimulator<T>> simulator; 
+
 public:
 
     DiffusionMixingModel();
 
+    // DiffusionMixingModel(lbmSimulator<T>* simulator);
+
+    // /**
+    //  * @brief Set the simulator for the hybrid mixing simulation.
+    // */
+    // void setSimulator(lbmSimulator<T>* simulator_);
+
     /**
      * @brief Create and/or propagate mixtures into channels downstream.
      * @param[in] timeStep the current timestep size.
@@ -290,13 +310,19 @@ 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);
+    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalFunction(T channelLength, T channelWidth, int numFourierTerms, 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);
+    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalFunction(T channelLength, T channelWidth, int numFourierTerms, 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, 
+    std::tuple<std::function<T(T)>,std::vector<T>, T> getAnalyticalSolution(T channelLength, T currChannelFlowRate, T channelWidth, int numFourierTerms, int specieId, T pecletNr, 
         const std::vector<FlowSection<T>>& flowSections, std::unordered_map<int, std::unique_ptr<Mixture<T>>>& diffusiveMixtures);
 
+    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);
+
+    std::vector<T> defineConcentrationNodeFieldForCfdInput(int resolutionIntoCfd, int specieId, int channelId, T channelWidth, std::unordered_map<int, std::unique_ptr<Mixture<T>>>& Mixtures, int numFourierTerms);
+
+    std::tuple<std::function<T(T)>, std::vector<T>, T> getAnalyticalSolutionHybridInput(T channelLength, T channelWidth, int numFourierTerms, T pecletNr, const std::vector<FlowSectionInput<T>>& parameters);
+
     void clean(arch::Network<T>* network);
 
     void printMixturesInNetwork();
@@ -309,6 +335,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] newNumFourierTerms The new number of Fourier terms.
+     */
+    void resetNumFourierTerms(int newNumFourierTerms) {
+        numFourierTerms = newNumFourierTerms;
+    }
+
+    /**
+     * @brief Return the number of Fourier terms used in the analytical solution.
+     * @param[out] numFourierTerms The current number of Fourier terms.
+     */
+    int getNumFourierTerms() const {
+        return numFourierTerms;
+    }
+
 };
 
 }   // namespace sim