🎨 Add PPM image output

cda-tum · Sep 6, 2024 · 6bef548 · 6bef548
1 parent ec84c7c
commit 6bef548
Show file tree

Hide file tree

Showing 5 changed files with 144 additions and 2 deletions.
diff --git a/src/simulation/Simulation.h b/src/simulation/Simulation.h
@@ -147,6 +147,7 @@ class Simulation {
     int continuousPhase = 0;                                                            ///< Fluid of the continuous phase.
     int iteration = 0;
     int maxIterations = 1e5;
+    int globalSteps = 0;
     T maximalAdaptiveTimeStep = 0;                                                      ///< Maximal adaptive time step that is applied when droplets change the channel.
     T time = 0.0;                                                                       ///< Current time of the simulation.
     T dt = 0.01;
@@ -599,6 +600,10 @@ class Simulation {
      */
     result::SimulationResult<T>* getSimulationResults();
 
+    std::tuple<T, T> getGlobalPressureBounds();
+
+    std::tuple<T, T> getGlobalVelocityBounds();
+
     /**
      * @brief Creates a new fluid out of two existing fluids.
      * @param fluid0Id Id of the first fluid.
@@ -629,6 +634,10 @@ class Simulation {
      * @brief Print the results as pressure at the nodes and flow rates at the channels
      */
     void printResults();
+
+    void writePressurePpm(std::tuple<T, T> bounds, int resolution=600);
+
+    void writeVelocityPpm(std::tuple<T, T> bounds, int resolution=600);
 };
 
 }   // namespace sim
diff --git a/src/simulation/Simulation.hh b/src/simulation/Simulation.hh
@@ -488,6 +488,38 @@ namespace sim {
         return simulationResult.get();
     }
 
