From 18fb32d23409c3a97b0e910eb183a479a078a27b Mon Sep 17 00:00:00 2001
From: "LEGOLAS\\lenz" <lenz@irmb.tu-bs.de>
Date: Fri, 26 Jul 2019 13:49:20 +0200
Subject: [PATCH] modifies PointTimeSeriesAnalyzer

---
 .../Analyzer/PointTimeSeriesAnalyzer.cu       | 14 ++--
 src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.h | 14 ++--
 .../Analyzer/PointTimeSeriesCollector.cpp     | 80 +++++++++++++++++++
 .../Analyzer/PointTimeseriesCollector.h       | 39 +++++++++
 .../GKS/RoomFireExtended/RoomFireExtended.cpp | 72 ++++++++++-------
 5 files changed, 177 insertions(+), 42 deletions(-)
 create mode 100644 src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp
 create mode 100644 src/GksGpu/Analyzer/PointTimeseriesCollector.h

diff --git a/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu b/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu
index 90f2c805e..b86dc72e5 100644
--- a/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu
+++ b/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.cu
@@ -16,6 +16,8 @@
 
 #include "Core/Logger/Logger.h"
 
+#include "GksMeshAdapter/GksMeshAdapter.h"
+
 #include "DataBase/DataBase.h"
 
 #include "Parameters/Parameters.h"
