From 43ff5b02d1de112627dc6a3f8a5c9ca167cda6b1 Mon Sep 17 00:00:00 2001
From: Anna Wellmann <a.wellmann@tu-braunschweig.de>
Date: Thu, 24 Aug 2023 09:44:56 +0000
Subject: [PATCH] First try of interpolation rotating to static
---
src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh | 3 +
.../interpolateRotatingToStatic.cu | 264 +++++++++++++++++-
.../interpolateStaticToRotating.cu | 16 +-
src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu | 3 +
.../CoordinateTransformation.h | 1 +
5 files changed, 279 insertions(+), 8 deletions(-)
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
index 98d15c770..d9535d7af 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
+++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
@@ -1920,6 +1920,9 @@ __global__ void interpolateRotatingToStatic(
const uint *neighborYrotating,
const uint *neighborZrotating,
const uint *neighborMMMrotating,
+ const uint *neighborXstatic,
+ const uint *neighborYstatic,
+ const uint *neighborZstatic,
real centerCoordX,
real centerCoordY,
real centerCoordZ,
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
index c4e66a9ee..4dbbb1619 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
@@ -43,7 +43,6 @@
#include <lbm/refinement/InterpolationCF.h>
using namespace vf::lbm;
-using namespace vf::gpu;
__device__ __inline__ void calculateRotationalForces(real &forceX, real &forceY, real &forceZ, real angularVelocityX,
real angularVelocityY, real angularVelocityZ, real coordSourceXlocal,
@@ -159,6 +158,9 @@ __global__ void interpolateRotatingToStatic(
const uint *neighborYrotating,
const uint *neighborZrotating,
const uint *neighborMMMrotating,
+ const uint *neighborXstatic,
+ const uint *neighborYstatic,
+ const uint *neighborZstatic,
real centerCoordX,
real centerCoordY,
real centerCoordZ,
@@ -276,7 +278,7 @@ __global__ void interpolateRotatingToStatic(
isEvenTimestep, forcesVertexX, forcesVertexY, forcesVertexZ);
// 2. calculate the coefficients for the interpolation
- // For this the relative coordinates of the destination pooint inside the source cell are needed
+ // For this the relative coordinates of the destination point inside the source cell are needed
// this cell coordinate system is centered at the middle of the source cell. The source cell has a side lenght of one
// add tangential velocity to velocity for coefficient computaion
const real cellCoordDestinationX = (coordSourceX[sourceIndex] + rotatedCoordDestinationX) / dx - c1o2;
@@ -286,7 +288,265 @@ __global__ void interpolateRotatingToStatic(
vf::lbm::Coefficients coefficients;
momentsSet.calculateCoefficients(coefficients, cellCoordDestinationX, cellCoordDestinationY, cellCoordDestinationZ);
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Set all moments for the destination cell to zero
+ //!
+ real m111 = c0o1;
+ real m211 = c0o1;
+ real m011 = c0o1;
+ real m121 = c0o1;
+ real m101 = c0o1;
+ real m112 = c0o1;
+ real m110 = c0o1;
+ real m221 = c0o1;
+ real m001 = c0o1;
+ real m201 = c0o1;
+ real m021 = c0o1;
+ real m212 = c0o1;
+ real m010 = c0o1;
+ real m210 = c0o1;
+ real m012 = c0o1;
+ real m122 = c0o1;
+ real m100 = c0o1;
+ real m120 = c0o1;
+ real m102 = c0o1;
+ real m222 = c0o1;
+ real m022 = c0o1;
+ real m202 = c0o1;
+ real m002 = c0o1;
+ real m220 = c0o1;
+ real m020 = c0o1;
+ real m200 = c0o1;
+ real m000 = c0o1;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Define aliases to use the same variable for the distributions (f's):
+ //!
+ real& f000 = m111;
+ real& fP00 = m211;
+ real& fM00 = m011;
+ real& f0P0 = m121;
+ real& f0M0 = m101;
+ real& f00P = m112;
+ real& f00M = m110;
+ real& fPP0 = m221;
+ real& fMM0 = m001;
+ real& fPM0 = m201;
+ real& fMP0 = m021;
+ real& fP0P = m212;
+ real& fM0M = m010;
+ real& fP0M = m210;
+ real& fM0P = m012;
+ real& f0PP = m122;
+ real& f0MM = m100;
+ real& f0PM = m120;
+ real& f0MP = m102;
+ real& fPPP = m222;
+ real& fMPP = m022;
+ real& fPMP = m202;
+ real& fMMP = m002;
+ real& fPPM = m220;
+ real& fMPM = m020;
+ real& fPMM = m200;
+ real& fMMM = m000;
+
+ ////////////////////////////////////////////////////////////////////////////////
+ //! - Set macroscopic values on destination node (zeroth and first order moments)
+ //!
+ real vvx, vvy, vvz, vvxTemp, vvyTemp, vvzTemp;
+
+ m000 = coefficients.d000; // m000 is press, if drho is interpolated directly
+ vvx = coefficients.a000;
+ vvy = coefficients.b000;
+ vvz = coefficients.c000;
+ vvxTemp = vvx;
+ vvyTemp = vvy;
+ vvzTemp = vvz;
+
+ if (angleX != 0) {
+ vvyTemp = vvy * cos(angleX) + vvz * sin(angleX);
+ vvzTemp = -vvy * sin(angleX) + vvz * cos(angleX);
+ } else if (angleY != 0) {
+ // rotate in y
+ vvxTemp = vvx * cos(angleY) - vvz * sin(angleY);
+ vvzTemp = vvx * sin(angleY) + vvz * cos(angleY);
+ } else if (angleZ != 0) {
+ // rotate in z
+ vvxTemp = vvx * cos(angleZ) + vvy * sin(angleZ);
+ vvyTemp = -vvx * sin(angleZ) + vvy * cos(angleZ);
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////
+ // calculate the squares of the velocities
+ //
+ real vxsq = vvx * vvx;
+ real vysq = vvy * vvy;
+ real vzsq = vvz * vvz;
+
+ ////////////////////////////////////////////////////////////////////////////////
+ //! - Set moments (second to sixth order) on destination node
+ //!
+ // linear combinations for second order moments
+ // The second order moments have to be rotated. First they are computed in a static frame of reference (SFOR)
+
+ real mxxPyyPzz = m000;
+
+ real mxxMyySFOR = -c2o3 * ((coefficients.a100 - coefficients.b010)) / omegaRotating * (c1o1 + m000);
+ real mxxMzzSFOR = -c2o3 * ((coefficients.a100 - coefficients.c001)) / omegaRotating * (c1o1 + m000);
+
+ real m011SFOR = -c1o3 * ((coefficients.b001 + coefficients.c010)) / omegaRotating * (c1o1 + m000);
+ real m101SFOR = -c1o3 * ((coefficients.a001 + coefficients.c100)) / omegaRotating * (c1o1 + m000);
+ real m110SFOR = -c1o3 * ((coefficients.a010 + coefficients.b100)) / omegaRotating * (c1o1 + m000);
+
+ // rotate some second order moments
+ real mxxMyy = mxxMyySFOR;
+ real mxxMzz = mxxMzzSFOR;
+ m011 = m011SFOR;
+ m101 = m101SFOR;
+ m110 = m110SFOR;
+ if (angleX != 0) {
+ mxxMyy = (mxxMyySFOR + mxxMzzSFOR + (mxxMyySFOR - mxxMzzSFOR) * cos(angleX * c2o1) -
+ c2o1 * m011SFOR * sin(angleX * c2o1)) *
+ c1o2;
+ mxxMzz = (mxxMyySFOR + mxxMzzSFOR + (-mxxMyySFOR + mxxMzzSFOR) * cos(angleX * c2o1) +
+ c2o1 * m011SFOR * sin(angleX * c2o1)) *
+ c1o2;
+
+ m011 = m011SFOR * cos(angleX * c2o1) + (mxxMyySFOR - mxxMzzSFOR) * cos(angleX) * sin(angleX);
+ m101 = m101SFOR * cos(angleX) - m110SFOR * sin(angleX);
+ m110 = m110SFOR * cos(angleX) + m101SFOR * sin(angleX);
+ } else if (angleY != 0) {
+ mxxMyy = mxxMyySFOR - mxxMzzSFOR * c1o2 + c1o2 * mxxMzzSFOR * cos(angleY * c2o1) + m101SFOR * sin(angleY * c2o1);
+ mxxMzz = mxxMzzSFOR * cos(angleY * c2o1) + c2o1 * m101SFOR * sin(angleY * c2o1);
+
+ m011 = m011SFOR * cos(angleY) - m110SFOR * sin(angleY);
+ m101 = m101SFOR * cos(angleY * c2o1) - mxxMzzSFOR * cos(angleY) * sin(angleY);
+ m110 = m110SFOR * cos(angleY) + m011SFOR * sin(angleY);
+ } else if (angleZ != 0) {
+ mxxMyy = mxxMyySFOR * cos(angleY * c2o1) - c2o1 * m110SFOR * sin(angleY * c2o1);
+ mxxMzz = -c1o2 * mxxMyySFOR + mxxMzzSFOR + c1o2 * mxxMyySFOR * cos(angleY * c2o1) - m110SFOR * sin(angleY * c2o1);
+
+ m011 = m011SFOR * cos(angleZ) + m101SFOR * sin(angleZ);
+ m101 = m101SFOR * cos(angleZ) - m011SFOR * sin(angleZ);
+ m110 = m110SFOR * cos(angleZ * c2o1) + mxxMyySFOR * cos(angleZ) * sin(angleZ);
+ }
+
+ // calculate remaining second order moments from previously rotated moments
+ real useNEQ = c1o1; // zero; //one; //.... one = on ..... zero = off
+ m200 = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz) * useNEQ;
+ m020 = c1o3 * (-c2o1 * mxxMyy + mxxMzz + mxxPyyPzz) * useNEQ;
+ m002 = c1o3 * (mxxMyy - c2o1 * mxxMzz + mxxPyyPzz) * useNEQ;
+
+ // linear combinations for third order moments
+ real mxxyPyzz, mxxyMyzz, mxxzPyyz, mxxzMyyz, mxyyPxzz, mxyyMxzz;
+ m111 = c0o1;
+
+ mxxyPyzz = c0o1;
+ mxxyMyzz = c0o1;
+ mxxzPyyz = c0o1;
+ mxxzMyyz = c0o1;
+ mxyyPxzz = c0o1;
+ mxyyMxzz = c0o1;
+
+ m210 = ( mxxyMyzz + mxxyPyzz) * c1o2;
+ m012 = (-mxxyMyzz + mxxyPyzz) * c1o2;
+ m201 = ( mxxzMyyz + mxxzPyyz) * c1o2;
+ m021 = (-mxxzMyyz + mxxzPyyz) * c1o2;
+ m120 = ( mxyyMxzz + mxyyPxzz) * c1o2;
+ m102 = (-mxyyMxzz + mxyyPxzz) * c1o2;
+
+ // fourth order moments
+ m022 = m000 * c1o9;
+ m202 = m022;
+ m220 = m022;
+
+ // fifth order moments
+
+ // sixth order moments
+ m222 = m000 * c1o27;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from central moments to well conditioned distributions as defined in Appendix J in
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
+ //! DOI:10.1016/j.camwa.2015.05.001 ]</b></a> see also Eq. (88)-(96) in <a
+ //! href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040
+ //! ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ backwardInverseChimeraWithK(m000, m100, m200, vvx, vxsq, c1o1, c1o1);
+ backwardChimera( m010, m110, m210, vvx, vxsq);
+ backwardInverseChimeraWithK(m020, m120, m220, vvx, vxsq, c3o1, c1o3);
+ backwardChimera( m001, m101, m201, vvx, vxsq);
+ backwardChimera( m011, m111, m211, vvx, vxsq);
+ backwardChimera( m021, m121, m221, vvx, vxsq);
+ backwardInverseChimeraWithK(m002, m102, m202, vvx, vxsq, c3o1, c1o3);
+ backwardChimera( m012, m112, m212, vvx, vxsq);
+ backwardInverseChimeraWithK(m022, m122, m222, vvx, vxsq, c9o1, c1o9);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ backwardInverseChimeraWithK(m000, m010, m020, vvy, vysq, c6o1, c1o6);
+ backwardChimera( m001, m011, m021, vvy, vysq);
+ backwardInverseChimeraWithK(m002, m012, m022, vvy, vysq, c18o1, c1o18);
+ backwardInverseChimeraWithK(m100, m110, m120, vvy, vysq, c3o2, c2o3);
+ backwardChimera( m101, m111, m121, vvy, vysq);
+ backwardInverseChimeraWithK(m102, m112, m122, vvy, vysq, c9o2, c2o9);
+ backwardInverseChimeraWithK(m200, m210, m220, vvy, vysq, c6o1, c1o6);
+ backwardChimera( m201, m211, m221, vvy, vysq);
+ backwardInverseChimeraWithK(m202, m212, m222, vvy, vysq, c18o1, c1o18);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ backwardInverseChimeraWithK(m000, m001, m002, vvz, vzsq, c36o1, c1o36);
+ backwardInverseChimeraWithK(m010, m011, m012, vvz, vzsq, c9o1, c1o9);
+ backwardInverseChimeraWithK(m020, m021, m022, vvz, vzsq, c36o1, c1o36);
+ backwardInverseChimeraWithK(m100, m101, m102, vvz, vzsq, c9o1, c1o9);
+ backwardInverseChimeraWithK(m110, m111, m112, vvz, vzsq, c9o4, c4o9);
+ backwardInverseChimeraWithK(m120, m121, m122, vvz, vzsq, c9o1, c1o9);
+ backwardInverseChimeraWithK(m200, m201, m202, vvz, vzsq, c36o1, c1o36);
+ backwardInverseChimeraWithK(m210, m211, m212, vvz, vzsq, c9o1, c1o9);
+ backwardInverseChimeraWithK(m220, m221, m222, vvz, vzsq, c36o1, c1o36);
+
+
+ // ////////////////////////////////////////////////////////////////////////////////
+
+ real fStatic[27];
+ fStatic[DIR_000] = f000;
+ fStatic[DIR_P00] = fP00;
+ fStatic[DIR_M00] = fM00;
+ fStatic[DIR_0P0] = f0P0;
+ fStatic[DIR_0M0] = f0M0;
+ fStatic[DIR_00P] = f00P;
+ fStatic[DIR_00M] = f00M;
+ fStatic[DIR_PP0] = fPP0;
+ fStatic[DIR_MM0] = fMM0;
+ fStatic[DIR_PM0] = fPM0;
+ fStatic[DIR_MP0] = fMP0;
+ fStatic[DIR_P0P] = fP0P;
+ fStatic[DIR_M0M] = fM0M;
+ fStatic[DIR_P0M] = fP0M;
+ fStatic[DIR_M0P] = fM0P;
+ fStatic[DIR_0PP] = f0PP;
+ fStatic[DIR_0MM] = f0MM;
+ fStatic[DIR_0PM] = f0PM;
+ fStatic[DIR_0MP] = f0MP;
+ fStatic[DIR_PPP] = fPPP;
+ fStatic[DIR_MPP] = fMPP;
+ fStatic[DIR_PMP] = fPMP;
+ fStatic[DIR_MMP] = fMMP;
+ fStatic[DIR_PPM] = fPPM;
+ fStatic[DIR_MPM] = fMPM;
+ fStatic[DIR_PMM] = fPMM;
+ fStatic[DIR_MMM] = fMMM;
+
+ // get distribution pointers for destination node
+ Distributions27 distStatic;
+ vf::gpu::getPointersToDistributions(distStatic, distributionsStatic, numberOfLBNodesStatic, isEvenTimestep);
+
+ // write
+ vf::gpu::ListIndices writeIndicesStatic(destinationIndex, neighborXstatic, neighborYstatic, neighborZstatic);
+ vf::gpu::write(distStatic, writeIndicesStatic, fStatic);
}
__global__ void updateGlobalCoordinates(
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
index ee96e6b8b..70dd6ff82 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
@@ -239,18 +239,22 @@ __global__ void interpolateStaticToRotating(
//! - rotate the velocities
//!
+ real vvxTemp, vvyTemp, vvzTemp;
if (angleX != 0) {
- vvy = vvy * cos(angleX) - vvz * sin(angleX);
- vvz = vvy * sin(angleX) + vvz * cos(angleX);
+ vvyTemp = vvy * cos(angleX) - vvz * sin(angleX);
+ vvzTemp = vvy * sin(angleX) + vvz * cos(angleX);
} else if (angleY != 0) {
// rotate in y
- vvx = vvx * cos(angleY) + vvz * sin(angleY);
- vvz = -vvx * sin(angleY) + vvz * cos(angleY);
+ vvxTemp = vvx * cos(angleY) + vvz * sin(angleY);
+ vvzTemp = -vvx * sin(angleY) + vvz * cos(angleY);
} else if (angleZ != 0) {
// rotate in z
- vvx = vvx * cos(angleZ) - vvy * sin(angleZ);
- vvy = vvx * sin(angleZ) + vvy * cos(angleZ);
+ vvxTemp = vvx * cos(angleZ) - vvy * sin(angleZ);
+ vvyTemp = vvx * sin(angleZ) + vvy * cos(angleZ);
}
+ vvx = vvxTemp;
+ vvy = vvyTemp;
+ vvz = vvzTemp;
////////////////////////////////////////////////////////////////////////////////
// calculate the squares of the velocities
diff --git a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
index 54bad32e7..a0049db48 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -4130,6 +4130,9 @@ void InterpolateRotatingToStatic(
parameterDeviceR->neighborY,
parameterDeviceR->neighborZ,
parameterDeviceR->neighborInverse,
+ parameterDeviceS->neighborX,
+ parameterDeviceS->neighborY,
+ parameterDeviceS->neighborZ,
paraRotDevice->centerPoint[0],
paraRotDevice->centerPoint[1],
paraRotDevice->centerPoint[2],
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h
index 9a7d55842..14fc5794e 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/CoordinateTransformation.h
@@ -45,6 +45,7 @@ __inline__ __device__ void transformGlobalToRotating(real &rotatingX, real &rota
rotatingX = globalX;
rotatingY = globalY;
rotatingZ = globalZ;
+
// rotate
if (angleX != 0) {
// rotate in x
--
GitLab