From 35ac9fb780040adeea31f3c2c7cfc162fb5cb789 Mon Sep 17 00:00:00 2001
From: Anna Wellmann <a.wellmann@tu-braunschweig.de>
Date: Fri, 25 Aug 2023 15:33:56 +0000
Subject: [PATCH] Enable vtk output for rotating grid and refactor coordinate
transformations
---
apps/gpu/RotatingGrid/RotatingGrid.cpp | 40 ++++---
.../Calculation/UpdateGrid27.cpp | 46 +++++++
src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu | 51 ++++++++
src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h | 30 +++++
src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh | 62 ++++++++++
.../interpolateRotatingToStatic.cu | 69 ++++++++---
.../interpolateStaticToRotating.cu | 70 ++++++++---
src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu | 113 ++++++++++++++++++
.../K17/K17CompressibleNavierStokes_Device.cu | 6 +-
.../CoordinateTransformation.h | 51 ++++++++
src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp | 24 +++-
.../Parameter/ParameterRotatingGrid.h | 2 +-
12 files changed, 508 insertions(+), 56 deletions(-)
diff --git a/apps/gpu/RotatingGrid/RotatingGrid.cpp b/apps/gpu/RotatingGrid/RotatingGrid.cpp
index 99c58f2e4..ad3e6111f 100644
--- a/apps/gpu/RotatingGrid/RotatingGrid.cpp
+++ b/apps/gpu/RotatingGrid/RotatingGrid.cpp
@@ -91,8 +91,8 @@ int main()
const real velocityLB = 0.05; // LB units
const uint nx = 64;
- const uint timeStepOut = 1000;
- const uint timeStepEnd = 10000;
+ const uint timeStepOut = 1;
+ const uint timeStepEnd = 10;
const uint timeStepStartOutput = 0;
//////////////////////////////////////////////////////////////////////////
@@ -111,8 +111,8 @@ int main()
gridBuilder->addCoarseGrid(-0.5 * L, -0.5 * L, -0.5 * L, 2.0 * L, 0.5 * L, 0.5 * L, dx);
- if (rotOrInt == Rot) gridBuilder->addGridRotatingGrid(std::make_shared<Cylinder>(0.1, 0.1, 0.1, 0.25 * L, 1. * L, Axis::x));
- if (rotOrInt == Int) gridBuilder->addGrid(std::make_shared<Cylinder>(0.1, 0.1, 0.1, 0.25 * L, 0.8 * L, Axis::x), 1);
+ if (rotOrInt == Rot) gridBuilder->addGridRotatingGrid(std::make_shared<Cylinder>(0.2, 0.1, 0.1, 0.25 * L, 1. * 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);
@@ -158,8 +158,8 @@ int main()
gridBuilder->setNoSlipBoundaryCondition(SideType::PZ);
- // gridBuilder->setVelocityBoundaryCondition(SideType::MX, 0.0, 0.0, 0.0);
- gridBuilder->setVelocityBoundaryCondition(SideType::MX, velocityLB, 0.0, 0.0);
+ gridBuilder->setVelocityBoundaryCondition(SideType::MX, 0.0, 0.0, 0.0);
+ // gridBuilder->setVelocityBoundaryCondition(SideType::MX, velocityLB, 0.0, 0.0);
gridBuilder->setPressureBoundaryCondition(SideType::PX, 0.0);
BoundaryConditionFactory bcFactory;
@@ -173,14 +173,26 @@ int main()
// set initial condition
//////////////////////////////////////////////////////////////////////////
- // 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;
- // // if (coordX > -0.34 && coordY > 0.27 && coordZ > -0.1 && coordX < -0.32 && coordY < 0.3 && coordZ < 0.11) rho = 1e-5;
- // vx = 0.0;
- // vy = 0.0;
- // vz = 0.0;
- // });
+ auto setPressPoint = [](std::array<real, 3> coordPressPoint, real rhoPressPoint, real dx,
+ std::array<real, 3> coordinates) -> real {
+ auto isCoordinateInPressPoint = [&](uint dimension) -> bool {
+ return coordinates[dimension] > coordPressPoint[dimension] - dx &&
+ coordinates[dimension] < coordPressPoint[dimension] + dx;
+ };
+ if (isCoordinateInPressPoint(0) && isCoordinateInPressPoint(1) && isCoordinateInPressPoint(2))
+ return rhoPressPoint;
+ else
+ return (real)0.0;
+ };
+
+ 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;
+ vz = 0.0;
+ });
//////////////////////////////////////////////////////////////////////////
// set copy mesh to simulation
diff --git a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
index 0ac72a9c0..e6e948169 100644
--- a/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
+++ b/src/gpu/VirtualFluids_GPU/Calculation/UpdateGrid27.cpp
@@ -424,4 +424,50 @@ 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);
+
+ SPtr<ParameterRotatingGrid> paraRot = para->getRotatingGridParameter();
+
+ uint level = 0;
+ real finalRotation = (real)0.78539816339744830961566084581988;
+ real steps = 400;
+ paraRot->parameterRotHost->angularVelocity[0] = finalRotation / steps;
+ paraRot->parameterRotDevice->angularVelocity[0] = paraRot->parameterRotHost->angularVelocity[0];
+ VF_LOG_WARNING("angular velocity: {}",paraRot->parameterRotHost->angularVelocity[0] );
+ for (int i =0; i < steps; i++){
+ paraRot->parameterRotHost->gridAngle[0] += paraRot->parameterRotHost->angularVelocity[0];
+ paraRot->parameterRotHost->gridAngle[1] += paraRot->parameterRotHost->angularVelocity[1];
+ paraRot->parameterRotHost->gridAngle[2] += paraRot->parameterRotHost->angularVelocity[2];
+ paraRot->parameterRotDevice->gridAngle[0] += paraRot->parameterRotDevice->angularVelocity[0];
+ paraRot->parameterRotDevice->gridAngle[1] += paraRot->parameterRotDevice->angularVelocity[1];
+ paraRot->parameterRotDevice->gridAngle[2] += paraRot->parameterRotDevice->angularVelocity[2];
+
+ // base to nested
+ TraverseStaticToRotating(
+ para->getParD(level).get(),
+ para->getParD(level+1).get(),
+ para->getRotatingGridParameter()->parameterRotDevice.get(),
+ ¶->getParD(level)->coarseToFine,
+ para->getParD(level)->neighborCoarseToFine);
+
+ // nested to base
+ TraverseRotatingToStatic(
+ para->getParD(level).get(),
+ para->getParD(level+1).get(),
+ para->getRotatingGridParameter()->parameterRotDevice.get(),
+ ¶->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;
+
+ 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]);
}
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
index 8907e8467..63de33118 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
@@ -34,6 +34,8 @@
#include "lbm/constants/D3Q27.h"
#include "basics/constants/NumericConstants.h"
#include "lbm/MacroscopicQuantities.h"
+#include "LBM/GPUHelperFunctions/CoordinateTransformation.h"
+
#include "Kernel/Utilities/DistributionHelper.cuh"
@@ -309,7 +311,56 @@ __global__ void LBCalcMacCompSP27(
+__global__ void LBCalcMacCompSP27RotatingToStatic(
+ real *vxD,
+ real *vyD,
+ real *vzD,
+ real *rhoD,
+ real *pressD,
+ unsigned int *geoD,
+ unsigned int *neighborX,
+ unsigned int *neighborY,
+ unsigned int *neighborZ,
+ unsigned long long numberOfLBnodes,
+ real *distributions,
+ bool isEvenTimestep,
+ real angleX,
+ real angleY,
+ real angleZ
+ )
+{
+ ////////////////////////////////////////////////////////////////////////////////
+ //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
+ //!
+ const unsigned nodeIndex = getNodeIndex();
+ if(nodeIndex >= numberOfLBnodes)
+ return;
+
+ pressD[nodeIndex] = c0o1;
+ rhoD[nodeIndex] = c0o1;
+ vxD[nodeIndex] = c0o1;
+ vyD[nodeIndex] = c0o1;
+ vzD[nodeIndex] = c0o1;
+
+ if (!isValidFluidNode(geoD[nodeIndex]))
+ return;
+
+ DistributionWrapper distr_wrapper(distributions, numberOfLBnodes, isEvenTimestep, nodeIndex, neighborX, neighborY, neighborZ);
+ const auto &distribution = distr_wrapper.distribution;
+
+ rhoD[nodeIndex] = vf::lbm::getDensity(distribution.f);
+ real velocityRotatingX = vf::lbm::getCompressibleVelocityX1(distribution.f, rhoD[nodeIndex]);
+ real velocityRotatingY = vf::lbm::getCompressibleVelocityX2(distribution.f, rhoD[nodeIndex]);
+ real velocityRotatingZ = vf::lbm::getCompressibleVelocityX3(distribution.f, rhoD[nodeIndex]);
+ pressD[nodeIndex] = vf::lbm::getPressure(distribution.f, rhoD[nodeIndex], vxD[nodeIndex], vyD[nodeIndex], vzD[nodeIndex]);
+
+ rotateVelocityFromGlobalToRotating(velocityRotatingX, velocityRotatingY, velocityRotatingZ, angleX, angleY, angleZ);
+
+ vxD[nodeIndex] = velocityRotatingX;
+ vyD[nodeIndex] = velocityRotatingY;
+ vzD[nodeIndex] = velocityRotatingZ;
+}
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h b/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h
index e36a5c6e5..e1608160c 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h
+++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h
@@ -389,6 +389,22 @@ void CalcMacCompSP27(real* vxD,
real* DD,
bool isEvenTimestep);
+void CalcMacCompSP27RotatingToStatic(
+ real* vxD,
+ real* vyD,
+ real* vzD,
+ real* rhoD,
+ real* pressD,
+ unsigned int* geoD,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned long long numberOfLBnodes,
+ unsigned int numberOfThreads,
+ real* DD,
+ bool isEvenTimestep,
+ ParameterRotatingGridSimulation* parameterRotDevice);
+
void CalcMacThS7( real* Conc,
unsigned int* geoD,
unsigned int* neighborX,
@@ -1630,6 +1646,13 @@ void InterpolateStaticToRotating(
ICells *baseToNested,
ICellNeigh &neighborBaseToNested);
+void TraverseStaticToRotating(
+ LBMSimulationParameter *parameterDeviceS,
+ LBMSimulationParameter *parameterDeviceR,
+ ParameterRotatingGridSimulation *paraRotDevice,
+ ICells *baseToNested,
+ ICellNeigh &neighborBaseToNested);
+
void InterpolateRotatingToStatic(
LBMSimulationParameter *parameterDeviceS,
LBMSimulationParameter *parameterDeviceR,
@@ -1637,6 +1660,13 @@ void InterpolateRotatingToStatic(
ICells *nestedToBase,
ICellNeigh &neighborNestedToBase);
+void TraverseRotatingToStatic(
+ LBMSimulationParameter *parameterDeviceS,
+ LBMSimulationParameter *parameterDeviceR,
+ ParameterRotatingGridSimulation *paraRotDevice,
+ ICells *nestedToBase,
+ ICellNeigh &neighborNestedToBase);
+
void UpdateGlobalCoordinates(
LBMSimulationParameter *parameterDeviceR,
ParameterRotatingGridSimulation *paraRotDevice
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
index d9535d7af..941f83252 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
+++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
@@ -351,6 +351,24 @@ __global__ void LBCalcMacCompSP27( real* vxD,
real* DD,
bool isEvenTimestep);
+__global__ void LBCalcMacCompSP27RotatingToStatic(
+ real *vxD,
+ real *vyD,
+ real *vzD,
+ real *rhoD,
+ real *pressD,
+ unsigned int *geoD,
+ unsigned int *neighborX,
+ unsigned int *neighborY,
+ unsigned int *neighborZ,
+ unsigned long long numberOfLBnodes,
+ real *distributions,
+ bool isEvenTimestep,
+ real angleX,
+ real angleY,
+ real angleZ
+ );
+
__global__ void CalcConc7( real* Conc,
unsigned int* geoD,
unsigned int* neighborX,
@@ -1902,6 +1920,28 @@ __global__ void interpolateStaticToRotating(
bool isEvenTimestep,
real dx);
+__global__ void traverseStaticToRotating(
+ unsigned int numberOfInterfaceNodes,
+ unsigned int *indicesStaticCell,
+ const unsigned int *indicesRotating,
+ const real *coordDestinationX,
+ const real *coordDestinationY,
+ const real *coordDestinationZ,
+ const real *coordSourceX,
+ const real *coordSourceY,
+ const real *coordSourceZ,
+ const uint *neighborXstatic,
+ const uint *neighborYstatic,
+ const uint *neighborZstatic,
+ const uint *neighborMMMstatic,
+ real centerCoordX,
+ real centerCoordY,
+ real centerCoordZ,
+ real angleX,
+ real angleY,
+ real angleZ,
+ real dx);
+
__global__ void interpolateRotatingToStatic(
real *distributionsStatic,
real *distributionsRotating,
@@ -1936,6 +1976,28 @@ __global__ void interpolateRotatingToStatic(
bool isEvenTimestep,
real dx);
+__global__ void traverseRotatingToStatic(
+ unsigned int numberOfInterfaceNodes,
+ const unsigned int *indicesStatic,
+ unsigned int *indicesRotatingCell,
+ const real *coordDestinationX,
+ const real *coordDestinationY,
+ const real *coordDestinationZ,
+ const real *coordSourceX,
+ const real *coordSourceY,
+ const real *coordSourceZ,
+ const uint *neighborXrotating,
+ const uint *neighborYrotating,
+ const uint *neighborZrotating,
+ const uint *neighborMMMrotating,
+ real centerCoordX,
+ real centerCoordY,
+ real centerCoordZ,
+ real angleX,
+ real angleY,
+ real angleZ,
+ real dx);
+
__global__ void updateGlobalCoordinates(
unsigned int numberOfNodes,
real *globalX,
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
index bd870c837..54eae4389 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
@@ -42,6 +42,7 @@
#include <lbm/refinement/Coefficients.h>
#include <lbm/refinement/InterpolationCF.h>
+
using namespace vf::lbm;
__device__ __inline__ void calculateRotationalForces(real &forceX, real &forceY, real &forceZ, real angularVelocityX,
@@ -361,25 +362,11 @@ __global__ void interpolateRotatingToStatic(
vvx = coefficients.a000;
vvy = coefficients.b000;
vvz = coefficients.c000;
- vvxTemp = vvx;
- vvyTemp = vvy;
- vvzTemp = vvz;
- if (angleX != c0o1) {
- vvyTemp = vvy * cos(angleX) + vvz * sin(angleX);
- vvzTemp = -vvy * sin(angleX) + vvz * cos(angleX);
- } else if (angleY != c0o1) {
- // rotate in y
- vvxTemp = vvx * cos(angleY) - vvz * sin(angleY);
- vvzTemp = vvx * sin(angleY) + vvz * cos(angleY);
- } else if (angleZ != c0o1) {
- // rotate in z
- vvxTemp = vvx * cos(angleZ) + vvy * sin(angleZ);
- vvyTemp = -vvx * sin(angleZ) + vvy * cos(angleZ);
- }
- vvx = vvxTemp;
- vvy = vvyTemp;
- vvz = vvzTemp;
+ ////////////////////////////////////////////////////////////////////////////////
+ //! - rotate the velocities
+ //!
+ rotateVelocityFromGlobalToRotating(vvx, vvy, vvz, angleX, angleY, angleZ);
////////////////////////////////////////////////////////////////////////////////
// calculate the squares of the velocities
@@ -584,3 +571,49 @@ __global__ void updateGlobalCoordinates(
globalY[nodeIndex] = globalYtemp;
globalZ[nodeIndex] = globalZtemp;
}
+
+__global__ void traverseRotatingToStatic(
+ unsigned int numberOfInterfaceNodes,
+ const unsigned int *indicesStatic,
+ unsigned int *indicesRotatingCell,
+ const real *coordDestinationX,
+ const real *coordDestinationY,
+ const real *coordDestinationZ,
+ const real *coordSourceX,
+ const real *coordSourceY,
+ const real *coordSourceZ,
+ const uint *neighborXrotating,
+ const uint *neighborYrotating,
+ const uint *neighborZrotating,
+ const uint *neighborMMMrotating,
+ real centerCoordX,
+ real centerCoordY,
+ real centerCoordZ,
+ real angleX,
+ real angleY,
+ real angleZ,
+ real dx)
+{
+
+ const unsigned listIndex = vf::gpu::getNodeIndex();
+ if (listIndex >= numberOfInterfaceNodes) return;
+
+ const uint destinationIndex = indicesStatic[listIndex];
+ const uint previousSourceIndex = indicesRotatingCell[listIndex];
+ const uint indexNeighborMMMsource = neighborMMMrotating[previousSourceIndex];
+
+ real rotatedCoordDestinationX;
+ real rotatedCoordDestinationY;
+ real rotatedCoordDestinationZ;
+ transformGlobalToRotating(rotatedCoordDestinationX, rotatedCoordDestinationY, rotatedCoordDestinationZ,
+ coordDestinationX[destinationIndex], coordDestinationY[destinationIndex],
+ coordDestinationZ[destinationIndex], centerCoordX, centerCoordY, centerCoordZ, angleX, angleY,
+ angleZ);
+
+ const uint sourceIndex = traverseSourceCell(rotatedCoordDestinationX, rotatedCoordDestinationY, rotatedCoordDestinationZ,
+ indexNeighborMMMsource, coordSourceX[indexNeighborMMMsource],
+ coordSourceY[indexNeighborMMMsource], coordSourceZ[indexNeighborMMMsource],
+ neighborXrotating, neighborYrotating, neighborZrotating, dx);
+
+ indicesRotatingCell[listIndex] = sourceIndex;
+}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
index 7707a1515..0119aab6e 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
@@ -238,25 +238,7 @@ __global__ void interpolateStaticToRotating(
////////////////////////////////////////////////////////////////////////////////
//! - rotate the velocities
//!
-
- real vvxTemp = vvx;
- real vvyTemp = vvy;
- real vvzTemp = vvz;
- if (angleX != c0o1) {
- vvyTemp = vvy * cos(angleX) - vvz * sin(angleX);
- vvzTemp = vvy * sin(angleX) + vvz * cos(angleX);
- } else if (angleY != c0o1) {
- // rotate in y
- vvxTemp = vvx * cos(angleY) + vvz * sin(angleY);
- vvzTemp = -vvx * sin(angleY) + vvz * cos(angleY);
- } else if (angleZ != c0o1) {
- // rotate in z
- vvxTemp = vvx * cos(angleZ) - vvy * sin(angleZ);
- vvyTemp = vvx * sin(angleZ) + vvy * cos(angleZ);
- }
- vvx = vvxTemp;
- vvy = vvyTemp;
- vvz = vvzTemp;
+ rotateVelocityFromRotatingToGlobal(vvx, vvy, vvz, angleX, angleY, angleZ);
////////////////////////////////////////////////////////////////////////////////
// calculate the squares of the velocities
@@ -427,3 +409,53 @@ __global__ void interpolateStaticToRotating(
vf::gpu::ListIndices indicesRotatingForWriting(destinationIndex, neighborXrotating, neighborYrotating, neighborZrotating);
vf::gpu::write(distRoating, indicesRotatingForWriting, fRotating);
}
+
+__global__ void traverseStaticToRotating(
+ unsigned int numberOfInterfaceNodes,
+ unsigned int *indicesStaticCell,
+ const unsigned int *indicesRotating,
+ const real *coordDestinationX,
+ const real *coordDestinationY,
+ const real *coordDestinationZ,
+ const real *coordSourceX,
+ const real *coordSourceY,
+ const real *coordSourceZ,
+ const uint *neighborXstatic,
+ const uint *neighborYstatic,
+ const uint *neighborZstatic,
+ const uint *neighborMMMstatic,
+ real centerCoordX,
+ real centerCoordY,
+ real centerCoordZ,
+ real angleX,
+ real angleY,
+ real angleZ,
+ real dx)
+{
+ // Interpolate from a cell on the static grid (source cell) to a node on the rotating grid (destination cell)
+
+ // 1. calculate the indices of involved nodes
+ const unsigned listIndex = vf::gpu::getNodeIndex();
+ if (listIndex >= numberOfInterfaceNodes) return;
+ const uint destinationIndex = indicesRotating[listIndex];
+ const uint previousSourceIndex = indicesStaticCell[listIndex];
+ const uint indexNeighborMMMsource = neighborMMMstatic[previousSourceIndex];
+
+ // calculate the coordinates of the destination cell in the global coordinate system
+ real globalCoordDestinationX;
+ real globalCoordDestinationY;
+ real globalCoordDestinationZ;
+ transformRotatingToGlobal(globalCoordDestinationX, globalCoordDestinationY, globalCoordDestinationZ,
+ coordDestinationX[destinationIndex], coordDestinationY[destinationIndex],
+ coordDestinationZ[destinationIndex], centerCoordX, centerCoordY, centerCoordZ, angleX, angleY,
+ angleZ);
+
+ // find the new index of the source cell (a static cell) after the rotation
+ const uint sourceIndex =
+ traverseSourceCell(globalCoordDestinationX, globalCoordDestinationY, globalCoordDestinationZ, indexNeighborMMMsource,
+ coordSourceX[indexNeighborMMMsource], coordSourceY[indexNeighborMMMsource],
+ coordSourceZ[indexNeighborMMMsource], neighborXstatic, neighborYstatic, neighborZstatic, dx);
+
+ // write the new source index to the array
+ indicesStaticCell[listIndex] = sourceIndex;
+}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
index a0049db48..3dfb649ef 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -820,6 +820,46 @@ void CalcMacCompSP27(
isEvenTimestep);
getLastCudaError("LBCalcMacCompSP27 execution failed");
}
+
+
+void CalcMacCompSP27RotatingToStatic(
+ real* vxD,
+ real* vyD,
+ real* vzD,
+ real* rhoD,
+ real* pressD,
+ unsigned int* geoD,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned long long numberOfLBnodes,
+ unsigned int numberOfThreads,
+ real* DD,
+ bool isEvenTimestep,
+ ParameterRotatingGridSimulation* parameterRotDevice)
+{
+ vf::cuda::CudaGrid grid = vf::cuda::CudaGrid(numberOfThreads, numberOfLBnodes);
+
+ LBCalcMacCompSP27RotatingToStatic<<< grid.grid, grid.threads >>> (
+ vxD,
+ vyD,
+ vzD,
+ rhoD,
+ pressD,
+ geoD,
+ neighborX,
+ neighborY,
+ neighborZ,
+ numberOfLBnodes,
+ DD,
+ isEvenTimestep,
+ parameterRotDevice->gridAngle[0],
+ parameterRotDevice->gridAngle[1],
+ parameterRotDevice->gridAngle[2]
+ );
+ getLastCudaError("LBCalcMacCompSP27RotatingToStatic execution failed");
+}
+
//////////////////////////////////////////////////////////////////////////
void CalcMacThS7(
real* Conc,
@@ -4102,6 +4142,43 @@ void InterpolateStaticToRotating(
getLastCudaError("interpolateStaticToRotating execution failed");
}
+
+void TraverseStaticToRotating(
+ LBMSimulationParameter *parameterDeviceS,
+ LBMSimulationParameter *parameterDeviceR,
+ ParameterRotatingGridSimulation *paraRotDevice,
+ ICells *baseToNested,
+ ICellNeigh &neighborBaseToNested)
+{
+ dim3 grid = vf::cuda::getCudaGrid(parameterDeviceS->numberofthreads, baseToNested->numberOfCells);
+ dim3 threads(parameterDeviceS->numberofthreads, 1, 1);
+
+ traverseStaticToRotating<<<grid, threads, 0, CU_STREAM_LEGACY>>>(
+ baseToNested->numberOfCells,
+ baseToNested->coarseCellIndices,
+ baseToNested->fineCellIndices,
+ paraRotDevice->nestedCoordinatesX,
+ paraRotDevice->nestedCoordinatesY,
+ paraRotDevice->nestedCoordinatesZ,
+ parameterDeviceS->coordinateX,
+ parameterDeviceS->coordinateY,
+ parameterDeviceS->coordinateZ,
+ parameterDeviceS->neighborX,
+ parameterDeviceS->neighborY,
+ parameterDeviceS->neighborZ,
+ parameterDeviceS->neighborInverse,
+ paraRotDevice->centerPoint[0],
+ paraRotDevice->centerPoint[1],
+ paraRotDevice->centerPoint[2],
+ paraRotDevice->gridAngle[0],
+ paraRotDevice->gridAngle[1],
+ paraRotDevice->gridAngle[2],
+ parameterDeviceS->gridSpacing
+ );
+
+ getLastCudaError("traverseStaticToRotating execution failed");
+}
+
void InterpolateRotatingToStatic(
LBMSimulationParameter *parameterDeviceS,
LBMSimulationParameter *parameterDeviceR,
@@ -4150,6 +4227,42 @@ void InterpolateRotatingToStatic(
getLastCudaError("interpolateRotatingToStatic execution failed");
}
+void TraverseRotatingToStatic(
+ LBMSimulationParameter *parameterDeviceS,
+ LBMSimulationParameter *parameterDeviceR,
+ ParameterRotatingGridSimulation *paraRotDevice,
+ ICells *nestedToBase,
+ ICellNeigh &neighborNestedToBase)
+{
+ dim3 grid = vf::cuda::getCudaGrid(parameterDeviceS->numberofthreads, nestedToBase->numberOfCells);
+ dim3 threads(parameterDeviceS->numberofthreads, 1, 1);
+
+ traverseRotatingToStatic<<<grid, threads, 0, CU_STREAM_LEGACY>>>(
+ nestedToBase->numberOfCells,
+ nestedToBase->coarseCellIndices,
+ nestedToBase->fineCellIndices,
+ parameterDeviceS->coordinateX,
+ parameterDeviceS->coordinateY,
+ parameterDeviceS->coordinateZ,
+ paraRotDevice->nestedCoordinatesX,
+ paraRotDevice->nestedCoordinatesY,
+ paraRotDevice->nestedCoordinatesZ,
+ parameterDeviceR->neighborX,
+ parameterDeviceR->neighborY,
+ parameterDeviceR->neighborZ,
+ parameterDeviceR->neighborInverse,
+ paraRotDevice->centerPoint[0],
+ paraRotDevice->centerPoint[1],
+ paraRotDevice->centerPoint[2],
+ paraRotDevice->gridAngle[0],
+ paraRotDevice->gridAngle[1],
+ paraRotDevice->gridAngle[2],
+ parameterDeviceS->gridSpacing
+ );
+
+ getLastCudaError("traverseRotatingToStatic execution failed");
+}
+
void UpdateGlobalCoordinates(
LBMSimulationParameter *parameterDeviceR,
ParameterRotatingGridSimulation *paraRotDevice
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cu b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cu
index 732e96ce4..81e2579eb 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cu
@@ -209,9 +209,9 @@ __global__ void K17CompressibleNavierStokes_Device(
//! DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
//!
real factor = c1o1;
- for (size_t i = 1; i <= level; i++) {
- factor *= c2o1;
- }
+ // for (size_t i = 1; i <= level; i++) {
+ // factor *= c2o1;
+ // }
real fx = forces[0];
real fy = forces[1];
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h
index 14fc5794e..fea772f62 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h
@@ -3,6 +3,33 @@
#include <basics/DataTypes.h>
#include <math.h>
+#include <basics/constants/NumericConstants.h>
+
+using namespace vf::basics::constant;
+
+__inline__ __device__ void rotateVelocityFromRotatingToGlobal(real &velocityX, real &velocityY, real &velocityZ, real angleX,
+ real angleY, real angleZ)
+{
+ real velocityXTemp = velocityX;
+ real velocityYTemp = velocityY;
+ real velocityZTemp = velocityZ;
+ if (angleX != c0o1) {
+ velocityYTemp = velocityY * cos(angleX) - velocityZ * sin(angleX);
+ velocityZTemp = velocityY * sin(angleX) + velocityZ * cos(angleX);
+ } else if (angleY != c0o1) {
+ // rotate in y
+ velocityXTemp = velocityX * cos(angleY) + velocityZ * sin(angleY);
+ velocityZTemp = -velocityX * sin(angleY) + velocityZ * cos(angleY);
+ } else if (angleZ != c0o1) {
+ // rotate in z
+ velocityXTemp = velocityX * cos(angleZ) - velocityY * sin(angleZ);
+ velocityYTemp = velocityX * sin(angleZ) + velocityY * cos(angleZ);
+ }
+ velocityX = velocityXTemp;
+ velocityY = velocityYTemp;
+ velocityZ = velocityZTemp;
+}
+
__inline__ __device__ void transformRotatingToGlobal(real &globalX, real &globalY, real &globalZ, real localX, real localY,
real localZ, real centerCoordX, real centerCoordY, real centerCoordZ,
real angleX, real angleY, real angleZ)
@@ -32,6 +59,30 @@ __inline__ __device__ void transformRotatingToGlobal(real &globalX, real &global
globalZ += centerCoordZ;
}
+
+__inline__ __device__ void rotateVelocityFromGlobalToRotating(real &velocityX, real &velocityY, real &velocityZ, real angleX,
+ real angleY, real angleZ)
+{
+ real velocityXTemp = velocityX;
+ real velocityYTemp = velocityY;
+ real velocityZTemp = velocityZ;
+ if (angleX != c0o1) {
+ velocityYTemp = velocityY * cos(angleX) + velocityZ * sin(angleX);
+ velocityZTemp = -velocityY * sin(angleX) + velocityZ * cos(angleX);
+ } else if (angleY != c0o1) {
+ // rotate in y
+ velocityXTemp = velocityX * cos(angleY) - velocityZ * sin(angleY);
+ velocityZTemp = velocityX * sin(angleY) + velocityZ * cos(angleY);
+ } else if (angleZ != c0o1) {
+ // rotate in z
+ velocityXTemp = velocityX * cos(angleZ) + velocityY * sin(angleZ);
+ velocityYTemp = -velocityX * sin(angleZ) + velocityY * cos(angleZ);
+ }
+ velocityX = velocityXTemp;
+ velocityY = velocityYTemp;
+ velocityZ = velocityZTemp;
+}
+
__inline__ __device__ void transformGlobalToRotating(real &rotatingX, real &rotatingY, real &rotatingZ, real globalX,
real globalY, real globalZ, real centerCoordX, real centerCoordY,
real centerCoordZ, real angleX, real angleY, real angleZ)
diff --git a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
index 948f923c1..65999b0ae 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -768,6 +768,8 @@ void Simulation::readAndWriteFiles(uint timestep)
// para->getParD(lev)->d0SP.f[0],
// para->getParD(lev)->evenOrOdd);
// getLastCudaError("CalcMacSP27 execution failed");
+
+ if (lev == 0){
CalcMacCompSP27(para->getParD(lev)->velocityX,
para->getParD(lev)->velocityY,
para->getParD(lev)->velocityZ,
@@ -782,6 +784,26 @@ void Simulation::readAndWriteFiles(uint timestep)
para->getParD(lev)->distributions.f[0],
para->getParD(lev)->isEvenTimestep);
getLastCudaError("CalcMacSP27 execution failed");
+ }
+ if (lev == 1)
+ {
+ CalcMacCompSP27RotatingToStatic(para->getParD(lev)->velocityX,
+ para->getParD(lev)->velocityY,
+ para->getParD(lev)->velocityZ,
+ para->getParD(lev)->rho,
+ para->getParD(lev)->pressure,
+ para->getParD(lev)->typeOfGridNode,
+ para->getParD(lev)->neighborX,
+ para->getParD(lev)->neighborY,
+ para->getParD(lev)->neighborZ,
+ para->getParD(lev)->numberOfNodes,
+ para->getParD(lev)->numberofthreads,
+ para->getParD(lev)->distributions.f[0],
+ para->getParD(lev)->isEvenTimestep,
+ para->getRotatingGridParameter()->parameterRotDevice.get());
+ getLastCudaError("CalcMacSP27 execution failed");
+ }
+
// // overwrite with wall nodes
// SetOutputWallVelocitySP27( para->getParD(lev)->numberofthreads,
// para->getParD(lev)->velocityX,
@@ -971,7 +993,7 @@ void Simulation::readAndWriteFiles(uint timestep)
UpdateGlobalCoordinates(para->getParD(1).get(), para->getRotatingGridParameter()->parameterRotDevice.get());
cudaMemoryManager->cudaCopyCoordDeviceToHost(1);
- // cudaMemoryManager->cudaCopyCoordRotationDeviceToHost(1);
+
// cudaMemoryManager->cudaCopyInterfaceCFDeviceToHost(0);
// cudaMemoryManager->cudaCopyInterfaceFCDeviceToHost(0);
// InterfaceDebugWriter::writeInterfaceLinesDebugCF(para.get(), timestep);
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/ParameterRotatingGrid.h b/src/gpu/VirtualFluids_GPU/Parameter/ParameterRotatingGrid.h
index ac78941f8..c3484739c 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/ParameterRotatingGrid.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/ParameterRotatingGrid.h
@@ -16,7 +16,7 @@ struct ParameterRotatingGridSimulation {
uint sizeOfNestedCoordinates;
uint memorySizeOfNestedCoordinates;
std::array<real, 3> gridAngle = { 0.0, 0.0, 0.0 };
- std::array<real, 3> angularVelocity = { 0.02, 0.0, 0.0 };
+ std::array<real, 3> angularVelocity = { 0.0, 0.0, 0.0 };
};
struct ParameterRotatingGrid {
--
GitLab