@@ -32,7 +34,7 @@ __host__ __device__ inline void pointTimeSeriesFunction( DataBaseStruct dataBase
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-bool PointTimeSeriesAnalyzer::run(uint iter, Parameters parameters)
+void PointTimeSeriesAnalyzer::run(uint iter, Parameters parameters)
 {
 
     CudaUtility::CudaGrid grid( 1, 1 );
@@ -99,7 +101,7 @@ PointTimeSeriesAnalyzer::~PointTimeSeriesAnalyzer()
     this->free();
 }
 
-PointTimeSeriesAnalyzer::PointTimeSeriesAnalyzer(SPtr<DataBase> dataBase, GksMeshAdapter & adapter, Vec3 coordinate, char quantity, uint outputIter)
+PointTimeSeriesAnalyzer::PointTimeSeriesAnalyzer(SPtr<DataBase> dataBase, GksMeshAdapter & adapter, Vec3 coordinates, char quantity, uint outputIter)
     : dataBase(dataBase),
       deviceSeries(nullptr),
       counter(0),
@@ -108,7 +110,7 @@ PointTimeSeriesAnalyzer::PointTimeSeriesAnalyzer(SPtr<DataBase> dataBase, GksMes
 {
     this->allocate();
 
-    this->findCellIndex( adapter, coordinate );
+    this->findCellIndex( adapter, coordinates );
 }
 
 void PointTimeSeriesAnalyzer::free()
@@ -123,7 +125,7 @@ void PointTimeSeriesAnalyzer::allocate()
     checkCudaErrors( cudaMalloc ( &this->deviceSeries , sizeof(real) * this->outputIter ) );
 }
 
-void PointTimeSeriesAnalyzer::findCellIndex(GksMeshAdapter & adapter, Vec3 coordinate)
+void PointTimeSeriesAnalyzer::findCellIndex(GksMeshAdapter & adapter, Vec3 coordinates)
 {
     real minDistance = 1.0e99;
 
@@ -131,7 +133,7 @@ void PointTimeSeriesAnalyzer::findCellIndex(GksMeshAdapter & adapter, Vec3 coord
     {
         MeshCell& cell = adapter.cells[ cellIdx ];
 
-        Vec3 vec = cell.cellCenter - coordinate;
+        Vec3 vec = cell.cellCenter - coordinates;
 
         real distance = sqrt( vec.x*vec.x + vec.y*vec.y + vec.z*vec.z );
 
@@ -142,6 +144,8 @@ void PointTimeSeriesAnalyzer::findCellIndex(GksMeshAdapter & adapter, Vec3 coord
         }
     }
 
+    this->coordinates = adapter.cells[ this->cellIndex ].cellCenter;
+
     *logging::out << logging::Logger::INFO_INTERMEDIATE << "PointTimeSeriesAnalyzer::cellIndex = " << this->cellIndex << "\n";
 }
 
diff --git a/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.h b/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.h
index 661c4cac9..867ae1fd6 100644
--- a/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.h
+++ b/src/GksGpu/Analyzer/PointTimeSeriesAnalyzer.h
@@ -9,13 +9,11 @@
 
 #include "Core/PointerDefinitions.h"
 #include "Core/DataTypes.h"
-
-#include "GksMeshAdapter/GksMeshAdapter.h"
-
-#include "FlowStateData/FlowStateData.cuh"
+#include "Core/VectorTypes.h"
 
 struct DataBase;
 struct Parameters;
+class  GksMeshAdapter;
 
 struct PointTimeSeriesAnalyzerStruct
 {
@@ -46,19 +44,21 @@ public:
 
     std::vector<real> hostSeries;
 
+    Vec3 coordinates;
+
 public:
 
     ~PointTimeSeriesAnalyzer();
 
-    PointTimeSeriesAnalyzer( SPtr<DataBase> dataBase, GksMeshAdapter & adapter, Vec3 coordinate, char quantity, uint outputIter = 10000 );
+    PointTimeSeriesAnalyzer( SPtr<DataBase> dataBase, GksMeshAdapter & adapter, Vec3 coordinates, char quantity, uint outputIter = 10000 );
 
     void free();
 
     void allocate();
 
-    void findCellIndex( GksMeshAdapter & adapter, Vec3 coordinate );
+    void findCellIndex( GksMeshAdapter & adapter, Vec3 coordinates );
 
-    bool run( uint iter, Parameters parameters );
+    void run( uint iter, Parameters parameters );
 
     void writeToFile( std::string filename );
 
diff --git a/src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp b/src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp
new file mode 100644
index 000000000..0c2a8ece4
--- /dev/null
+++ b/src/GksGpu/Analyzer/PointTimeSeriesCollector.cpp
@@ -0,0 +1,80 @@
+#include "PointTimeseriesCollector.h"
+
+#include <iomanip>
+#include <fstream>
+
+#include "Core/Logger/Logger.h"
+
+#include "Analyzer/PointTimeseriesAnalyzer.h"
+
+#include "Parameters/Parameters.h"
+
+PointTimeSeriesCollector::~PointTimeSeriesCollector()
+{
+}
+
+PointTimeSeriesCollector::PointTimeSeriesCollector()
+{
+}
+
+void PointTimeSeriesCollector::addAnalyzer(SPtr<DataBase> dataBase, GksMeshAdapter & adapter, Vec3 coordinate, char quantity, uint outputIter)
+{
+    auto pointTimeSeriesAnalyzer = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, adapter, coordinate, quantity, outputIter );
+
+    this->analyzerList.push_back( pointTimeSeriesAnalyzer );
+}
+
+void PointTimeSeriesCollector::run(uint iter, Parameters parameters)
+{
+    for( auto analyzer : this->analyzerList )
+        analyzer->run(iter, parameters);
+}
+
+void PointTimeSeriesCollector::writeToFile(std::string filename)
+{
+    *logging::out << logging::Logger::INFO_INTERMEDIATE << "PointTimeSeriesCollector::writeToFile( " << filename << " )" << "\n";
+
+    if( this->analyzerList.size() == 0 )
+    {
+        *logging::out << logging::Logger::WARNING << "empty!\n";
+        return;
+    }
+
+    std::ofstream file;
+
+    file.open(filename + ".dat" );
+
+    //////////////////////////////////////////////////////////////////////////
+
+    file << "Number of Points = " << this->analyzerList.size() << "\n";
+
+    for( uint j = 0; j < this->analyzerList.size(); j++ )
+    {
+        file << "Point " << j << ", ";
+        file << "Quantity = "     << this->analyzerList[j]->quantity << ", ";
+        file << "Coordinates = ( " << this->analyzerList[j]->coordinates.x << ", "
+                                   << this->analyzerList[j]->coordinates.y << ", "
+                                   << this->analyzerList[j]->coordinates.z << " )";
+        file << "\n";
+    }
+
+    //////////////////////////////////////////////////////////////////////////
+
+    uint numberOfTimeSteps = this->analyzerList[0]->hostSeries.size();
+
+    for( uint i = 0; i < numberOfTimeSteps; i++ )
+    {
+        for( uint j = 0; j < this->analyzerList.size(); j++ )
+        {
+            file << std::setprecision(15) << this->analyzerList[j]->hostSeries[i] << ", ";
+        }
+
+        file << "\n";
+    }
+
+    //////////////////////////////////////////////////////////////////////////
+
+    file.close();
+
+    *logging::out << logging::Logger::INFO_INTERMEDIATE << "done!\n";
+}
diff --git a/src/GksGpu/Analyzer/PointTimeseriesCollector.h b/src/GksGpu/Analyzer/PointTimeseriesCollector.h
new file mode 100644
index 000000000..4ed4c6ab8
--- /dev/null
+++ b/src/GksGpu/Analyzer/PointTimeseriesCollector.h
@@ -0,0 +1,39 @@
+#ifndef  PointTimeSeriesCelloctor_H
+#define  PointTimeSeriesCelloctor_H
+
+
+#include <vector>
+#include <string>
+#include <memory>
+
+#include "VirtualFluidsDefinitions.h"
+
+#include "Core/PointerDefinitions.h"
+#include "Core/DataTypes.h"
+#include "Core/VectorTypes.h"
+
+class  PointTimeSeriesAnalyzer;
+struct DataBase;
+struct Parameters;
+class  GksMeshAdapter;
+
+class VF_PUBLIC PointTimeSeriesCollector
+{
+public:
+
+    std::vector< SPtr<PointTimeSeriesAnalyzer> > analyzerList;
+
+public:
+
+    ~PointTimeSeriesCollector();
+
+    PointTimeSeriesCollector(  );
+
+    void addAnalyzer( SPtr<DataBase> dataBase, GksMeshAdapter & adapter, Vec3 coordinate, char quantity, uint outputIter = 10000 );
+
+    void run( uint iter, Parameters parameters );
+
+    void writeToFile( std::string filename );
+};
+
+#endif
diff --git a/targets/apps/GKS/RoomFireExtended/RoomFireExtended.cpp b/targets/apps/GKS/RoomFireExtended/RoomFireExtended.cpp
index da290ce18..cc327901f 100644
--- a/targets/apps/GKS/RoomFireExtended/RoomFireExtended.cpp
+++ b/targets/apps/GKS/RoomFireExtended/RoomFireExtended.cpp
@@ -50,7 +50,7 @@
 #include "GksGpu/Analyzer/CupsAnalyzer.h"
 #include "GksGpu/Analyzer/ConvergenceAnalyzer.h"
 #include "GksGpu/Analyzer/TurbulenceAnalyzer.h"
-#include "GksGpu/Analyzer/PointTimeseriesAnalyzer.h"
+#include "GksGpu/Analyzer/PointTimeseriesCollector.h"
 
 #include "GksGpu/Restart/Restart.h"
 
@@ -82,11 +82,19 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
     real cs  = sqrt( ( ( K + 5.0 ) / ( K + 3.0 ) ) / ( 2.0 * prim.lambda ) );
 
     real mu      = 1.8e-5;
-    real U       = 0.0125;       // 900 kW on top
+    //real U       = 0.025;       // 750 kW on top
     //real U       = 0.015;       // 900 kW on top
     //real U       = 0.005;       // 900 kW all around
     real rhoFuel = 0.5405;
 
+    real heatOfReaction = real(8000.0); // J / mol 
+
+    real specificHeatOfReaction = heatOfReaction / 0.016;
+
+    real HRR = 750.0; // kW
+
+    real U = HRR * 1000.0 / ( rhoFuel * LBurner * LBurner * (specificHeatOfReaction * 100.0) );
+
     real CFL = 0.125;
 
     real dt  = CFL * ( dx / ( ( U + cs ) * ( c1o1 + ( c2o1 * mu ) / ( U * dx * rho ) ) ) );
@@ -96,7 +104,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
     *logging::out << logging::Logger::INFO_HIGH << "cs = " << cs << " m/s\n";
     *logging::out << logging::Logger::INFO_HIGH << "mu = " << mu << " kg/sm\n";
 
-    *logging::out << logging::Logger::INFO_HIGH << "HRR = " << U * rho * LBurner * LBurner * 800000.0 / 0.016 / 1000.0 << " kW\n";
+    *logging::out << logging::Logger::INFO_HIGH << "HRR = " << U * rhoFuel * LBurner * LBurner * (heatOfReaction * 100.0) / 0.016 / 1000.0 << " kW\n";
 
     //////////////////////////////////////////////////////////////////////////
 
@@ -122,6 +130,8 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
 
     parameters.enableReaction = true;
 
+    parameters.heatOfReaction = heatOfReaction;
+
     parameters.useHeatReleaseRateLimiter = true;
     parameters.useTemperatureLimiter     = true;
     parameters.usePassiveScalarLimiter   = true;
@@ -145,8 +155,11 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
-    gridBuilder->addCoarseGrid(-2.1, -1.6, -0.1,  
-                                2.1,  6.0,  5.0, dx);
+    //gridBuilder->addCoarseGrid(-2.1, -1.6, -0.1,  
+    //                            2.1,  6.0,  5.0, dx);
+
+    gridBuilder->addCoarseGrid(-2.1, -0.6, -0.1,  
+                                2.1,  0.6,  5.0, dx);
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
@@ -241,7 +254,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
 
     //////////////////////////////////////////////////////////////////////////
 
-    SPtr<BoundaryCondition> bcBurner = std::make_shared<CreepingMassFlux>( dataBase, rho, U, prim.lambda );
+    SPtr<BoundaryCondition> bcBurner = std::make_shared<CreepingMassFlux>( dataBase, rhoFuel, U, prim.lambda );
 
     bcBurner->findBoundaryCells( meshAdapter, false, [&](Vec3 center){ 
 
@@ -271,13 +284,24 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
 
     //////////////////////////////////////////////////////////////////////////
 
-    auto pointTimeSeriesAnalyzer_P1 = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, meshAdapter, Vec3(-1.5, 0.0, 2.5999), 'T' );
-    auto pointTimeSeriesAnalyzer_P2 = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, meshAdapter, Vec3(-1.0, 0.0, 2.5999), 'T' );
-    auto pointTimeSeriesAnalyzer_P3 = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, meshAdapter, Vec3(-0.5, 0.0, 2.5999), 'T' );
-    auto pointTimeSeriesAnalyzer_P4 = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, meshAdapter, Vec3( 0.0, 0.0, 2.5999), 'T' );
-    auto pointTimeSeriesAnalyzer_P5 = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, meshAdapter, Vec3( 0.5, 0.0, 2.5999), 'T' );
-    auto pointTimeSeriesAnalyzer_P6 = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, meshAdapter, Vec3( 1.0, 0.0, 2.5999), 'T' );
-    auto pointTimeSeriesAnalyzer_P7 = std::make_shared<PointTimeSeriesAnalyzer>( dataBase, meshAdapter, Vec3( 1.5, 0.0, 2.5999), 'T' );
+    auto pointTimeSeriesCollector = std::make_shared<PointTimeSeriesCollector>();
+
+    for( real x = -0.0001; x < 2; x += 0.5 )
+    {
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x, -1.4999, 2.9999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x, -1.0,    2.9999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x, -0.5,    2.9999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x, -0.2001, 2.9999), 'T' );
+
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x, -0.2001, 2.5999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x,  0.0,    2.5999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x,  0.2001, 2.5999), 'T' );
+        
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x,  0.2001, 2.9999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x,  0.5,    2.9999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x,  1.0,    2.9999), 'T' );
+        pointTimeSeriesCollector->addAnalyzer( dataBase, meshAdapter, Vec3(x,  1.4999, 2.9999), 'T' );
+    }
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -339,13 +363,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
 
         TimeStepping::nestedTimeStep(dataBase, parameters, 0);
 
