From 2a1b2f3dd18d57e0a1fe59babb16826bfbe73638 Mon Sep 17 00:00:00 2001
From: Anna Wellmann <a.wellmann@tu-braunschweig.de>
Date: Fri, 10 Nov 2023 13:09:25 +0000
Subject: [PATCH] Move initial rotation into separate function

---
 apps/gpu/RotatingGrid/RotatingGrid.cpp        | 56 ++++++++++++-------
 src/gpu/core/Calculation/UpdateGrid27.cpp     | 45 ++++++++-------
 src/gpu/core/Calculation/UpdateGrid27.h       |  1 +
 src/gpu/core/GPU/GPU_Interface.h              | 14 ++---
 src/gpu/core/GPU/LBMKernel.cu                 | 12 ++--
 src/gpu/core/LBM/Simulation.cpp               |  4 ++
 .../core/Parameter/ParameterRotatingGrid.h    | 22 +++++---
 7 files changed, 91 insertions(+), 63 deletions(-)

diff --git a/apps/gpu/RotatingGrid/RotatingGrid.cpp b/apps/gpu/RotatingGrid/RotatingGrid.cpp
index 84b96bc81..9498bc629 100644
--- a/apps/gpu/RotatingGrid/RotatingGrid.cpp
+++ b/apps/gpu/RotatingGrid/RotatingGrid.cpp
@@ -73,17 +73,28 @@
 int main()
 {
     try {
-         vf::logging::Logger::initializeLogger();
+        vf::logging::Logger::initializeLogger();
         //////////////////////////////////////////////////////////////////////////
         // Simulation parameters
         //////////////////////////////////////////////////////////////////////////
-        enum RotationOrInterpolation {Rot, Int};
+        enum RotationOrInterpolation { Rot, Int };
+        enum InitialRotation { Zero, TwoPi, Pi, PiHalf };
+        using namespace vf::basics::constant;
+        std::map<InitialRotation, real> rotations = { { Zero, 0.0 }, { TwoPi, c2Pi }, { Pi, cPi }, { PiHalf, PiHalf } };
+        std::map<InitialRotation, std::string> rotationStrings = {
+            { Zero, "NoRotation" }, { TwoPi, "2Pi" }, { Pi, "Pi" }, { PiHalf, "PiHalf" }
+        };
+
         const RotationOrInterpolation rotOrInt = Rot;
-        if (rotOrInt == Int) VF_LOG_INFO("Use interpolation.");
-        if (rotOrInt == Rot) VF_LOG_INFO("Use rotation.");
+        const InitialRotation initialRot = PiHalf;
+
+        if (rotOrInt == Int)
+            VF_LOG_INFO("Use interpolation.");
+        if (rotOrInt == Rot)
+            VF_LOG_INFO("Use rotation.");
 
         const std::string path("./output/RotatingGrid");
-        const std::string simulationName = rotOrInt == Int ? "RotatingGridInterpolationTest" : "RotatingGrid"; //PiHalbe";
+        const std::string simulationName = (rotOrInt == Int ? "RotatingGridInterpolationTest" : "RotatingGrid") + rotationStrings[initialRot];
 
         const real L = 1.0;
         const real Re = 2000.0;
@@ -100,7 +111,7 @@ int main()
         //////////////////////////////////////////////////////////////////////////
 
         const real dx = L / real(nx);
-        const real dt  = velocityLB / velocity * dx;
+        const real dt = velocityLB / velocity * dx;
 
         const real viscosityLB = nx * velocityLB / Re; // LB units
 
@@ -111,8 +122,10 @@ int main()
 
         gridBuilder->addCoarseGrid(-0.5 * L, -0.5 * L, -0.5 * L, 0.5 * L, 0.5 * L, 0.5 * L, dx);
 
-        if (rotOrInt == Rot) gridBuilder->addGridRotatingGrid(std::make_shared<Cylinder>(0.0, 0.0, 0.0, 0.25 * L, 0.75 * L, Axis::z));
-        if (rotOrInt == Int) gridBuilder->addGrid(std::make_shared<Cylinder>(0.2, 0.1, 0.1, 0.25 * L, 0.8 * L, Axis::x), 1);
+        if (rotOrInt == Rot)
+            gridBuilder->addGridRotatingGrid(std::make_shared<Cylinder>(0.0, 0.0, 0.0, 0.25 * L, 0.75 * L, Axis::x));
+        if (rotOrInt == Int)
+            gridBuilder->addGrid(std::make_shared<Cylinder>(0.2, 0.1, 0.1, 0.25 * L, 0.8 * L, Axis::x), 1);
 
         GridScalingFactory scalingFactory = GridScalingFactory();
         scalingFactory.setScalingFactory(GridScalingFactory::GridScaling::ScaleCompressible);
@@ -152,13 +165,17 @@ int main()
         // gridBuilder->setSlipBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
         // gridBuilder->setSlipBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0);
 
-        gridBuilder->setNoSlipBoundaryCondition(SideType::MY);
-        gridBuilder->setNoSlipBoundaryCondition(SideType::PY);
         gridBuilder->setNoSlipBoundaryCondition(SideType::MZ);
         gridBuilder->setNoSlipBoundaryCondition(SideType::PZ);
 
-        gridBuilder->setVelocityBoundaryCondition(SideType::PX, -velocityLB, 0.0, 0.0);
-        gridBuilder->setPressureBoundaryCondition(SideType::MX, 0.0);
+        // gridBuilder->setVelocityBoundaryCondition(SideType::PX, -velocityLB, 0.0, 0.0);
+        // gridBuilder->setPressureBoundaryCondition(SideType::MX, 0.0);
+
+        // rot around x-Axis
+        gridBuilder->setNoSlipBoundaryCondition(SideType::MX);
+        gridBuilder->setNoSlipBoundaryCondition(SideType::PX);
+        gridBuilder->setVelocityBoundaryCondition(SideType::PY, 0.0, -velocityLB, 0.0);
+        gridBuilder->setPressureBoundaryCondition(SideType::MY, 0.0);
 
         BoundaryConditionFactory bcFactory;
 
@@ -185,10 +202,8 @@ int main()
 
         // para->setInitialCondition([&](real coordX, real coordY, real coordZ, real &rho, real &vx, real &vy, real &vz) {
         //     // rho = (real) 0.0;
-        //     // if (coordX > -0.52 && coordY > 0.27 && coordZ > -0.1 && coordX < -0.49 && coordY < 0.3 && coordZ < 0.11) rho = 1e-5;
-        //     rho =  setPressPoint({-0.2, 0.24, 0.0}, 1e-5, dx, {coordX, coordY, coordZ});
-        //     vx = 0.0;
-        //     vy = 0.0;
+        //     // if (coordX > -0.52 && coordY > 0.27 && coordZ > -0.1 && coordX < -0.49 && coordY < 0.3 && coordZ < 0.11)
+        //     rho = 1e-5; rho =  setPressPoint({-0.2, 0.24, 0.0}, 1e-5, dx, {coordX, coordY, coordZ}); vx = 0.0; vy = 0.0;
         //     vz = 0.0;
         // });
 
@@ -196,7 +211,7 @@ int main()
         // set copy mesh to simulation
         //////////////////////////////////////////////////////////////////////////
 
-        vf::parallel::Communicator &communicator = *vf::parallel::MPICommunicator::getInstance();
+        vf::parallel::Communicator& communicator = *vf::parallel::MPICommunicator::getInstance();
 
         auto cudaMemoryManager = std::make_shared<CudaMemoryManager>(para);
         SPtr<GridProvider> gridGenerator =
@@ -238,13 +253,14 @@ int main()
         // NeighborDebugWriter::writeNeighborLinkLinesDebug(para.get());
 
         SPtr<ParameterRotatingGrid> paraRot = para->getRotatingGridParameter();
+        paraRot->initialGridRotation = rotations[initialRot];
         sim.run();
 
-    } catch (const spdlog::spdlog_ex &ex) {
+    } catch (const spdlog::spdlog_ex& ex) {
         std::cout << "Log initialization failed: " << ex.what() << std::endl;
-    } catch (const std::bad_alloc &e) {
+    } catch (const std::bad_alloc& e) {
         VF_LOG_CRITICAL("Bad Alloc: {}", e.what());
-    } catch (const std::exception &e) {
+    } catch (const std::exception& e) {
         VF_LOG_CRITICAL("exception: {}", e.what());
     } catch (...) {
         VF_LOG_CRITICAL("Unknown exception!");
diff --git a/src/gpu/core/Calculation/UpdateGrid27.cpp b/src/gpu/core/Calculation/UpdateGrid27.cpp
index 95ac0106f..243642679 100644
--- a/src/gpu/core/Calculation/UpdateGrid27.cpp
+++ b/src/gpu/core/Calculation/UpdateGrid27.cpp
@@ -427,18 +427,24 @@ UpdateGrid27::UpdateGrid27(SPtr<Parameter> para, vf::parallel::Communicator &com
     this->bcKernelManager = std::make_shared<BCKernelManager>(para, bcFactory);
     this->adKernelManager = std::make_shared<ADKernelManager>(para);
     this->gridScalingKernelManager = std::make_shared<GridScalingKernelManager>(para, scalingFactory);
+}
 
+void UpdateGrid27::rotateGridInInitializationProcess(int level)
+{
+    real steps = 800; // to few steps may break the connectivity of the interpolation cells
     SPtr<ParameterRotatingGrid> paraRot = para->getRotatingGridParameter();
 
-    uint level = 0;
-    real finalRotation = (real) (vf::basics::constant::cPi);
-    VF_LOG_INFO("intended rotation: {} around center point ({}, {}, {})", finalRotation,
-                paraRot->parameterRotHost->centerPoint[0], paraRot->parameterRotHost->centerPoint[1],
-                paraRot->parameterRotHost->centerPoint[2]);
-    real steps = 800;
-    paraRot->parameterRotHost->angularVelocity[0] = finalRotation / steps;
-    paraRot->parameterRotDevice->angularVelocity[0] = paraRot->parameterRotHost->angularVelocity[0];
-    VF_LOG_INFO("angular velocity: {}", paraRot->parameterRotHost->angularVelocity[0]);
+    real intendedRotation = (real)paraRot->initialGridRotation;
+    auto axis = paraRot->rotationalAxis;
+    auto angularVelocityInSimulation = paraRot->parameterRotHost->angularVelocity;
+    paraRot->parameterRotHost->angularVelocity[axis] = intendedRotation / steps;
+    paraRot->parameterRotDevice->angularVelocity[axis] = paraRot->parameterRotHost->angularVelocity[axis];
+
+
+    VF_LOG_INFO("Intended initial rotation for the rotating grid: {}. Center point for the rotation: ({}, {}, {}), Axis : {}, Angular Velocity: {}.",
+                intendedRotation, paraRot->parameterRotHost->centerPoint[0], paraRot->parameterRotHost->centerPoint[1],
+                paraRot->parameterRotHost->centerPoint[2], axis, paraRot->parameterRotHost->angularVelocity[axis]);
+
     for (int i = 0; i < steps; i++) {
         paraRot->parameterRotHost->gridAngle[0] += paraRot->parameterRotHost->angularVelocity[0];
         paraRot->parameterRotHost->gridAngle[1] += paraRot->parameterRotHost->angularVelocity[1];
@@ -450,7 +456,7 @@ UpdateGrid27::UpdateGrid27(SPtr<Parameter> para, vf::parallel::Communicator &com
         // base to nested
         TraverseStaticToRotating(
             para->getParD(level).get(),
-            para->getParD(level+1).get(),
+            para->getParD(level + 1).get(),
             para->getRotatingGridParameter()->parameterRotDevice.get(),
             &para->getParD(level)->coarseToFine,
             para->getParD(level)->neighborCoarseToFine);
@@ -458,22 +464,19 @@ UpdateGrid27::UpdateGrid27(SPtr<Parameter> para, vf::parallel::Communicator &com
         // nested to base
         TraverseRotatingToStatic(
             para->getParD(level).get(),
-            para->getParD(level+1).get(),
+            para->getParD(level + 1).get(),
             para->getRotatingGridParameter()->parameterRotDevice.get(),
             &para->getParD(level)->fineToCoarse,
             para->getParD(level)->neighborFineToCoarse);
     }
-    paraRot->parameterRotHost->angularVelocity[0] = 0.0;
-    paraRot->parameterRotHost->angularVelocity[1] = 0.0;
-    paraRot->parameterRotHost->angularVelocity[2] = 0.0;
-    paraRot->parameterRotDevice->angularVelocity[0] = 0.0;
-    paraRot->parameterRotDevice->angularVelocity[1] = 0.0;
-    paraRot->parameterRotDevice->angularVelocity[2] = 0.0;
+
+    paraRot->parameterRotHost->angularVelocity[axis] = angularVelocityInSimulation[axis];
+    paraRot->parameterRotDevice->angularVelocity[axis] = angularVelocityInSimulation[axis];
 
     VF_LOG_INFO("Angle x {}", paraRot->parameterRotHost->gridAngle[0]);
     VF_LOG_INFO("Angle y {}", paraRot->parameterRotHost->gridAngle[1]);
     VF_LOG_INFO("Angle z {}", paraRot->parameterRotHost->gridAngle[2]);
-    VF_LOG_INFO("Angular velocity x {}", paraRot->parameterRotHost->angularVelocity[0]);
-    VF_LOG_INFO("Angular velocity y {}", paraRot->parameterRotHost->angularVelocity[1]);
-    VF_LOG_INFO("Angular velocity z {}", paraRot->parameterRotHost->angularVelocity[2]);
-}
+    VF_LOG_INFO("In simulation: Angular velocity x {}", paraRot->parameterRotHost->angularVelocity[0]);
+    VF_LOG_INFO("In simulation: Angular velocity y {}", paraRot->parameterRotHost->angularVelocity[1]);
+    VF_LOG_INFO("In simulation: Angular velocity z {}", paraRot->parameterRotHost->angularVelocity[2]);
+}
\ No newline at end of file
diff --git a/src/gpu/core/Calculation/UpdateGrid27.h b/src/gpu/core/Calculation/UpdateGrid27.h
index 5218fdf80..f6ba49678 100644
--- a/src/gpu/core/Calculation/UpdateGrid27.h
+++ b/src/gpu/core/Calculation/UpdateGrid27.h
@@ -33,6 +33,7 @@ public:
     void updateGrid(int level, unsigned int t);
     void exchangeData(int level);
 
+    void rotateGridInInitializationProcess(int baseLevelForInterpolation);
 private:
     void collisionAllNodes(int level, unsigned int t);
     void collisionUsingIndices(int level, unsigned int t, uint *taggedFluidNodeIndices = nullptr, uint numberOfTaggedFluidNodes = 0, CollisionTemplate collisionTemplate = CollisionTemplate::Default, CudaStreamIndex streamIndex=CudaStreamIndex::Legacy);
diff --git a/src/gpu/core/GPU/GPU_Interface.h b/src/gpu/core/GPU/GPU_Interface.h
index e1608160c..ba9b1e4cf 100644
--- a/src/gpu/core/GPU/GPU_Interface.h
+++ b/src/gpu/core/GPU/GPU_Interface.h
@@ -25,7 +25,7 @@
 #endif
 
 struct LBMSimulationParameter;
-struct ParameterRotatingGridSimulation;
+struct ParameterRotatingGridHostDevice;
 class Parameter;
 
 //////////////////////////////////////////////////////////////////////////
@@ -403,7 +403,7 @@ void CalcMacCompSP27RotatingToStatic(
     unsigned int numberOfThreads,
     real* DD,
     bool isEvenTimestep,
-    ParameterRotatingGridSimulation* parameterRotDevice);
+    ParameterRotatingGridHostDevice* parameterRotDevice);
 
 void CalcMacThS7(  real* Conc,
                               unsigned int* geoD,
@@ -1642,34 +1642,34 @@ template<bool hasTurbulentViscosity> void ScaleCF_compressible(LBMSimulationPara
 void InterpolateStaticToRotating(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *baseToNested,
     ICellNeigh &neighborBaseToNested);
 
 void TraverseStaticToRotating(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *baseToNested,
     ICellNeigh &neighborBaseToNested);
 
 void InterpolateRotatingToStatic(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *nestedToBase,
     ICellNeigh &neighborNestedToBase);
 
 void TraverseRotatingToStatic(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *nestedToBase,
     ICellNeigh &neighborNestedToBase);
 
 void UpdateGlobalCoordinates(
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice
+    ParameterRotatingGridHostDevice *paraRotDevice
 );
 
 void ScaleCF_RhoSq_3rdMom_comp_27( real* DC, 
diff --git a/src/gpu/core/GPU/LBMKernel.cu b/src/gpu/core/GPU/LBMKernel.cu
index 5bb4ada3a..bd4feab17 100644
--- a/src/gpu/core/GPU/LBMKernel.cu
+++ b/src/gpu/core/GPU/LBMKernel.cu
@@ -836,7 +836,7 @@ void CalcMacCompSP27RotatingToStatic(
     unsigned int numberOfThreads,
     real* DD,
     bool isEvenTimestep,
-    ParameterRotatingGridSimulation* parameterRotDevice)
+    ParameterRotatingGridHostDevice* parameterRotDevice)
 {
     vf::cuda::CudaGrid grid = vf::cuda::CudaGrid(numberOfThreads, numberOfLBnodes);
 
@@ -4097,7 +4097,7 @@ template void ScaleCF_compressible<false>(LBMSimulationParameter * parameterDevi
 void InterpolateStaticToRotating(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *baseToNested,
     ICellNeigh &neighborBaseToNested)
 {
@@ -4146,7 +4146,7 @@ void InterpolateStaticToRotating(
 void TraverseStaticToRotating(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *baseToNested,
     ICellNeigh &neighborBaseToNested)
 {
@@ -4182,7 +4182,7 @@ void TraverseStaticToRotating(
 void InterpolateRotatingToStatic(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *nestedToBase,
     ICellNeigh &neighborNestedToBase)
 {
@@ -4230,7 +4230,7 @@ void InterpolateRotatingToStatic(
 void TraverseRotatingToStatic(
     LBMSimulationParameter *parameterDeviceS,
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice,
+    ParameterRotatingGridHostDevice *paraRotDevice,
     ICells *nestedToBase,
     ICellNeigh &neighborNestedToBase)
 {
@@ -4265,7 +4265,7 @@ void TraverseRotatingToStatic(
 
 void UpdateGlobalCoordinates(
     LBMSimulationParameter *parameterDeviceR,
-    ParameterRotatingGridSimulation *paraRotDevice
+    ParameterRotatingGridHostDevice *paraRotDevice
 )
 {
     dim3 grid = vf::cuda::getCudaGrid(parameterDeviceR->numberofthreads, parameterDeviceR->numberOfNodes);
diff --git a/src/gpu/core/LBM/Simulation.cpp b/src/gpu/core/LBM/Simulation.cpp
index 1cf5dcf46..8b68f7bdb 100644
--- a/src/gpu/core/LBM/Simulation.cpp
+++ b/src/gpu/core/LBM/Simulation.cpp
@@ -449,6 +449,10 @@ void Simulation::allocNeighborsOffsetsScalesAndBoundaries(GridProvider &gridProv
 
 void Simulation::run()
 {
+    // #TODO: do not hard code level and move to appropriate place
+    if (para->getRotatingGridParameter() != nullptr)
+        updateGrid27->rotateGridInInitializationProcess(0);
+
     this->initTimers();
 
     //////////////////////////////////////////////////////////////////////////
diff --git a/src/gpu/core/Parameter/ParameterRotatingGrid.h b/src/gpu/core/Parameter/ParameterRotatingGrid.h
index c3484739c..e2ebd37ac 100644
--- a/src/gpu/core/Parameter/ParameterRotatingGrid.h
+++ b/src/gpu/core/Parameter/ParameterRotatingGrid.h
@@ -8,25 +8,29 @@
 #include <basics/DataTypes.h>
 #include <basics/geometry3d/Axis.h>
 
-struct ParameterRotatingGridSimulation {
+struct ParameterRotatingGridHostDevice
+{
     std::array<real, 3> centerPoint;
-    real *nestedCoordinatesX = nullptr; // in local coordinate system of rotating grid
-    real *nestedCoordinatesY = nullptr;
-    real *nestedCoordinatesZ = nullptr;
+    real* nestedCoordinatesX = nullptr; // in local coordinate system of rotating grid
+    real* nestedCoordinatesY = nullptr;
+    real* nestedCoordinatesZ = nullptr;
     uint sizeOfNestedCoordinates;
     uint memorySizeOfNestedCoordinates;
     std::array<real, 3> gridAngle = { 0.0, 0.0, 0.0 };
     std::array<real, 3> angularVelocity = { 0.0, 0.0, 0.0 };
 };
 
-struct ParameterRotatingGrid {
+class ParameterRotatingGrid
+{
 public:
-    ParameterRotatingGrid(const std::array<real, 3> &centerPoint, const Axis &rotationalAxis, uint sizeOfNestedCoordinates);
-    void fillNestedCoordinateVectorsOnHost(const std::array<real *, 3> &globalCoordinates);
+    ParameterRotatingGrid(const std::array<real, 3>& centerPoint, const Axis& rotationalAxis, uint sizeOfNestedCoordinates);
+    void fillNestedCoordinateVectorsOnHost(const std::array<real*, 3>& globalCoordinates);
 
     const Axis rotationalAxis;
-    SPtr<ParameterRotatingGridSimulation> parameterRotHost = std::make_shared<ParameterRotatingGridSimulation>();
-    SPtr<ParameterRotatingGridSimulation> parameterRotDevice = std::make_shared<ParameterRotatingGridSimulation>();
+    SPtr<ParameterRotatingGridHostDevice> parameterRotHost = std::make_shared<ParameterRotatingGridHostDevice>();
+    SPtr<ParameterRotatingGridHostDevice> parameterRotDevice = std::make_shared<ParameterRotatingGridHostDevice>();
+
+    real initialGridRotation = 0.0; // for debugging purposes
 };
 
 #endif
-- 
GitLab