From f70453e2af16a14d8553de134b4a1f481fb5b695 Mon Sep 17 00:00:00 2001
From: Anna Wellmann <a.wellmann@tu-bs.de>
Date: Mon, 11 Jul 2022 10:15:13 +0000
Subject: [PATCH] Add new boundary conditions to simulation and kernelManager

---
 .../GridReaderGenerator/GridGenerator.cpp     | 117 +++++++++++++++++-
 .../KernelManager/BCKernelManager.cpp         |   1 +
 src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp  |  16 ++-
 3 files changed, 126 insertions(+), 8 deletions(-)

diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
index 732482ede..59f45be60 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
@@ -90,6 +90,23 @@ void GridGenerator::allocArrays_BoundaryValues()
 {
     std::cout << "------read BoundaryValues------" << std::endl;
     int blocks = 0;
+
+
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    const auto numberOfPressureValues = int(builder->getPressureSize(0));
+    std::cout << "size pressure: " << numberOfPressureValues << std::endl;
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    blocks                                      = (numberOfPressureValues / para->getParH()->numberofthreads) + 1;
+    para->getParH()->pressureBC.numberOfBCnodes = blocks * para->getParH()->numberofthreads;
+    para->getParD()->pressureBC.numberOfBCnodes = para->getParH()->pressureBC.numberOfBCnodes;
+    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    if (numberOfPressureValues > 1)
+        {
+            cudaMemoryManager->cudaAllocPress();
+            builder->getPressureValues(para->getParH()->pressureBC.RhoBC, para->getParH()->pressureBC.k, para->getParH()->pressureBC.kN, 0);
+            cudaMemoryManager->cudaCopyPress();
+        }
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -137,13 +154,55 @@ void GridGenerator::allocArrays_BoundaryValues()
 
 void GridGenerator::allocArrays_BoundaryQs()
 {
+    
     std::cout << "------read BoundaryQs-------" << std::endl;
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+    
+    const auto numberOfPressureValues = (int)builder->getPressureSize(0);
+    if (numberOfPressureValues > 0)
+    {
+        std::cout << "size Pressure: " << numberOfPressureValues << std::endl;
+        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //preprocessing
+        real* QQ = para->getParH()->pressureBC.q27[0];
+        unsigned int sizeQ = para->getParH()->pressureBC.numberOfBCnodes;
+        QforBoundaryConditions Q;
+        Q.q27[E] = &QQ[E   *sizeQ];
+        Q.q27[W] = &QQ[W   *sizeQ];
+        Q.q27[N] = &QQ[N   *sizeQ];
+        Q.q27[S] = &QQ[S   *sizeQ];
+        Q.q27[T] = &QQ[T   *sizeQ];
+        Q.q27[B] = &QQ[B   *sizeQ];
+        Q.q27[NE] = &QQ[NE  *sizeQ];
+        Q.q27[SW] = &QQ[SW  *sizeQ];
+        Q.q27[SE] = &QQ[SE  *sizeQ];
+        Q.q27[NW] = &QQ[NW  *sizeQ];
+        Q.q27[TE] = &QQ[TE  *sizeQ];
+        Q.q27[BW] = &QQ[BW  *sizeQ];
+        Q.q27[BE] = &QQ[BE  *sizeQ];
+        Q.q27[TW] = &QQ[TW  *sizeQ];
+        Q.q27[TN] = &QQ[TN  *sizeQ];
+        Q.q27[BS] = &QQ[BS  *sizeQ];
+        Q.q27[BN] = &QQ[BN  *sizeQ];
+        Q.q27[TS] = &QQ[TS  *sizeQ];
+        Q.q27[REST] = &QQ[REST*sizeQ];
+        Q.q27[TNE] = &QQ[TNE *sizeQ];
+        Q.q27[TSW] = &QQ[TSW *sizeQ];
+        Q.q27[TSE] = &QQ[TSE *sizeQ];
+        Q.q27[TNW] = &QQ[TNW *sizeQ];
+        Q.q27[BNE] = &QQ[BNE *sizeQ];
+        Q.q27[BSW] = &QQ[BSW *sizeQ];
+        Q.q27[BSE] = &QQ[BSE *sizeQ];
+        Q.q27[BNW] = &QQ[BNW *sizeQ];
+
+        builder->getPressureQs(Q.q27, 0);
+
+        cudaMemoryManager->cudaCopyPress();
+    }
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     const auto numberOfSlipNodes = int(builder->getSlipSize(0));
     if (numberOfSlipNodes > 0) {
-        std::cout << "size velocity: " << numberOfSlipNodes << std::endl;
+        std::cout << "size slip: " << numberOfSlipNodes << std::endl;
         ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
         // preprocessing
         real *QQ           = para->getParH()->slipBC.q27[0];
@@ -181,7 +240,6 @@ void GridGenerator::allocArrays_BoundaryQs()
 
         cudaMemoryManager->cudaCopySlipBC();
     }
-    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     const auto numberOfVelocityNodes = int(builder->getVelocitySize(0));
@@ -226,6 +284,59 @@ void GridGenerator::allocArrays_BoundaryQs()
     }
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
+    const auto numberOfGeometryNodes = builder->getGeometrySize(0);
+    std::cout << "size geometry: " << numberOfGeometryNodes << std::endl;
+    para->getParH()->geometryBC.numberOfBCnodes = numberOfGeometryNodes;
+    para->getParD()->geometryBC.numberOfBCnodes = para->getParH()->geometryBC.numberOfBCnodes;
+    if (numberOfGeometryNodes > 0)
+    {
+        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        cudaMemoryManager->cudaAllocGeomBC();
+        ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+        //////////////////////////////////////////////////////////////////////////
+        //Indexarray
+        builder->getGeometryIndices(para->getParH()->geometryBC.k, 0);
+        //////////////////////////////////////////////////////////////////////////
+        //preprocessing
+        real* QQ = para->getParH()->geometryBC.q27[0];
+        unsigned int sizeQ = para->getParH()->geometryBC.numberOfBCnodes;
+        QforBoundaryConditions Q;
+        Q.q27[E] = &QQ[E   *sizeQ];
+        Q.q27[W] = &QQ[W   *sizeQ];
+        Q.q27[N] = &QQ[N   *sizeQ];
+        Q.q27[S] = &QQ[S   *sizeQ];
+        Q.q27[T] = &QQ[T   *sizeQ];
+        Q.q27[B] = &QQ[B   *sizeQ];
+        Q.q27[NE] = &QQ[NE  *sizeQ];
+        Q.q27[SW] = &QQ[SW  *sizeQ];
+        Q.q27[SE] = &QQ[SE  *sizeQ];
+        Q.q27[NW] = &QQ[NW  *sizeQ];
+        Q.q27[TE] = &QQ[TE  *sizeQ];
+        Q.q27[BW] = &QQ[BW  *sizeQ];
+        Q.q27[BE] = &QQ[BE  *sizeQ];
+        Q.q27[TW] = &QQ[TW  *sizeQ];
+        Q.q27[TN] = &QQ[TN  *sizeQ];
+        Q.q27[BS] = &QQ[BS  *sizeQ];
+        Q.q27[BN] = &QQ[BN  *sizeQ];
+        Q.q27[TS] = &QQ[TS  *sizeQ];
+        Q.q27[REST] = &QQ[REST*sizeQ];
+        Q.q27[TNE] = &QQ[TNE *sizeQ];
+        Q.q27[TSW] = &QQ[TSW *sizeQ];
+        Q.q27[TSE] = &QQ[TSE *sizeQ];
+        Q.q27[TNW] = &QQ[TNW *sizeQ];
+        Q.q27[BNE] = &QQ[BNE *sizeQ];
+        Q.q27[BSW] = &QQ[BSW *sizeQ];
+        Q.q27[BSE] = &QQ[BSE *sizeQ];
+        Q.q27[BNW] = &QQ[BNW *sizeQ];
+        
+        builder->getGeometryQs(Q.q27,0);
+        
+        for (uint node_i = 0; node_i < numberOfGeometryNodes; node_i++)
+        {
+            Q.q27[REST][node_i] = 0.0f;
+        }
+        cudaMemoryManager->cudaCopyGeomBC();
+    }
     ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
     std::cout << "-----finish BoundaryQs------" << std::endl;
 }