-        pointTimeSeriesAnalyzer_P1->run(iter, parameters);
-        pointTimeSeriesAnalyzer_P2->run(iter, parameters);
-        pointTimeSeriesAnalyzer_P3->run(iter, parameters);
-        pointTimeSeriesAnalyzer_P4->run(iter, parameters);
-        pointTimeSeriesAnalyzer_P5->run(iter, parameters);
-        pointTimeSeriesAnalyzer_P6->run(iter, parameters);
-        pointTimeSeriesAnalyzer_P7->run(iter, parameters);
+        pointTimeSeriesCollector->run(iter, parameters);
 
         int crashCellIndex = dataBase->getCrashCellIndex();
         if( crashCellIndex >= 0 )
@@ -357,27 +375,21 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
             break;
         }
 
-        if( iter % 5000 == 0 )
+        if( iter % 10000 == 0 )
         {
             dataBase->copyDataDeviceToHost();
 
             writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) );
         }
 
-        if( iter % 5000 == 0 )
+        if( iter % 10000 == 0 )
         {
             Restart::writeRestart( dataBase, path + simulationName + "_" + std::to_string( iter ), iter );
         }
 
-        if( iter % 5000 == 0 )
+        if( iter % 10000 == 0 )
         {
-            pointTimeSeriesAnalyzer_P1->writeToFile(path + simulationName + "_P1_TimeSeries_" + std::to_string( iter ));
-            pointTimeSeriesAnalyzer_P2->writeToFile(path + simulationName + "_P2_TimeSeries_" + std::to_string( iter ));
-            pointTimeSeriesAnalyzer_P3->writeToFile(path + simulationName + "_P3_TimeSeries_" + std::to_string( iter ));
-            pointTimeSeriesAnalyzer_P4->writeToFile(path + simulationName + "_P4_TimeSeries_" + std::to_string( iter ));
-            pointTimeSeriesAnalyzer_P5->writeToFile(path + simulationName + "_P5_TimeSeries_" + std::to_string( iter ));
-            pointTimeSeriesAnalyzer_P6->writeToFile(path + simulationName + "_P6_TimeSeries_" + std::to_string( iter ));
-            pointTimeSeriesAnalyzer_P7->writeToFile(path + simulationName + "_P7_TimeSeries_" + std::to_string( iter ));
+            pointTimeSeriesCollector->writeToFile(path + simulationName + "_TimeSeries_" + std::to_string( iter ));
         }
 
         //turbulenceAnalyzer->run( iter, parameters );
-- 
GitLab