+    template<typename T>
+    std::tuple<T, T> Simulation<T>::getGlobalPressureBounds() {
+        T minP = std::numeric_limits<T>::max();
+        T maxP = 0.0;
+        for (auto& [key, simulator] : cfdSimulators) {
+            auto localBounds = simulator->getPressureBounds();
+            if (std::get<0>(localBounds) < minP) {
+                minP = std::get<0>(localBounds);
+            }
+            if (std::get<1>(localBounds) > maxP) {
+                maxP = std::get<1>(localBounds);
+            }
+        }
+        return std::tuple<T, T> {minP, maxP};
+    }
+
+    template<typename T>
+    std::tuple<T, T> Simulation<T>::getGlobalVelocityBounds() {
+        T minVel = std::numeric_limits<T>::max();
+        T maxVel = 0.0;
+        for (auto& [key, simulator] : cfdSimulators) {
+            auto localBounds = simulator->getVelocityBounds();
+            if (std::get<0>(localBounds) < minVel) {
+                minVel = std::get<0>(localBounds);
+            }
+            if (std::get<1>(localBounds) > maxVel) {
+                maxVel = std::get<1>(localBounds);
+            }
+        }
+        return std::tuple<T, T> {minVel, maxVel};
+    }
+
     template<typename T>
     Fluid<T>* Simulation<T>::mixFluids(int fluid0Id, T volume0, int fluid1Id, T volume1) {
         // check if fluids are identically (no merging needed) and if they exist
@@ -593,6 +625,9 @@ namespace sim {
 
             }
 
+            writePressurePpm(getGlobalPressureBounds());
+            writeVelocityPpm(getGlobalVelocityBounds());
+
             #ifdef VERBOSE     
                 if (pressureConverged && allConverged) {
                     std::cout << "[Simulation] All pressures have converged." << std::endl;
@@ -1256,4 +1291,18 @@ namespace sim {
         return events;
     }
 
+    template<typename T>
+    void Simulation<T>::writePressurePpm(std::tuple<T, T> bounds, int resolution) {
+        for (auto& [key, simulator] : cfdSimulators) {
+            simulator->writePressurePpm(0.98*std::get<0>(bounds), 1.02*std::get<1>(bounds), resolution);
+        }
+    }
+
+    template<typename T>
+    void Simulation<T>::writeVelocityPpm(std::tuple<T, T> bounds, int resolution) {
+        for (auto& [key, simulator] : cfdSimulators) {
+            simulator->writeVelocityPpm(0.98*std::get<0>(bounds), 1.02*std::get<1>(bounds), resolution);
+        }
+    }
+
 }   /// namespace sim
diff --git a/src/simulation/simulators/cfdSimulator.h b/src/simulation/simulators/cfdSimulator.h
@@ -269,6 +269,34 @@ class CFDSimulator {
         throw std::runtime_error("The function prepareLattice is undefined for this CFD simulator.");
     }
 
+    /**
+     * @brief Write the vtk file with results of the CFD simulation to file system.
+     * @param[in] iT Iteration step.
+    */
+    virtual void writePressurePpm (T min, T max, int imgResolution=600)
+    {
+        throw std::runtime_error("The function prepareLattice is undefined for this CFD simulator.");
+    }
+
+    /**
+     * @brief Write the vtk file with results of the CFD simulation to file system.
+     * @param[in] iT Iteration step.
+    */
+    virtual void writeVelocityPpm (T min, T max, int imgResolution=600)
+    {
+        throw std::runtime_error("The function prepareLattice is undefined for this CFD simulator.");
+    }
+
+    virtual std::tuple<T, T> getPressureBounds()
+    {
+        throw std::runtime_error("The function prepareLattice is undefined for this CFD simulator.");
+    }
+
+    virtual std::tuple<T, T> getVelocityBounds()
+    {
+        throw std::runtime_error("The function prepareLattice is undefined for this CFD simulator.");
+    }
+
     /**
      * @brief Update the values at the module nodes based on the simulation result after stepIter iterations.
      * @param[in] iT Iteration step.

diff --git a/src/simulation/simulators/olbContinuous.h b/src/simulation/simulators/olbContinuous.h
@@ -184,7 +184,11 @@ using BounceBack = olb::BounceBack<T,DESCRIPTOR>;
      * @brief Write the vtk file with results of the CFD simulation to file system.
      * @param[in] iT Iteration step.
     */
-    void writeVTK(int iT);
+    void writeVTK(int iT) override;
+
+    void writePressurePpm (T min, T max, int imgResolution) override;
+
+    void writeVelocityPpm (T min, T max, int imgResolution) override;
 
     /**
      * @brief Store the abstract pressures at the nodes on the module boundary in the simulator.
@@ -269,6 +273,10 @@ using BounceBack = olb::BounceBack<T,DESCRIPTOR>;
     T getEpsilon() const { 
         return epsilon; 
     };
+
+    std::tuple<T, T> getPressureBounds() override;
+
+    std::tuple<T, T> getVelocityBounds() override;
 };
 
 }   // namespace arch
diff --git a/src/simulation/simulators/olbContinuous.hh b/src/simulation/simulators/olbContinuous.hh
@@ -136,12 +136,33 @@ void lbmSimulator<T>::writeVTK (int iT) {
 
 }
 
+template<typename T>
+void lbmSimulator<T>::writePressurePpm (T min, T max, int imgResolution) {
+    // Color map options are 'earth'|'water'|'air'|'fire'|'leeloo'
+    std::string colorMap = "leeloo";
+    olb::SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure(getLattice(), getConverter());
+    olb::BlockReduction2D2D<T> planeReductionP(pressure, imgResolution);
+    olb::BlockGifWriter<T> ppmWriter(colorMap);
+    ppmWriter.write(planeReductionP, min, max, step, this->name+"_pressure_");
+}
+
+template<typename T>
+void lbmSimulator<T>::writeVelocityPpm (T min, T max, int imgResolution) {
+    // Color map options are 'earth'|'water'|'air'|'fire'|'leeloo'
+    std::string colorMap = "leeloo";
+    olb::SuperLatticePhysVelocity2D<T,DESCRIPTOR> velocity(getLattice(), getConverter());
+    olb::SuperEuklidNorm2D<T, DESCRIPTOR> normVel( velocity );   
+    olb::BlockReduction2D2D<T> planeReductionVel(normVel, imgResolution);
+    olb::BlockGifWriter<T> ppmWriter(colorMap);
+    ppmWriter.write(planeReductionVel, min, max, step, this->name+"_velocity_");
+}
+
 template<typename T>
 void lbmSimulator<T>::solve() {
     int theta = this->updateScheme->getTheta(this->cfdModule->getId());
     this->setBoundaryValues(step);
     for (int iT = 0; iT < theta; ++iT){    
-        writeVTK(step);            
+        writeVTK(step);       
         lattice->collideAndStream();
         step += 1;
     }
@@ -414,4 +435,31 @@ void lbmSimulator<T>::storeCfdResults (int iT) {
     }
 }
 
+template<typename T>
+std::tuple<T, T> lbmSimulator<T>::getPressureBounds() {
+    olb::SuperLatticePhysPressure2D<T,DESCRIPTOR> pressure(getLattice(), getConverter());
+    olb::SuperMin2D<T> minPressureF( pressure, getGeometry(), 1);
+    olb::SuperMax2D<T> maxPressureF( pressure, getGeometry(), 1);
+    int input[0];
+    T minPressure[1];
+    T maxPressure[1];
+    minPressureF( minPressure, input );
+    maxPressureF( maxPressure, input );
+    return std::tuple<T, T> {minPressure[0], maxPressure[0]};
+}
+
+template<typename T>
+std::tuple<T, T> lbmSimulator<T>::getVelocityBounds() {
+    olb::SuperLatticePhysVelocity2D<T,DESCRIPTOR> velocity(getLattice(), getConverter());
+    olb::SuperEuklidNorm2D<T, DESCRIPTOR> normVel( velocity );   
+    olb::SuperMin2D<T> minVelF( normVel, getGeometry(), 1);
+    olb::SuperMax2D<T> maxVelF( normVel, getGeometry(), 1);
+    int input[0];
+    T minVel[1];
+    T maxVel[1];
+    minVelF( minVel, input );
+    maxVelF( maxVel, input );
+    return std::tuple<T, T> {minVel[0], maxVel[0]};
+}
+
 }   // namespace arch