diff --git a/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp b/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp
index b58c59376..586e817e0 100644
--- a/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp
+++ b/src/gpu/VirtualFluids_GPU/KernelManager/BCKernelManager.cpp
@@ -32,6 +32,7 @@
 //=======================================================================================
 #include <cuda_runtime.h>
 #include <helper_cuda.h>
+#include <iostream>
 
 #include "BCKernelManager.h"
 #include "Parameter/Parameter.h"
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index d9d8f0e5b..77a4eec97 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -146,14 +146,20 @@ void Simulation::run()
         lbKernelManager->runLBMKernel(para);
 
         ////////////////////////////////////////////////////////////////////////////////
-        // velocity boundary condition
+        // run post-collision boundary conditions
         this->bcKernelManager->runVelocityBCKernelPost(level);
+        this->bcKernelManager->runNoSlipBCKernelPost(level);
+        this->bcKernelManager->runSlipBCKernelPost(level);
+        this->bcKernelManager->runGeoBCKernelPost(level);
 
         ////////////////////////////////////////////////////////////////////////////////
-        if (para->getParD()->isEvenTimestep)
-            para->getParD()->isEvenTimestep = false;
-        else
-            para->getParD()->isEvenTimestep = true;
+        // swap between even and odd timestep
+        para->getParD()->isEvenTimestep = ! para->getParD()->isEvenTimestep;
+
+        ////////////////////////////////////////////////////////////////////////////////
+        //  run pre-collision boundary conditions
+        this->bcKernelManager->runPressureBCKernelPre(level);
+
         ////////////////////////////////////////////////////////////////////////////////
         // File IO and calculation of MNUPS(million node updates per second)
         if (para->getTimestepOut() > 0 && timestep % para->getTimestepOut() == 0 &&
-- 
GitLab