From 1a60ef65757500f2f0d081121ef6da42d656ccbf Mon Sep 17 00:00:00 2001
From: Martin Schoenherr <m.schoenherr@tu-braunschweig.de>
Date: Tue, 14 Nov 2023 17:13:40 +0100
Subject: [PATCH] fix init Advection Diffusion

---
 src/gpu/core/GPU/GPU_Interface.h              |  47 --
 src/gpu/core/GPU/GPU_Kernels.cuh              |  42 --
 src/gpu/core/GPU/Init27.cu                    | 512 -----------------
 src/gpu/core/GPU/InitAdvectionDiffusion27.cu  | 529 ------------------
 src/gpu/core/GPU/LBMKernel.cu                 | 105 ----
 src/gpu/core/Init/InitLattice.cpp             |   5 +-
 src/gpu/core/Init/InitLattice.h               |   2 +-
 .../core/KernelManager/ADKernelManager.cpp    |  50 --
 src/gpu/core/KernelManager/ADKernelManager.h  |   3 -
 src/gpu/core/LBM/Simulation.cpp               |   4 +-
 src/gpu/core/LBM/Simulation.h                 |   1 +
 11 files changed, 8 insertions(+), 1292 deletions(-)
 delete mode 100644 src/gpu/core/GPU/Init27.cu
 delete mode 100644 src/gpu/core/GPU/InitAdvectionDiffusion27.cu

diff --git a/src/gpu/core/GPU/GPU_Interface.h b/src/gpu/core/GPU/GPU_Interface.h
index 00894db28..72e839bcf 100644
--- a/src/gpu/core/GPU/GPU_Interface.h
+++ b/src/gpu/core/GPU/GPU_Interface.h
@@ -27,53 +27,6 @@
 struct LBMSimulationParameter;
 class Parameter;
 
-//////////////////////////////////////////////////////////////////////////
-//Kernel
-//////////////////////////////////////////////////////////////////////////
-void Init27(int myid,
-                       int numprocs,
-                       real u0,
-                       unsigned int* geoD,
-                       unsigned int* neighborX,
-                       unsigned int* neighborY,
-                       unsigned int* neighborZ,
-                       real* vParab,
-                       unsigned long long numberOfLBnodes,
-                       unsigned int grid_nx, 
-                       unsigned int grid_ny, 
-                       unsigned int grid_nz, 
-                       real* DD,
-                       int level,
-                       int maxlevel);
-
-void InitNonEqPartSP27(unsigned int numberOfThreads,
-                                  unsigned int* neighborX,
-                                  unsigned int* neighborY,
-                                  unsigned int* neighborZ,
-                                  unsigned int* neighborWSB,
-                                  unsigned int* geoD,
-                                  real* rho,
-                                  real* ux,
-                                  real* uy,
-                                  real* uz,
-                                  unsigned long long numberOfLBnodes,
-                                  real* DD,
-                                  real omega,
-                                  bool EvenOrOdd);
-
-
-void InitADDev27( unsigned int numberOfThreads,
-                           unsigned int* neighborX,
-                           unsigned int* neighborY,
-                           unsigned int* neighborZ,
-                           unsigned int* geoD,
-                           real* Conc,
-                           real* ux,
-                           real* uy,
-                           real* uz,
-                           unsigned long long numberOfLBnodes,
-                           real* DD27,
-                           bool EvenOrOdd);
 
 void CalcMac27( real* vxD,
                           real* vyD,
diff --git a/src/gpu/core/GPU/GPU_Kernels.cuh b/src/gpu/core/GPU/GPU_Kernels.cuh
index 2fd05b9e7..ab8ab0fa3 100644
--- a/src/gpu/core/GPU/GPU_Kernels.cuh
+++ b/src/gpu/core/GPU/GPU_Kernels.cuh
@@ -17,48 +17,6 @@
 #include "LBM/LB.h"
 
 
-__global__ void LBInit27( int myid,
-                                     int numprocs,
-                                     real u0,
-                                     unsigned int* geoD,
-                                     unsigned int* neighborX,
-                                     unsigned int* neighborY,
-                                     unsigned int* neighborZ,
-                                     real* vParabel,
-                                     unsigned long long numberOfLBnodes,
-                                     unsigned int grid_nx,
-                                     unsigned int grid_ny,
-                                     unsigned int grid_nz,
-                                     real* DD,
-                                     int lev,
-                                     int maxlev);
-
-__global__ void LBInitNonEqPartSP27(unsigned int* neighborX,
-                                               unsigned int* neighborY,
-                                               unsigned int* neighborZ,
-                                               unsigned int* neighborWSB,
-                                               unsigned int* geoD,
-                                               real* rho,
-                                               real* ux,
-                                               real* uy,
-                                               real* uz,
-                                               unsigned long long numberOfLBnodes,
-                                               real* DD,
-                                               real omega,
-                                               bool EvenOrOdd);
-
-__global__ void InitAD27(unsigned int* neighborX,
-                                       unsigned int* neighborY,
-                                       unsigned int* neighborZ,
-                                       unsigned int* geoD,
-                                       real* Conc,
-                                       real* ux,
-                                       real* uy,
-                                       real* uz,
-                                       unsigned long long numberOfLBnodes,
-                                       real* DD27,
-                                       bool EvenOrOdd);
-
 __global__ void LBCalcMac27( real* vxD,
                                         real* vyD,
                                         real* vzD,
diff --git a/src/gpu/core/GPU/Init27.cu b/src/gpu/core/GPU/Init27.cu
deleted file mode 100644
index 961045d3c..000000000
--- a/src/gpu/core/GPU/Init27.cu
+++ /dev/null
@@ -1,512 +0,0 @@
-/* Device code */
-#include "LBM/LB.h" 
-#include "lbm/constants/D3Q27.h"
-#include <basics/constants/NumericConstants.h>
-
-using namespace vf::basics::constant;
-using namespace vf::lbm::dir;
-
-////////////////////////////////////////////////////////////////////////////////
-__global__ void LBInit27( int myid,
-                                     int numprocs,
-                                     real u0,
-                                     unsigned int* geoD,
-                                     unsigned int* neighborX,
-                                     unsigned int* neighborY,
-                                     unsigned int* neighborZ,
-                                     real* vParabel,
-                                     unsigned long long numberOfLBnodes,
-                                     unsigned int grid_nx, 
-                                     unsigned int grid_ny, 
-                                     unsigned int grid_nz, 
-                                     real* DD,
-                                     int lev,
-                                     int maxlev)
-{
-   Distributions27 D;
-   D.f[dP00] = &DD[dP00 * numberOfLBnodes];
-   D.f[dM00] = &DD[dM00 * numberOfLBnodes];
-   D.f[d0P0] = &DD[d0P0 * numberOfLBnodes];
-   D.f[d0M0] = &DD[d0M0 * numberOfLBnodes];
-   D.f[d00P] = &DD[d00P * numberOfLBnodes];
-   D.f[d00M] = &DD[d00M * numberOfLBnodes];
-   D.f[dPP0] = &DD[dPP0 * numberOfLBnodes];
-   D.f[dMM0] = &DD[dMM0 * numberOfLBnodes];
-   D.f[dPM0] = &DD[dPM0 * numberOfLBnodes];
-   D.f[dMP0] = &DD[dMP0 * numberOfLBnodes];
-   D.f[dP0P] = &DD[dP0P * numberOfLBnodes];
-   D.f[dM0M] = &DD[dM0M * numberOfLBnodes];
-   D.f[dP0M] = &DD[dP0M * numberOfLBnodes];
-   D.f[dM0P] = &DD[dM0P * numberOfLBnodes];
-   D.f[d0PP] = &DD[d0PP * numberOfLBnodes];
-   D.f[d0MM] = &DD[d0MM * numberOfLBnodes];
-   D.f[d0PM] = &DD[d0PM * numberOfLBnodes];
-   D.f[d0MP] = &DD[d0MP * numberOfLBnodes];
-   D.f[d000] = &DD[d000 * numberOfLBnodes];
-   D.f[dPPP] = &DD[dPPP * numberOfLBnodes];
-   D.f[dMMP] = &DD[dMMP * numberOfLBnodes];
-   D.f[dPMP] = &DD[dPMP * numberOfLBnodes];
-   D.f[dMPP] = &DD[dMPP * numberOfLBnodes];
-   D.f[dPPM] = &DD[dPPM * numberOfLBnodes];
-   D.f[dMMM] = &DD[dMMM * numberOfLBnodes];
-   D.f[dPMM] = &DD[dPMM * numberOfLBnodes];
-   D.f[dMPM] = &DD[dMPM * numberOfLBnodes];
-   ////////////////////////////////////////////////////////////////////////////////
-   unsigned int  k;                   // Zugriff auf arrays im device
-   //
-   unsigned int tx = threadIdx.x;     // Thread index = lokaler i index
-   unsigned int by = blockIdx.x;      // Block index x
-   unsigned int bz = blockIdx.y;      // Block index y
-   unsigned int  x = tx + STARTOFFX;  // Globaler x-Index 
-   unsigned int  y = by + STARTOFFY;  // Globaler y-Index 
-   unsigned int  z = bz + STARTOFFZ;  // Globaler z-Index 
-
-   const unsigned sizeX = blockDim.x;
-   const unsigned sizeY = gridDim.x;
-   const unsigned nx = sizeX + 2 * STARTOFFX;
-   const unsigned ny = sizeY + 2 * STARTOFFY;
-
-   k = nx*(ny*z + y) + x;
-   //////////////////////////////////////////////////////////////////////////
-   geoD[k] = GEO_FLUID;
-   if (lev==0)
-   {
-      if( by == 0 || by == grid_ny-1 || tx == 0 || tx == grid_nx-1 )             
-         geoD[k] = GEO_SOLID;
-      else if( bz == grid_nz-1 && myid == numprocs - 1 && geoD[k] != GEO_SOLID )
-         geoD[k] = GEO_PRESS;                 
-      else if( bz == 0 && myid == 0 && geoD[k] != GEO_SOLID)
-         geoD[k] = GEO_SOLID;//GEO_VELO;
-   }
-   else if (lev==maxlev-1)
-   {
-      unsigned int centerX = grid_nx / 2;
-      unsigned int centerY = grid_ny / 2;
-      unsigned int centerZ = grid_nz / 2;
-      real        radius  = grid_ny / 2.56;
-
-      unsigned int distSq = (centerX-tx)*(centerX-tx)+(centerY-by)*(centerY-by)+(centerZ-bz)*(centerZ-bz);
-      real radiSq = radius*radius;
-
-      if( distSq < radiSq)        geoD[k] = GEO_SOLID;
-   }
-   //////////////////////////////////////////////////////////////////////////
-   real drho = c0o1;
-   real  vx1 = c0o1;
-   real  vx2 = c0o1;
-   real  vx3 = u0;
-   vParabel[k] = vx3;
-   ////////////////////////////////////////////////////////////////////////////////
-   //index
-   unsigned int nxny = nx*ny;
-   ////////////////////////////////////////////////////////////////////////////////
-   //neighborX[k]      = k+1;
-   //neighborY[k+1]    = k+nx+1;
-   //neighborZ[k+1]    = k+nxny+1;
-   //neighborY[k]      = k+nx;
-   //neighborX[k+nx]   = k+nx+1;
-   //neighborZ[k+nx]   = k+nx+nxny;
-   //neighborZ[k]      = k+nxny;
-   //neighborX[k+nxny] = k+nxny+1;
-   //neighborY[k+nxny] = k+nxny+nx;
-   ////////////////////////////////////////////////////////////////////////////////
-   unsigned int kzero= k;
-   unsigned int ke   = k;
-   unsigned int kw   = k + 1;
-   unsigned int kn   = k;
-   unsigned int ks   = k + nx;
-   unsigned int kt   = k;
-   unsigned int kb   = k + nxny;
-   unsigned int ksw  = k + nx + 1;
-   unsigned int kne  = k;
-   unsigned int kse  = k + nx;
-   unsigned int knw  = k + 1;
-   unsigned int kbw  = k + nxny + 1;
-   unsigned int kte  = k;
-   unsigned int kbe  = k + nxny;
-   unsigned int ktw  = k + 1;
-   unsigned int kbs  = k + nxny + nx;
-   unsigned int ktn  = k;
-   unsigned int kbn  = k + nxny;
-   unsigned int kts  = k + nx;
-   unsigned int ktse = k + nx;
-   unsigned int kbnw = k + nxny + 1;
-   unsigned int ktnw = k + 1;
-   unsigned int kbse = k + nxny + nx;
-   unsigned int ktsw = k + nx + 1;
-   unsigned int kbne = k + nxny;
-   unsigned int ktne = k;
-   unsigned int kbsw = k + nxny + nx + 1;
-   //////////////////////////////////////////////////////////////////////////
-
-   real cu_sq=c3o2*(vx1*vx1+vx2*vx2+vx3*vx3);
-
-   (D.f[d000])[kzero] =   c8o27* (drho-cu_sq);
-   (D.f[dP00])[ke   ] =   c2o27* (drho+c3o1*( vx1        )+c9o2*( vx1        )*( vx1        )-cu_sq);
-   (D.f[dM00])[kw   ] =   c2o27* (drho+c3o1*(-vx1        )+c9o2*(-vx1        )*(-vx1        )-cu_sq);
-   (D.f[d0P0])[kn   ] =   c2o27* (drho+c3o1*(    vx2     )+c9o2*(     vx2    )*(     vx2    )-cu_sq);
-   (D.f[d0M0])[ks   ] =   c2o27* (drho+c3o1*(   -vx2     )+c9o2*(    -vx2    )*(    -vx2    )-cu_sq);
-   (D.f[d00P])[kt   ] =   c2o27* (drho+c3o1*(         vx3)+c9o2*(         vx3)*(         vx3)-cu_sq);
-   (D.f[d00M])[kb   ] =   c2o27* (drho+c3o1*(        -vx3)+c9o2*(        -vx3)*(        -vx3)-cu_sq);
-   (D.f[dPP0])[kne  ] =   c1o54* (drho+c3o1*( vx1+vx2    )+c9o2*( vx1+vx2    )*( vx1+vx2    )-cu_sq);
-   (D.f[dMM0])[ksw  ] =   c1o54* (drho+c3o1*(-vx1-vx2    )+c9o2*(-vx1-vx2    )*(-vx1-vx2    )-cu_sq);
-   (D.f[dPM0])[kse  ] =   c1o54* (drho+c3o1*( vx1-vx2    )+c9o2*( vx1-vx2    )*( vx1-vx2    )-cu_sq);
-   (D.f[dMP0])[knw  ] =   c1o54* (drho+c3o1*(-vx1+vx2    )+c9o2*(-vx1+vx2    )*(-vx1+vx2    )-cu_sq);
-   (D.f[dP0P])[kte  ] =   c1o54* (drho+c3o1*( vx1    +vx3)+c9o2*( vx1    +vx3)*( vx1    +vx3)-cu_sq);
-   (D.f[dM0M])[kbw  ] =   c1o54* (drho+c3o1*(-vx1    -vx3)+c9o2*(-vx1    -vx3)*(-vx1    -vx3)-cu_sq);
-   (D.f[dP0M])[kbe  ] =   c1o54* (drho+c3o1*( vx1    -vx3)+c9o2*( vx1    -vx3)*( vx1    -vx3)-cu_sq);
-   (D.f[dM0P])[ktw  ] =   c1o54* (drho+c3o1*(-vx1    +vx3)+c9o2*(-vx1    +vx3)*(-vx1    +vx3)-cu_sq);
-   (D.f[d0PP])[ktn  ] =   c1o54* (drho+c3o1*(     vx2+vx3)+c9o2*(     vx2+vx3)*(     vx2+vx3)-cu_sq);
-   (D.f[d0MM])[kbs  ] =   c1o54* (drho+c3o1*(    -vx2-vx3)+c9o2*(    -vx2-vx3)*(    -vx2-vx3)-cu_sq);
-   (D.f[d0PM])[kbn  ] =   c1o54* (drho+c3o1*(     vx2-vx3)+c9o2*(     vx2-vx3)*(     vx2-vx3)-cu_sq);
-   (D.f[d0MP])[kts  ] =   c1o54* (drho+c3o1*(    -vx2+vx3)+c9o2*(    -vx2+vx3)*(    -vx2+vx3)-cu_sq);
-   (D.f[dPPP])[ktne ] =   c1o216*(drho+c3o1*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
-   (D.f[dMMM])[kbsw ] =   c1o216*(drho+c3o1*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
-   (D.f[dPPM])[kbne ] =   c1o216*(drho+c3o1*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
-   (D.f[dMMP])[ktsw ] =   c1o216*(drho+c3o1*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
-   (D.f[dPMP])[ktse ] =   c1o216*(drho+c3o1*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
-   (D.f[dMPM])[kbnw ] =   c1o216*(drho+c3o1*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
-   (D.f[dPMM])[kbse ] =   c1o216*(drho+c3o1*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
-   (D.f[dMPP])[ktnw ] =   c1o216*(drho+c3o1*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
-
-}
-////////////////////////////////////////////////////////////////////////////////
-
-
-
-
-
-
-
-
-
-
-////////////////////////////////////////////////////////////////////////////////
-__global__ void LBInitNonEqPartSP27( unsigned int* neighborX,
-                                                unsigned int* neighborY,
-                                                unsigned int* neighborZ,
-                                                unsigned int* neighborWSB,
-                                                unsigned int* geoD,
-                                                real* rho,
-                                                real* ux,
-                                                real* uy,
-                                                real* uz,
-                                                unsigned long long numberOfLBnodes,
-                                                real* DD,
-                                                real omega,
-                                                bool EvenOrOdd)
-{
-    ////////////////////////////////////////////////////////////////////////////////
-    const unsigned  x = threadIdx.x;  // Globaler x-Index 
-    const unsigned  y = blockIdx.x;   // Globaler y-Index 
-    const unsigned  z = blockIdx.y;   // Globaler z-Index 
-    
-    const unsigned nx = blockDim.x;
-    const unsigned ny = gridDim.x;
-    
-    const unsigned k = nx*(ny*z + y) + x;
-    //////////////////////////////////////////////////////////////////////////
-    
-    if(k<numberOfLBnodes)
-    {
-        ////////////////////////////////////////////////////////////////////////////////
-        unsigned int BC;
-        BC = geoD[k];
-
-        if( BC != GEO_SOLID &&  BC != GEO_VOID)
-        {
-            Distributions27 D;
-            if (EvenOrOdd==true)
-            {
-                D.f[dP00] = &DD[dP00 * numberOfLBnodes];
-                D.f[dM00] = &DD[dM00 * numberOfLBnodes];
-                D.f[d0P0] = &DD[d0P0 * numberOfLBnodes];
-                D.f[d0M0] = &DD[d0M0 * numberOfLBnodes];
-                D.f[d00P] = &DD[d00P * numberOfLBnodes];
-                D.f[d00M] = &DD[d00M * numberOfLBnodes];
-                D.f[dPP0] = &DD[dPP0 * numberOfLBnodes];
-                D.f[dMM0] = &DD[dMM0 * numberOfLBnodes];
-                D.f[dPM0] = &DD[dPM0 * numberOfLBnodes];
-                D.f[dMP0] = &DD[dMP0 * numberOfLBnodes];
-                D.f[dP0P] = &DD[dP0P * numberOfLBnodes];
-                D.f[dM0M] = &DD[dM0M * numberOfLBnodes];
-                D.f[dP0M] = &DD[dP0M * numberOfLBnodes];
-                D.f[dM0P] = &DD[dM0P * numberOfLBnodes];
-                D.f[d0PP] = &DD[d0PP * numberOfLBnodes];
-                D.f[d0MM] = &DD[d0MM * numberOfLBnodes];
-                D.f[d0PM] = &DD[d0PM * numberOfLBnodes];
-                D.f[d0MP] = &DD[d0MP * numberOfLBnodes];
-                D.f[d000] = &DD[d000 * numberOfLBnodes];
-                D.f[dPPP] = &DD[dPPP * numberOfLBnodes];
-                D.f[dMMP] = &DD[dMMP * numberOfLBnodes];
-                D.f[dPMP] = &DD[dPMP * numberOfLBnodes];
-                D.f[dMPP] = &DD[dMPP * numberOfLBnodes];
-                D.f[dPPM] = &DD[dPPM * numberOfLBnodes];
-                D.f[dMMM] = &DD[dMMM * numberOfLBnodes];
-                D.f[dPMM] = &DD[dPMM * numberOfLBnodes];
-                D.f[dMPM] = &DD[dMPM * numberOfLBnodes];
-            }
-            else
-            {
-                D.f[dM00] = &DD[dP00 * numberOfLBnodes];
-                D.f[dP00] = &DD[dM00 * numberOfLBnodes];
-                D.f[d0M0] = &DD[d0P0 * numberOfLBnodes];
-                D.f[d0P0] = &DD[d0M0 * numberOfLBnodes];
-                D.f[d00M] = &DD[d00P * numberOfLBnodes];
-                D.f[d00P] = &DD[d00M * numberOfLBnodes];
-                D.f[dMM0] = &DD[dPP0 * numberOfLBnodes];
-                D.f[dPP0] = &DD[dMM0 * numberOfLBnodes];
-                D.f[dMP0] = &DD[dPM0 * numberOfLBnodes];
-                D.f[dPM0] = &DD[dMP0 * numberOfLBnodes];
-                D.f[dM0M] = &DD[dP0P * numberOfLBnodes];
-                D.f[dP0P] = &DD[dM0M * numberOfLBnodes];
-                D.f[dM0P] = &DD[dP0M * numberOfLBnodes];
-                D.f[dP0M] = &DD[dM0P * numberOfLBnodes];
-                D.f[d0MM] = &DD[d0PP * numberOfLBnodes];
-                D.f[d0PP] = &DD[d0MM * numberOfLBnodes];
-                D.f[d0MP] = &DD[d0PM * numberOfLBnodes];
-                D.f[d0PM] = &DD[d0MP * numberOfLBnodes];
-                D.f[d000] = &DD[d000 * numberOfLBnodes];
-                D.f[dMMM] = &DD[dPPP * numberOfLBnodes];
-                D.f[dPPM] = &DD[dMMP * numberOfLBnodes];
-                D.f[dMPM] = &DD[dPMP * numberOfLBnodes];
-                D.f[dPMM] = &DD[dMPP * numberOfLBnodes];
-                D.f[dMMP] = &DD[dPPM * numberOfLBnodes];
-                D.f[dPPP] = &DD[dMMM * numberOfLBnodes];
-                D.f[dMPP] = &DD[dPMM * numberOfLBnodes];
-                D.f[dPMP] = &DD[dMPM * numberOfLBnodes];
-            }
-            //////////////////////////////////////////////////////////////////////////
-            real drho = rho[k];//0.0f;//
-            real  vx1 = ux[k]; //0.0f;//
-            real  vx2 = uy[k]; //0.0f;//
-            real  vx3 = uz[k]; //0.0f;//
-            //////////////////////////////////////////////////////////////////////////
-            //index
-            //////////////////////////////////////////////////////////////////////////
-            unsigned int kzero= k;
-            unsigned int ke   = k;
-            unsigned int kw   = neighborX[k];
-            unsigned int kn   = k;
-            unsigned int ks   = neighborY[k];
-            unsigned int kt   = k;
-            unsigned int kb   = neighborZ[k];
-            unsigned int ksw  = neighborY[kw];
-            unsigned int kne  = k;
-            unsigned int kse  = ks;
-            unsigned int knw  = kw;
-            unsigned int kbw  = neighborZ[kw];
-            unsigned int kte  = k;
-            unsigned int kbe  = kb;
-            unsigned int ktw  = kw;
-            unsigned int kbs  = neighborZ[ks];
-            unsigned int ktn  = k;
-            unsigned int kbn  = kb;
-            unsigned int kts  = ks;
-            unsigned int ktse = ks;
-            unsigned int kbnw = kbw;
-            unsigned int ktnw = kw;
-            unsigned int kbse = kbs;
-            unsigned int ktsw = ksw;
-            unsigned int kbne = kb;
-            unsigned int ktne = k;
-            unsigned int kbsw = neighborZ[ksw];
-            //////////////////////////////////////////////////////////////////////////////
-            //neighbor index
-            uint kPx   = neighborX[k];
-            uint kPy   = neighborY[k];
-            uint kPz   = neighborZ[k];
-            uint kMxyz = neighborWSB[k];
-            uint kMx   = neighborZ[neighborY[kMxyz]];
-            uint kMy   = neighborZ[neighborX[kMxyz]];
-            uint kMz   = neighborY[neighborX[kMxyz]];
-            //////////////////////////////////////////////////////////////////////////
-            //getVeloX//
-            real vx1NeighborPx = ux[kPx];
-            real vx1NeighborMx = ux[kMx];
-            real vx1NeighborPy = ux[kPy];
-            real vx1NeighborMy = ux[kMy];
-            real vx1NeighborPz = ux[kPz];
-            real vx1NeighborMz = ux[kMz];
-            //getVeloY//
-            real vx2NeighborPx = uy[kPx];
-            real vx2NeighborMx = uy[kMx];
-            real vx2NeighborPy = uy[kPy];
-            real vx2NeighborMy = uy[kMy];
-            real vx2NeighborPz = uy[kPz];
-            real vx2NeighborMz = uy[kMz];
-            //getVeloZ//
-            real vx3NeighborPx = uz[kPx];
-            real vx3NeighborMx = uz[kMx];
-            real vx3NeighborPy = uz[kPy];
-            real vx3NeighborMy = uz[kMy];
-            real vx3NeighborPz = uz[kPz];
-            real vx3NeighborMz = uz[kMz];
-            //////////////////////////////////////////////////////////////////////////
-
-            real dvx1dx = (vx1NeighborPx - vx1NeighborMx) / c2o1;
-            real dvx1dy = (vx1NeighborPy - vx1NeighborMy) / c2o1;
-            real dvx1dz = (vx1NeighborPz - vx1NeighborMz) / c2o1;
-
-            real dvx2dx = (vx2NeighborPx - vx2NeighborMx) / c2o1;
-            real dvx2dy = (vx2NeighborPy - vx2NeighborMy) / c2o1;
-            real dvx2dz = (vx2NeighborPz - vx2NeighborMz) / c2o1;
-
-            real dvx3dx = (vx3NeighborPx - vx3NeighborMx) / c2o1;
-            real dvx3dy = (vx3NeighborPy - vx3NeighborMy) / c2o1;
-            real dvx3dz = (vx3NeighborPz - vx3NeighborMz) / c2o1;
-
-            //////////////////////////////////////////////////////////////////////////
-
-            // the following code is copy and pasted from VirtualFluidsVisitors/InitDistributionsBlockVisitor.cpp
-            // i.e. Konstantins code
-
-            real ax = dvx1dx;
-            real ay = dvx1dy;
-            real az = dvx1dz;
-
-            real bx = dvx2dx;
-            real by = dvx2dy;
-            real bz = dvx2dz;
-
-            real cx = dvx3dx;
-            real cy = dvx3dy;
-            real cz = dvx3dz;
-
-            real eps_new = c1o1;
-            real op      = c1o1;
-            real o       = omega;
-
-            real f_E    = eps_new *((5.*ax*omega + 5.*by*o + 5.*cz*o - 8.*ax*op + 4.*by*op + 4.*cz*op)/(54.*o*op));
-
-            real f_N    =    f_E   + eps_new *((2.*(ax - by))/(9.*o));
-            real f_T    =    f_E   + eps_new *((2.*(ax - cz))/(9.*o));
-            real f_NE   =            eps_new *(-(5.*cz*o + 3.*(ay + bx)*op - 2.*cz*op + ax*(5.*o + op) + by*(5.*o + op))/(54.*o*op));
-            real f_SE   =    f_NE  + eps_new *((  ay + bx )/(9.*o));
-            real f_TE   =            eps_new *(-(5.*cz*o + by*(5.*o - 2.*op) + 3.*(az + cx)*op + cz*op + ax*(5.*o + op))/(54.*o*op));
-            real f_BE   =    f_TE  + eps_new *((  az + cx )/(9.*o));
-            real f_TN   =            eps_new *(-(5.*ax*o + 5.*by*o + 5.*cz*o - 2.*ax*op + by*op + 3.*bz*op + 3.*cy*op + cz*op)/(54.*o*op));
-            real f_BN   =    f_TN  + eps_new *((  bz + cy )/(9.*o));
-            real f_ZERO =            eps_new *((5.*(ax + by + cz))/(9.*op));
-            real f_TNE  =            eps_new *(-(ay + az + bx + bz + cx + cy)/(72.*o));
-            real f_TSW  =  - f_TNE - eps_new *((ay + bx)/(36.*o));
-            real f_TSE  =  - f_TNE - eps_new *((az + cx)/(36.*o));
-            real f_TNW  =  - f_TNE - eps_new *((bz + cy)/(36.*o));
-
-            //////////////////////////////////////////////////////////////////////////
-            real cu_sq=c3o2*(vx1*vx1+vx2*vx2+vx3*vx3);
-            
-            (D.f[d000])[kzero] =   c8o27* (drho-cu_sq);
-            (D.f[dP00])[ke   ] =   c2o27* (drho+c3o1*( vx1        )+c9o2*( vx1        )*( vx1        )-cu_sq);
-            (D.f[dM00])[kw   ] =   c2o27* (drho+c3o1*(-vx1        )+c9o2*(-vx1        )*(-vx1        )-cu_sq);
-            (D.f[d0P0])[kn   ] =   c2o27* (drho+c3o1*(    vx2     )+c9o2*(     vx2    )*(     vx2    )-cu_sq);
-            (D.f[d0M0])[ks   ] =   c2o27* (drho+c3o1*(   -vx2     )+c9o2*(    -vx2    )*(    -vx2    )-cu_sq);
-            (D.f[d00P])[kt   ] =   c2o27* (drho+c3o1*(         vx3)+c9o2*(         vx3)*(         vx3)-cu_sq);
-            (D.f[d00M])[kb   ] =   c2o27* (drho+c3o1*(        -vx3)+c9o2*(        -vx3)*(        -vx3)-cu_sq);
-            (D.f[dPP0])[kne  ] =   c1o54* (drho+c3o1*( vx1+vx2    )+c9o2*( vx1+vx2    )*( vx1+vx2    )-cu_sq);
-            (D.f[dMM0])[ksw  ] =   c1o54* (drho+c3o1*(-vx1-vx2    )+c9o2*(-vx1-vx2    )*(-vx1-vx2    )-cu_sq);
-            (D.f[dPM0])[kse  ] =   c1o54* (drho+c3o1*( vx1-vx2    )+c9o2*( vx1-vx2    )*( vx1-vx2    )-cu_sq);
-            (D.f[dMP0])[knw  ] =   c1o54* (drho+c3o1*(-vx1+vx2    )+c9o2*(-vx1+vx2    )*(-vx1+vx2    )-cu_sq);
-            (D.f[dP0P])[kte  ] =   c1o54* (drho+c3o1*( vx1    +vx3)+c9o2*( vx1    +vx3)*( vx1    +vx3)-cu_sq);
-            (D.f[dM0M])[kbw  ] =   c1o54* (drho+c3o1*(-vx1    -vx3)+c9o2*(-vx1    -vx3)*(-vx1    -vx3)-cu_sq);
-            (D.f[dP0M])[kbe  ] =   c1o54* (drho+c3o1*( vx1    -vx3)+c9o2*( vx1    -vx3)*( vx1    -vx3)-cu_sq);
-            (D.f[dM0P])[ktw  ] =   c1o54* (drho+c3o1*(-vx1    +vx3)+c9o2*(-vx1    +vx3)*(-vx1    +vx3)-cu_sq);
-            (D.f[d0PP])[ktn  ] =   c1o54* (drho+c3o1*(     vx2+vx3)+c9o2*(     vx2+vx3)*(     vx2+vx3)-cu_sq);
-            (D.f[d0MM])[kbs  ] =   c1o54* (drho+c3o1*(    -vx2-vx3)+c9o2*(    -vx2-vx3)*(    -vx2-vx3)-cu_sq);
-            (D.f[d0PM])[kbn  ] =   c1o54* (drho+c3o1*(     vx2-vx3)+c9o2*(     vx2-vx3)*(     vx2-vx3)-cu_sq);
-            (D.f[d0MP])[kts  ] =   c1o54* (drho+c3o1*(    -vx2+vx3)+c9o2*(    -vx2+vx3)*(    -vx2+vx3)-cu_sq);
-            (D.f[dPPP])[ktne ] =   c1o216*(drho+c3o1*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
-            (D.f[dMMM])[kbsw ] =   c1o216*(drho+c3o1*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
-            (D.f[dPPM])[kbne ] =   c1o216*(drho+c3o1*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
-            (D.f[dMMP])[ktsw ] =   c1o216*(drho+c3o1*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
-            (D.f[dPMP])[ktse ] =   c1o216*(drho+c3o1*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
-            (D.f[dMPM])[kbnw ] =   c1o216*(drho+c3o1*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
-            (D.f[dPMM])[kbse ] =   c1o216*(drho+c3o1*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
-            (D.f[dMPP])[ktnw ] =   c1o216*(drho+c3o1*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
-
-            //////////////////////////////////////////////////////////////////////////
-
-            (D.f[d000])[kzero] += f_ZERO;
-            (D.f[dP00])[ke   ] += f_E   ;
-            (D.f[dM00])[kw   ] += f_E   ;
-            (D.f[d0P0])[kn   ] += f_N   ;
-            (D.f[d0M0])[ks   ] += f_N   ;
-            (D.f[d00P])[kt   ] += f_T   ;
-            (D.f[d00M])[kb   ] += f_T   ;
-            (D.f[dPP0])[kne  ] += f_NE  ;
-            (D.f[dMM0])[ksw  ] += f_NE  ;
-            (D.f[dPM0])[kse  ] += f_SE  ;
-            (D.f[dMP0])[knw  ] += f_SE  ;
-            (D.f[dP0P])[kte  ] += f_TE  ;
-            (D.f[dM0M])[kbw  ] += f_TE  ;
-            (D.f[dP0M])[kbe  ] += f_BE  ;
-            (D.f[dM0P])[ktw  ] += f_BE  ;
-            (D.f[d0PP])[ktn  ] += f_TN  ;
-            (D.f[d0MM])[kbs  ] += f_TN  ;
-            (D.f[d0PM])[kbn  ] += f_BN  ;
-            (D.f[d0MP])[kts  ] += f_BN  ;
-            (D.f[dPPP])[ktne ] += f_TNE ;
-            (D.f[dMMM])[kbsw ] += f_TNE ;
-            (D.f[dPPM])[kbne ] += f_TSW ;
-            (D.f[dMMP])[ktsw ] += f_TSW ;
-            (D.f[dPMP])[ktse ] += f_TSE ;
-            (D.f[dMPM])[kbnw ] += f_TSE ;
-            (D.f[dPMM])[kbse ] += f_TNW ;
-            (D.f[dMPP])[ktnw ] += f_TNW ;
-
-            //////////////////////////////////////////////////////////////////////////
-        }
-        else
-        {
-            //////////////////////////////////////////////////////////////////////////
-            Distributions27 D;
-            D.f[d000] = &DD[d000 * numberOfLBnodes];
-            //////////////////////////////////////////////////////////////////////////
-            (D.f[d000])[k] = c96o1;
-            //////////////////////////////////////////////////////////////////////////
-        }
-   }
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/src/gpu/core/GPU/InitAdvectionDiffusion27.cu b/src/gpu/core/GPU/InitAdvectionDiffusion27.cu
deleted file mode 100644
index 8e2118e25..000000000
--- a/src/gpu/core/GPU/InitAdvectionDiffusion27.cu
+++ /dev/null
@@ -1,529 +0,0 @@
-//=======================================================================================
-// ____          ____    __    ______     __________   __      __       __        __
-// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
-//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
-//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
-//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
-//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
-//      \    \  |    |   ________________________________________________________________
-//       \    \ |    |  |  ______________________________________________________________|
-//        \    \|    |  |  |         __          __     __     __     ______      _______
-//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
-//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
-//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
-//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
-//
-//  This file is part of VirtualFluids. VirtualFluids is free software: you can
-//  redistribute it and/or modify it under the terms of the GNU General Public
-//  License as published by the Free Software Foundation, either version 3 of
-//  the License, or (at your option) any later version.
-//
-//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
-//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-//  for more details.
-//
-//  You should have received a copy of the GNU General Public License along
-//  with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
-//
-//! \file InitAdvectionDiffusion.cu
-//! \ingroup GPU
-//! \author Martin Schoenherr
-//=======================================================================================
-/* Device code */
-#include "LBM/LB.h"
-#include "lbm/constants/D3Q27.h"
-#include <basics/constants/NumericConstants.h>
-
-using namespace vf::basics::constant;
-using namespace vf::lbm::dir;
-
-__global__ void InitAD27(
-    uint* neighborX,
-    uint* neighborY,
-    uint* neighborZ,
-    uint* typeOfGridNode,
-    real* concentration,
-    real* velocityX,
-    real* velocityY,
-    real* velocityZ,
-    unsigned long long numberOfLBnodes,
-    real* distributionsAD,
-    bool isEvenTimestep)
-{
-    //////////////////////////////////////////////////////////////////////////
-    //! The initialization is executed in the following steps
-    //!
-    ////////////////////////////////////////////////////////////////////////////////
-    //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
-    //!
-    const unsigned  x = threadIdx.x;  // Globaler x-Index
-    const unsigned  y = blockIdx.x;   // Globaler y-Index
-    const unsigned  z = blockIdx.y;   // Globaler z-Index
-
-    const unsigned nx = blockDim.x;
-    const unsigned ny = gridDim.x;
-
-    const unsigned k = nx*(ny*z + y) + x;
-
-    //////////////////////////////////////////////////////////////////////////
-    // run for all indices in size_Mat and fluid nodes
-    if ((k < numberOfLBnodes) && (typeOfGridNode[k] == GEO_FLUID))
-    {
-        //////////////////////////////////////////////////////////////////////////
-        //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on timestep is based on the esoteric twist algorithm \ref
-        //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
-        //!
-        Distributions27 distAD;
-        if (isEvenTimestep)
-        {
-            distAD.f[dP00] = &distributionsAD[dP00 * numberOfLBnodes];
-            distAD.f[dM00] = &distributionsAD[dM00 * numberOfLBnodes];
-            distAD.f[d0P0] = &distributionsAD[d0P0 * numberOfLBnodes];
-            distAD.f[d0M0] = &distributionsAD[d0M0 * numberOfLBnodes];
-            distAD.f[d00P] = &distributionsAD[d00P * numberOfLBnodes];
-            distAD.f[d00M] = &distributionsAD[d00M * numberOfLBnodes];
-            distAD.f[dPP0] = &distributionsAD[dPP0 * numberOfLBnodes];
-            distAD.f[dMM0] = &distributionsAD[dMM0 * numberOfLBnodes];
-            distAD.f[dPM0] = &distributionsAD[dPM0 * numberOfLBnodes];
-            distAD.f[dMP0] = &distributionsAD[dMP0 * numberOfLBnodes];
-            distAD.f[dP0P] = &distributionsAD[dP0P * numberOfLBnodes];
-            distAD.f[dM0M] = &distributionsAD[dM0M * numberOfLBnodes];
-            distAD.f[dP0M] = &distributionsAD[dP0M * numberOfLBnodes];
-            distAD.f[dM0P] = &distributionsAD[dM0P * numberOfLBnodes];
-            distAD.f[d0PP] = &distributionsAD[d0PP * numberOfLBnodes];
-            distAD.f[d0MM] = &distributionsAD[d0MM * numberOfLBnodes];
-            distAD.f[d0PM] = &distributionsAD[d0PM * numberOfLBnodes];
-            distAD.f[d0MP] = &distributionsAD[d0MP * numberOfLBnodes];
-            distAD.f[d000] = &distributionsAD[d000 * numberOfLBnodes];
-            distAD.f[dPPP] = &distributionsAD[dPPP * numberOfLBnodes];
-            distAD.f[dMMP] = &distributionsAD[dMMP * numberOfLBnodes];
-            distAD.f[dPMP] = &distributionsAD[dPMP * numberOfLBnodes];
-            distAD.f[dMPP] = &distributionsAD[dMPP * numberOfLBnodes];
-            distAD.f[dPPM] = &distributionsAD[dPPM * numberOfLBnodes];
-            distAD.f[dMMM] = &distributionsAD[dMMM * numberOfLBnodes];
-            distAD.f[dPMM] = &distributionsAD[dPMM * numberOfLBnodes];
-            distAD.f[dMPM] = &distributionsAD[dMPM * numberOfLBnodes];
-        }
-        else
-        {
-            distAD.f[dM00] = &distributionsAD[dP00 * numberOfLBnodes];
-            distAD.f[dP00] = &distributionsAD[dM00 * numberOfLBnodes];
-            distAD.f[d0M0] = &distributionsAD[d0P0 * numberOfLBnodes];
-            distAD.f[d0P0] = &distributionsAD[d0M0 * numberOfLBnodes];
-            distAD.f[d00M] = &distributionsAD[d00P * numberOfLBnodes];
-            distAD.f[d00P] = &distributionsAD[d00M * numberOfLBnodes];
-            distAD.f[dMM0] = &distributionsAD[dPP0 * numberOfLBnodes];
-            distAD.f[dPP0] = &distributionsAD[dMM0 * numberOfLBnodes];
-            distAD.f[dMP0] = &distributionsAD[dPM0 * numberOfLBnodes];
-            distAD.f[dPM0] = &distributionsAD[dMP0 * numberOfLBnodes];
-            distAD.f[dM0M] = &distributionsAD[dP0P * numberOfLBnodes];
-            distAD.f[dP0P] = &distributionsAD[dM0M * numberOfLBnodes];
-            distAD.f[dM0P] = &distributionsAD[dP0M * numberOfLBnodes];
-            distAD.f[dP0M] = &distributionsAD[dM0P * numberOfLBnodes];
-            distAD.f[d0MM] = &distributionsAD[d0PP * numberOfLBnodes];
-            distAD.f[d0PP] = &distributionsAD[d0MM * numberOfLBnodes];
-            distAD.f[d0MP] = &distributionsAD[d0PM * numberOfLBnodes];
-            distAD.f[d0PM] = &distributionsAD[d0MP * numberOfLBnodes];
-            distAD.f[d000] = &distributionsAD[d000 * numberOfLBnodes];
-            distAD.f[dMMM] = &distributionsAD[dPPP * numberOfLBnodes];
-            distAD.f[dPPM] = &distributionsAD[dMMP * numberOfLBnodes];
-            distAD.f[dMPM] = &distributionsAD[dPMP * numberOfLBnodes];
-            distAD.f[dPMM] = &distributionsAD[dMPP * numberOfLBnodes];
-            distAD.f[dMMP] = &distributionsAD[dPPM * numberOfLBnodes];
-            distAD.f[dPPP] = &distributionsAD[dMMM * numberOfLBnodes];
-            distAD.f[dMPP] = &distributionsAD[dPMM * numberOfLBnodes];
-            distAD.f[dPMP] = &distributionsAD[dMPM * numberOfLBnodes];
-        }
-        //////////////////////////////////////////////////////////////////////////
-        //! - Set local velocities and concetration
-        //!
-        real conc = concentration[k];
-        real  vx1 = velocityX[k];
-        real  vx2 = velocityY[k];
-        real  vx3 = velocityZ[k];
-        //////////////////////////////////////////////////////////////////////////
-        //! - Set neighbor indices (necessary for indirect addressing)
-        //!
-        uint kzero = k;
-        uint ke    = k;
-        uint kw    = neighborX[k];
-        uint kn    = k;
-        uint ks    = neighborY[k];
-        uint kt    = k;
-        uint kb    = neighborZ[k];
-        uint ksw   = neighborY[kw];
-        uint kne   = k;
-        uint kse   = ks;
-        uint knw   = kw;
-        uint kbw   = neighborZ[kw];
-        uint kte   = k;
-        uint kbe   = kb;
-        uint ktw   = kw;
-        uint kbs   = neighborZ[ks];
-        uint ktn   = k;
-        uint kbn   = kb;
-        uint kts   = ks;
-        uint ktse  = ks;
-        uint kbnw  = kbw;
-        uint ktnw  = kw;
-        uint kbse  = kbs;
-        uint ktsw  = ksw;
-        uint kbne  = kb;
-        uint ktne  = k;
-        uint kbsw  = neighborZ[ksw];
-        //////////////////////////////////////////////////////////////////////////
-        //! - Calculate the equilibrium and set the distributions
-        //!
-        real cu_sq = c3o2*(vx1*vx1 + vx2*vx2 + vx3*vx3);
-
-        (distAD.f[d000])[kzero] = c8o27  * conc * (c1o1 - cu_sq);
-        (distAD.f[dP00])[ke   ] = c2o27  * conc * (c1o1 + c3o1 * ( vx1            ) + c9o2 * ( vx1            ) * ( vx1            ) - cu_sq);
-        (distAD.f[dM00])[kw   ] = c2o27  * conc * (c1o1 + c3o1 * (-vx1            ) + c9o2 * (-vx1            ) * (-vx1            ) - cu_sq);
-        (distAD.f[d0P0])[kn   ] = c2o27  * conc * (c1o1 + c3o1 * (       vx2      ) + c9o2 * (       vx2      ) * (       vx2      ) - cu_sq);
-        (distAD.f[d0M0])[ks   ] = c2o27  * conc * (c1o1 + c3o1 * (     - vx2      ) + c9o2 * (     - vx2      ) * (     - vx2      ) - cu_sq);
-        (distAD.f[d00P])[kt   ] = c2o27  * conc * (c1o1 + c3o1 * (             vx3) + c9o2 * (             vx3) * (             vx3) - cu_sq);
-        (distAD.f[d00M])[kb   ] = c2o27  * conc * (c1o1 + c3o1 * (           - vx3) + c9o2 * (           - vx3) * (           - vx3) - cu_sq);
-        (distAD.f[dPP0])[kne  ] = c1o54  * conc * (c1o1 + c3o1 * ( vx1 + vx2      ) + c9o2 * ( vx1 + vx2      ) * ( vx1 + vx2      ) - cu_sq);
-        (distAD.f[dMM0])[ksw  ] = c1o54  * conc * (c1o1 + c3o1 * (-vx1 - vx2      ) + c9o2 * (-vx1 - vx2      ) * (-vx1 - vx2      ) - cu_sq);
-        (distAD.f[dPM0])[kse  ] = c1o54  * conc * (c1o1 + c3o1 * ( vx1 - vx2      ) + c9o2 * ( vx1 - vx2      ) * ( vx1 - vx2      ) - cu_sq);
-        (distAD.f[dMP0])[knw  ] = c1o54  * conc * (c1o1 + c3o1 * (-vx1 + vx2      ) + c9o2 * (-vx1 + vx2      ) * (-vx1 + vx2      ) - cu_sq);
-        (distAD.f[dP0P])[kte  ] = c1o54  * conc * (c1o1 + c3o1 * ( vx1       + vx3) + c9o2 * ( vx1       + vx3) * ( vx1       + vx3) - cu_sq);
-        (distAD.f[dM0M])[kbw  ] = c1o54  * conc * (c1o1 + c3o1 * (-vx1       - vx3) + c9o2 * (-vx1       - vx3) * (-vx1       - vx3) - cu_sq);
-        (distAD.f[dP0M])[kbe  ] = c1o54  * conc * (c1o1 + c3o1 * ( vx1       - vx3) + c9o2 * ( vx1       - vx3) * ( vx1       - vx3) - cu_sq);
-        (distAD.f[dM0P])[ktw  ] = c1o54  * conc * (c1o1 + c3o1 * (-vx1       + vx3) + c9o2 * (-vx1       + vx3) * (-vx1       + vx3) - cu_sq);
-        (distAD.f[d0PP])[ktn  ] = c1o54  * conc * (c1o1 + c3o1 * (       vx2 + vx3) + c9o2 * (       vx2 + vx3) * (       vx2 + vx3) - cu_sq);
-        (distAD.f[d0MM])[kbs  ] = c1o54  * conc * (c1o1 + c3o1 * (     - vx2 - vx3) + c9o2 * (     - vx2 - vx3) * (     - vx2 - vx3) - cu_sq);
-        (distAD.f[d0PM])[kbn  ] = c1o54  * conc * (c1o1 + c3o1 * (       vx2 - vx3) + c9o2 * (       vx2 - vx3) * (       vx2 - vx3) - cu_sq);
-        (distAD.f[d0MP])[kts  ] = c1o54  * conc * (c1o1 + c3o1 * (     - vx2 + vx3) + c9o2 * (     - vx2 + vx3) * (     - vx2 + vx3) - cu_sq);
-        (distAD.f[dPPP])[ktne ] = c1o216 * conc * (c1o1 + c3o1 * ( vx1 + vx2 + vx3) + c9o2 * ( vx1 + vx2 + vx3) * ( vx1 + vx2 + vx3) - cu_sq);
-        (distAD.f[dMMM])[kbsw ] = c1o216 * conc * (c1o1 + c3o1 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq);
-        (distAD.f[dPPM])[kbne ] = c1o216 * conc * (c1o1 + c3o1 * ( vx1 + vx2 - vx3) + c9o2 * ( vx1 + vx2 - vx3) * ( vx1 + vx2 - vx3) - cu_sq);
-        (distAD.f[dMMP])[ktsw ] = c1o216 * conc * (c1o1 + c3o1 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq);
-        (distAD.f[dPMP])[ktse ] = c1o216 * conc * (c1o1 + c3o1 * ( vx1 - vx2 + vx3) + c9o2 * ( vx1 - vx2 + vx3) * ( vx1 - vx2 + vx3) - cu_sq);
-        (distAD.f[dMPM])[kbnw ] = c1o216 * conc * (c1o1 + c3o1 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq);
-        (distAD.f[dPMM])[kbse ] = c1o216 * conc * (c1o1 + c3o1 * ( vx1 - vx2 - vx3) + c9o2 * ( vx1 - vx2 - vx3) * ( vx1 - vx2 - vx3) - cu_sq);
-        (distAD.f[dMPP])[ktnw ] = c1o216 * conc * (c1o1 + c3o1 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
-    }
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-// DEPRECATED (2022)
-
-// ////////////////////////////////////////////////////////////////////////////////
-// __global__ void InitAD27(unsigned int* neighborX,
-//                                        unsigned int* neighborY,
-//                                        unsigned int* neighborZ,
-//                                        unsigned int* geoD,
-//                                        real* Conc,
-//                                        real* ux,
-//                                        real* uy,
-//                                        real* uz,
-//                                        unsigned int size_Mat,
-//                                        real* DD27,
-//                                        bool EvenOrOdd)
-// {
-//    ////////////////////////////////////////////////////////////////////////////////
-//    const unsigned  x = threadIdx.x;  // Globaler x-Index 
-//    const unsigned  y = blockIdx.x;   // Globaler y-Index 
-//    const unsigned  z = blockIdx.y;   // Globaler z-Index 
-
-//    const unsigned nx = blockDim.x;
-//    const unsigned ny = gridDim.x;
-
-//    const unsigned k = nx*(ny*z + y) + x;
-//    //////////////////////////////////////////////////////////////////////////
-
-//    if(k<size_Mat)
-//    {
-//       ////////////////////////////////////////////////////////////////////////////////
-//       unsigned int BC;
-//       BC        =   geoD[k];
-
-//       if( BC != GEO_SOLID && BC != GEO_VOID)
-//       {
-//          Distributions27 D27;
-//          if (EvenOrOdd==true)
-//          {
-//             D27.f[dP00] = &DD27[dP00 * size_Mat];
-//             D27.f[dM00] = &DD27[dM00 * size_Mat];
-//             D27.f[d0P0] = &DD27[d0P0 * size_Mat];
-//             D27.f[d0M0] = &DD27[d0M0 * size_Mat];
-//             D27.f[d00P] = &DD27[d00P * size_Mat];
-//             D27.f[d00M] = &DD27[d00M * size_Mat];
-//             D27.f[dPP0] = &DD27[dPP0 * size_Mat];
-//             D27.f[dMM0] = &DD27[dMM0 * size_Mat];
-//             D27.f[dPM0] = &DD27[dPM0 * size_Mat];
-//             D27.f[dMP0] = &DD27[dMP0 * size_Mat];
-//             D27.f[dP0P] = &DD27[dP0P * size_Mat];
-//             D27.f[dM0M] = &DD27[dM0M * size_Mat];
-//             D27.f[dP0M] = &DD27[dP0M * size_Mat];
-//             D27.f[dM0P] = &DD27[dM0P * size_Mat];
-//             D27.f[d0PP] = &DD27[d0PP * size_Mat];
-//             D27.f[d0MM] = &DD27[d0MM * size_Mat];
-//             D27.f[d0PM] = &DD27[d0PM * size_Mat];
-//             D27.f[d0MP] = &DD27[d0MP * size_Mat];
-//             D27.f[d000] = &DD27[d000 * size_Mat];
-//             D27.f[dPPP] = &DD27[dPPP * size_Mat];
-//             D27.f[dMMP] = &DD27[dMMP * size_Mat];
-//             D27.f[dPMP] = &DD27[dPMP * size_Mat];
-//             D27.f[dMPP] = &DD27[dMPP * size_Mat];
-//             D27.f[dPPM] = &DD27[dPPM * size_Mat];
-//             D27.f[dMMM] = &DD27[dMMM * size_Mat];
-//             D27.f[dPMM] = &DD27[dPMM * size_Mat];
-//             D27.f[dMPM] = &DD27[dMPM * size_Mat];
-//          }
-//          else
-//          {
-//             D27.f[dM00] = &DD27[dP00 * size_Mat];
-//             D27.f[dP00] = &DD27[dM00 * size_Mat];
-//             D27.f[d0M0] = &DD27[d0P0 * size_Mat];
-//             D27.f[d0P0] = &DD27[d0M0 * size_Mat];
-//             D27.f[d00M] = &DD27[d00P * size_Mat];
-//             D27.f[d00P] = &DD27[d00M * size_Mat];
-//             D27.f[dMM0] = &DD27[dPP0 * size_Mat];
-//             D27.f[dPP0] = &DD27[dMM0 * size_Mat];
-//             D27.f[dMP0] = &DD27[dPM0 * size_Mat];
-//             D27.f[dPM0] = &DD27[dMP0 * size_Mat];
-//             D27.f[dM0M] = &DD27[dP0P * size_Mat];
-//             D27.f[dP0P] = &DD27[dM0M * size_Mat];
-//             D27.f[dM0P] = &DD27[dP0M * size_Mat];
-//             D27.f[dP0M] = &DD27[dM0P * size_Mat];
-//             D27.f[d0MM] = &DD27[d0PP * size_Mat];
-//             D27.f[d0PP] = &DD27[d0MM * size_Mat];
-//             D27.f[d0MP] = &DD27[d0PM * size_Mat];
-//             D27.f[d0PM] = &DD27[d0MP * size_Mat];
-//             D27.f[d000] = &DD27[d000 * size_Mat];
-//             D27.f[dMMM] = &DD27[dPPP * size_Mat];
-//             D27.f[dPPM] = &DD27[dMMP * size_Mat];
-//             D27.f[dMPM] = &DD27[dPMP * size_Mat];
-//             D27.f[dPMM] = &DD27[dMPP * size_Mat];
-//             D27.f[dMMP] = &DD27[dPPM * size_Mat];
-//             D27.f[dPPP] = &DD27[dMMM * size_Mat];
-//             D27.f[dMPP] = &DD27[dPMM * size_Mat];
-//             D27.f[dPMP] = &DD27[dMPM * size_Mat];
-//          }
-//          //////////////////////////////////////////////////////////////////////////
-//          real ConcD = Conc[k];
-//          real   vx1 = ux[k];
-//          real   vx2 = uy[k];
-//          real   vx3 = uz[k];
-//          //real lambdaD     = -three + sqrt(three);
-//          //real Diffusivity = c1o20;
-//          //real Lam         = -(c1o2+one/lambdaD);
-//          //real nue_d       = Lam/three;
-//          //real ae          = Diffusivity/nue_d - one;
-//          //real ux_sq       = vx1 * vx1;
-//          //real uy_sq       = vx2 * vx2;
-//          //real uz_sq       = vx3 * vx3;
-//          ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//          //D3Q7
-//          ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//          //index
-//          //unsigned int kzero= k;
-//          //unsigned int ke   = k;
-//          //unsigned int kw   = neighborX[k];
-//          //unsigned int kn   = k;
-//          //unsigned int ks   = neighborY[k];
-//          //unsigned int kt   = k;
-//          //unsigned int kb   = neighborZ[k];
-//          //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//          //(D7.f[0])[kzero] = ConcD*(c1o3*(ae*(-three))-(ux_sq+uy_sq+uz_sq));
-//          //(D7.f[1])[ke   ] = ConcD*(c1o6*(ae+one)+c1o2*(ux_sq)+vx1*c1o2);
-//          //(D7.f[2])[kw   ] = ConcD*(c1o6*(ae+one)+c1o2*(ux_sq)-vx1*c1o2);
-//          //(D7.f[3])[kn   ] = ConcD*(c1o6*(ae+one)+c1o2*(uy_sq)+vx2*c1o2);
-//          //(D7.f[4])[ks   ] = ConcD*(c1o6*(ae+one)+c1o2*(uy_sq)-vx2*c1o2);
-//          //(D7.f[5])[kt   ] = ConcD*(c1o6*(ae+one)+c1o2*(uz_sq)+vx3*c1o2);
-//          //(D7.f[6])[kb   ] = ConcD*(c1o6*(ae+one)+c1o2*(uz_sq)-vx3*c1o2);
-//          ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-
-
-//          ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//          //D3Q27
-//          ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//          //index
-//          unsigned int kzero= k;
-//          unsigned int ke   = k;
-//          unsigned int kw   = neighborX[k];
-//          unsigned int kn   = k;
-//          unsigned int ks   = neighborY[k];
-//          unsigned int kt   = k;
-//          unsigned int kb   = neighborZ[k];
-//          unsigned int ksw  = neighborY[kw];
-//          unsigned int kne  = k;
-//          unsigned int kse  = ks;
-//          unsigned int knw  = kw;
-//          unsigned int kbw  = neighborZ[kw];
-//          unsigned int kte  = k;
-//          unsigned int kbe  = kb;
-//          unsigned int ktw  = kw;
-//          unsigned int kbs  = neighborZ[ks];
-//          unsigned int ktn  = k;
-//          unsigned int kbn  = kb;
-//          unsigned int kts  = ks;
-//          unsigned int ktse = ks;
-//          unsigned int kbnw = kbw;
-//          unsigned int ktnw = kw;
-//          unsigned int kbse = kbs;
-//          unsigned int ktsw = ksw;
-//          unsigned int kbne = kb;
-//          unsigned int ktne = k;
-//          unsigned int kbsw = neighborZ[ksw];
-//          ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//          real cu_sq=c3o2*(vx1*vx1+vx2*vx2+vx3*vx3);
-
-//          (D27.f[d000])[kzero] =   c8o27* ConcD*(c1o1-cu_sq);
-//          (D27.f[dP00])[ke   ] =   c2o27* ConcD*(c1o1+c3o1*( vx1        )+c9o2*( vx1        )*( vx1        )-cu_sq);
-//          (D27.f[dM00])[kw   ] =   c2o27* ConcD*(c1o1+c3o1*(-vx1        )+c9o2*(-vx1        )*(-vx1        )-cu_sq);
-//          (D27.f[d0P0])[kn   ] =   c2o27* ConcD*(c1o1+c3o1*(    vx2     )+c9o2*(     vx2    )*(     vx2    )-cu_sq);
-//          (D27.f[d0M0])[ks   ] =   c2o27* ConcD*(c1o1+c3o1*(   -vx2     )+c9o2*(    -vx2    )*(    -vx2    )-cu_sq);
-//          (D27.f[d00P])[kt   ] =   c2o27* ConcD*(c1o1+c3o1*(         vx3)+c9o2*(         vx3)*(         vx3)-cu_sq);
-//          (D27.f[d00M])[kb   ] =   c2o27* ConcD*(c1o1+c3o1*(        -vx3)+c9o2*(        -vx3)*(        -vx3)-cu_sq);
-//          (D27.f[dPP0])[kne  ] =   c1o54* ConcD*(c1o1+c3o1*( vx1+vx2    )+c9o2*( vx1+vx2    )*( vx1+vx2    )-cu_sq);
-//          (D27.f[dMM0])[ksw  ] =   c1o54* ConcD*(c1o1+c3o1*(-vx1-vx2    )+c9o2*(-vx1-vx2    )*(-vx1-vx2    )-cu_sq);
-//          (D27.f[dPM0])[kse  ] =   c1o54* ConcD*(c1o1+c3o1*( vx1-vx2    )+c9o2*( vx1-vx2    )*( vx1-vx2    )-cu_sq);
-//          (D27.f[dMP0])[knw  ] =   c1o54* ConcD*(c1o1+c3o1*(-vx1+vx2    )+c9o2*(-vx1+vx2    )*(-vx1+vx2    )-cu_sq);
-//          (D27.f[dP0P])[kte  ] =   c1o54* ConcD*(c1o1+c3o1*( vx1    +vx3)+c9o2*( vx1    +vx3)*( vx1    +vx3)-cu_sq);
-//          (D27.f[dM0M])[kbw  ] =   c1o54* ConcD*(c1o1+c3o1*(-vx1    -vx3)+c9o2*(-vx1    -vx3)*(-vx1    -vx3)-cu_sq);
-//          (D27.f[dP0M])[kbe  ] =   c1o54* ConcD*(c1o1+c3o1*( vx1    -vx3)+c9o2*( vx1    -vx3)*( vx1    -vx3)-cu_sq);
-//          (D27.f[dM0P])[ktw  ] =   c1o54* ConcD*(c1o1+c3o1*(-vx1    +vx3)+c9o2*(-vx1    +vx3)*(-vx1    +vx3)-cu_sq);
-//          (D27.f[d0PP])[ktn  ] =   c1o54* ConcD*(c1o1+c3o1*(     vx2+vx3)+c9o2*(     vx2+vx3)*(     vx2+vx3)-cu_sq);
-//          (D27.f[d0MM])[kbs  ] =   c1o54* ConcD*(c1o1+c3o1*(    -vx2-vx3)+c9o2*(    -vx2-vx3)*(    -vx2-vx3)-cu_sq);
-//          (D27.f[d0PM])[kbn  ] =   c1o54* ConcD*(c1o1+c3o1*(     vx2-vx3)+c9o2*(     vx2-vx3)*(     vx2-vx3)-cu_sq);
-//          (D27.f[d0MP])[kts  ] =   c1o54* ConcD*(c1o1+c3o1*(    -vx2+vx3)+c9o2*(    -vx2+vx3)*(    -vx2+vx3)-cu_sq);
-//          (D27.f[dPPP])[ktne ] =   c1o216*ConcD*(c1o1+c3o1*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
-//          (D27.f[dMMM])[kbsw ] =   c1o216*ConcD*(c1o1+c3o1*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
-//          (D27.f[dPPM])[kbne ] =   c1o216*ConcD*(c1o1+c3o1*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
-//          (D27.f[dMMP])[ktsw ] =   c1o216*ConcD*(c1o1+c3o1*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
-//          (D27.f[dPMP])[ktse ] =   c1o216*ConcD*(c1o1+c3o1*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
-//          (D27.f[dMPM])[kbnw ] =   c1o216*ConcD*(c1o1+c3o1*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
-//          (D27.f[dPMM])[kbse ] =   c1o216*ConcD*(c1o1+c3o1*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
-//          (D27.f[dMPP])[ktnw ] =   c1o216*ConcD*(c1o1+c3o1*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
-//          ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
-//       }
-//    }
-// }
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-////////////////////////////////////////////////////////////////////////////////
-__global__ void InitAD7( unsigned int* neighborX,
-                                    unsigned int* neighborY,
-                                    unsigned int* neighborZ,
-                                    unsigned int* geoD,
-                                    real* Conc,
-                                    real* ux,
-                                    real* uy,
-                                    real* uz,
-                                    unsigned long long numberOfLBnodes,
-                                    real* DD7,
-                                    bool EvenOrOdd)
-{
-   ////////////////////////////////////////////////////////////////////////////////
-   const unsigned  x = threadIdx.x;  // Globaler x-Index 
-   const unsigned  y = blockIdx.x;   // Globaler y-Index 
-   const unsigned  z = blockIdx.y;   // Globaler z-Index 
-
-   const unsigned nx = blockDim.x;
-   const unsigned ny = gridDim.x;
-
-   const unsigned k = nx*(ny*z + y) + x;
-   //////////////////////////////////////////////////////////////////////////
-
-   if(k<numberOfLBnodes)
-   {
-      ////////////////////////////////////////////////////////////////////////////////
-      unsigned int BC;
-      BC        =   geoD[k];
-
-      if( BC != GEO_SOLID && BC != GEO_VOID)
-      {
-         Distributions7 D7;
-         if (EvenOrOdd==true)
-         {
-            D7.f[0] = &DD7[0*numberOfLBnodes];
-            D7.f[1] = &DD7[1*numberOfLBnodes];
-            D7.f[2] = &DD7[2*numberOfLBnodes];
-            D7.f[3] = &DD7[3*numberOfLBnodes];
-            D7.f[4] = &DD7[4*numberOfLBnodes];
-            D7.f[5] = &DD7[5*numberOfLBnodes];
-            D7.f[6] = &DD7[6*numberOfLBnodes];
-         }
-         else
-         {
-            D7.f[0] = &DD7[0*numberOfLBnodes];
-            D7.f[2] = &DD7[1*numberOfLBnodes];
-            D7.f[1] = &DD7[2*numberOfLBnodes];
-            D7.f[4] = &DD7[3*numberOfLBnodes];
-            D7.f[3] = &DD7[4*numberOfLBnodes];
-            D7.f[6] = &DD7[5*numberOfLBnodes];
-            D7.f[5] = &DD7[6*numberOfLBnodes];
-         }
-         //////////////////////////////////////////////////////////////////////////
-         real ConcD = Conc[k];
-         real   vx1 = ux[k];
-         real   vx2 = uy[k];
-         real   vx3 = uz[k];
-         real lambdaD     = -c3o1 + sqrt(c3o1);
-         real Diffusivity = c1o20;
-         real Lam         = -(c1o2+c1o1/lambdaD);
-         real nue_d       = Lam/c3o1;
-         real ae          = Diffusivity/nue_d - c1o1;
-         real ux_sq       = vx1 * vx1;
-         real uy_sq       = vx2 * vx2;
-         real uz_sq       = vx3 * vx3;
-         //////////////////////////////////////////////////////////////////////////
-         //index
-         //////////////////////////////////////////////////////////////////////////
-         unsigned int kzero= k;
-         unsigned int ke   = k;
-         unsigned int kw   = neighborX[k];
-         unsigned int kn   = k;
-         unsigned int ks   = neighborY[k];
-         unsigned int kt   = k;
-         unsigned int kb   = neighborZ[k];
-         //////////////////////////////////////////////////////////////////////////
-
-         (D7.f[0])[kzero] = ConcD*(c1o3*(ae*(-c3o1))-(ux_sq+uy_sq+uz_sq));
-         (D7.f[1])[ke   ] = ConcD*(c1o6*(ae+c1o1)+c1o2*(ux_sq)+vx1*c1o2);
-         (D7.f[2])[kw   ] = ConcD*(c1o6*(ae+c1o1)+c1o2*(ux_sq)-vx1*c1o2);
-         (D7.f[3])[kn   ] = ConcD*(c1o6*(ae+c1o1)+c1o2*(uy_sq)+vx2*c1o2);
-         (D7.f[4])[ks   ] = ConcD*(c1o6*(ae+c1o1)+c1o2*(uy_sq)-vx2*c1o2);
-         (D7.f[5])[kt   ] = ConcD*(c1o6*(ae+c1o1)+c1o2*(uz_sq)+vx3*c1o2);
-         (D7.f[6])[kb   ] = ConcD*(c1o6*(ae+c1o1)+c1o2*(uz_sq)-vx3*c1o2);
-      }
-   }
-}
\ No newline at end of file
diff --git a/src/gpu/core/GPU/LBMKernel.cu b/src/gpu/core/GPU/LBMKernel.cu
index 67aa3aead..403ed748f 100644
--- a/src/gpu/core/GPU/LBMKernel.cu
+++ b/src/gpu/core/GPU/LBMKernel.cu
@@ -18,111 +18,6 @@
 
 #include "Parameter/Parameter.h"
 //////////////////////////////////////////////////////////////////////////
-void Init27(
-    int myid,
-    int numprocs,
-    real u0,
-    unsigned int* geoD,
-    unsigned int* neighborX,
-    unsigned int* neighborY,
-    unsigned int* neighborZ,
-    real* vParab,
-    unsigned long long numberOfLBnodes,
-    unsigned int grid_nx,
-    unsigned int grid_ny,
-    unsigned int grid_nz,
-    real* DD,
-    int level,
-    int maxlevel)
-{
-    dim3 threads       ( grid_nx, 1, 1 );
-    dim3 grid          ( grid_ny, grid_nz );
-
-    LBInit27<<< grid, threads >>> (
-        myid,
-        numprocs,
-        u0,
-        geoD,
-        neighborX,
-        neighborY,
-        neighborZ,
-        vParab,
-        numberOfLBnodes,
-        grid_nx,
-        grid_ny,
-        grid_nz,
-        DD,
-        level,
-        maxlevel);
-    getLastCudaError("LBInit27 execution failed");
-}
-//////////////////////////////////////////////////////////////////////////
-void InitNonEqPartSP27(
-    unsigned int numberOfThreads,
-    unsigned int* neighborX,
-    unsigned int* neighborY,
-    unsigned int* neighborZ,
-    unsigned int* neighborWSB,
-    unsigned int* geoD,
-    real* rho,
-    real* ux,
-    real* uy,
-    real* uz,
-    unsigned long long numberOfLBnodes,
-    real* DD,
-    real omega,
-    bool EvenOrOdd)
-{
-    vf::cuda::CudaGrid grid = vf::cuda::CudaGrid(numberOfThreads, numberOfLBnodes);
-
-    LBInitNonEqPartSP27<<< grid.grid, grid.threads >>>(
-        neighborX,
-        neighborY,
-        neighborZ,
-        neighborWSB,
-        geoD,
-        rho,
-        ux,
-        uy,
-        uz,
-        numberOfLBnodes,
-        DD,
-        omega,
-        EvenOrOdd);
-    getLastCudaError("LBInitNonEqPartSP27 execution failed");
-}
-//////////////////////////////////////////////////////////////////////////
-void InitADDev27(
-    unsigned int numberOfThreads,
-    unsigned int* neighborX,
-    unsigned int* neighborY,
-    unsigned int* neighborZ,
-    unsigned int* geoD,
-    real* Conc,
-    real* ux,
-    real* uy,
-    real* uz,
-    unsigned long long numberOfLBnodes,
-    real* DD27,
-    bool EvenOrOdd)
-{
-    vf::cuda::CudaGrid grid = vf::cuda::CudaGrid(numberOfThreads, numberOfLBnodes);
-
-    InitAD27<<< grid.grid, grid.threads >>>(
-        neighborX,
-        neighborY,
-        neighborZ,
-        geoD,
-        Conc,
-        ux,
-        uy,
-        uz,
-        numberOfLBnodes,
-        DD27,
-        EvenOrOdd);
-    getLastCudaError("InitAD27 execution failed");
-}
-//////////////////////////////////////////////////////////////////////////
 void CalcMac27(
     real* vxD,
     real* vyD,
diff --git a/src/gpu/core/Init/InitLattice.cpp b/src/gpu/core/Init/InitLattice.cpp
index ba4a5cae5..f8e00006e 100644
--- a/src/gpu/core/Init/InitLattice.cpp
+++ b/src/gpu/core/Init/InitLattice.cpp
@@ -39,7 +39,7 @@
 #include "Temperature/FindTemperature.h"
 
 
-void initLattice(SPtr<Parameter> para, SPtr<PreProcessor> preProcessor, SPtr<CudaMemoryManager> cudaMemoryManager)
+void initLattice(SPtr<Parameter> para, SPtr<PreProcessor> preProcessor, SPtr<PreProcessor> preProcessorAD, SPtr<CudaMemoryManager> cudaMemoryManager)
 {
     for (int lev = para->getFine(); lev >= para->getCoarse(); lev--) {
         preProcessor->init(para, lev);
@@ -84,7 +84,8 @@ void initLattice(SPtr<Parameter> para, SPtr<PreProcessor> preProcessor, SPtr<Cud
             for (size_t index = 0; index < para->getParH(lev)->numberOfNodes; index++) {
                 para->getParH(lev)->concentration[index] = para->getTemperatureInit();
             }
-            initTemperatur(para.get(), cudaMemoryManager.get(), lev);
+
+            preProcessorAD->init(para, lev);
         }
     }
 }
diff --git a/src/gpu/core/Init/InitLattice.h b/src/gpu/core/Init/InitLattice.h
index 8de22f321..2e51f911c 100644
--- a/src/gpu/core/Init/InitLattice.h
+++ b/src/gpu/core/Init/InitLattice.h
@@ -39,6 +39,6 @@ class Parameter;
 class PreProcessor;
 class CudaMemoryManager;
 
-void initLattice(SPtr<Parameter> para, SPtr<PreProcessor> preProcessor, SPtr<CudaMemoryManager> cudaMemoryManager);
+void initLattice(SPtr<Parameter> para, SPtr<PreProcessor> preProcessor, SPtr<PreProcessor> preProcessorAD, SPtr<CudaMemoryManager> cudaMemoryManager);
 
 #endif
diff --git a/src/gpu/core/KernelManager/ADKernelManager.cpp b/src/gpu/core/KernelManager/ADKernelManager.cpp
index ac504a839..7b867dc84 100644
--- a/src/gpu/core/KernelManager/ADKernelManager.cpp
+++ b/src/gpu/core/KernelManager/ADKernelManager.cpp
@@ -38,56 +38,6 @@
 
 ADKernelManager::ADKernelManager(SPtr<Parameter> parameter, std::vector<SPtr<AdvectionDiffusionKernel>>& adkernels): para(parameter), adkernels(adkernels){}
 
-void ADKernelManager::initAD(const int level) const
-{
-    //////////////////////////////////////////////////////////////////////////
-    // calculation of omega for diffusivity
-    para->getParD(level)->omegaDiffusivity = (real)2.0 / ((real)6.0 * para->getParD(level)->diffusivity + (real)1.0);
-    //////////////////////////////////////////////////////////////////////////
-    para->getParD(level)->isEvenTimestep = true;
-    //////////////////////////////////////////////////////////////////////////
-    initAdvectionDiffusion(
-        para->getParD(level)->numberofthreads, 
-        para->getParD(level)->neighborX, 
-        para->getParD(level)->neighborY,
-        para->getParD(level)->neighborZ, 
-        para->getParD(level)->typeOfGridNode, 
-        para->getParD(level)->concentration,
-        para->getParD(level)->velocityX, 
-        para->getParD(level)->velocityY, 
-        para->getParD(level)->velocityZ,
-        para->getParD(level)->numberOfNodes, 
-        para->getParD(level)->distributionsAD.f[0],
-        para->getParD(level)->isEvenTimestep);
-    //////////////////////////////////////////////////////////////////////////
-    para->getParD(level)->isEvenTimestep = false;
-    //////////////////////////////////////////////////////////////////////////
-    initAdvectionDiffusion(
-        para->getParD(level)->numberofthreads, 
-        para->getParD(level)->neighborX, 
-        para->getParD(level)->neighborY,
-        para->getParD(level)->neighborZ, 
-        para->getParD(level)->typeOfGridNode, 
-        para->getParD(level)->concentration,
-        para->getParD(level)->velocityX, 
-        para->getParD(level)->velocityY, 
-        para->getParD(level)->velocityZ,
-        para->getParD(level)->numberOfNodes, 
-        para->getParD(level)->distributionsAD.f[0],
-        para->getParD(level)->isEvenTimestep);
-    //////////////////////////////////////////////////////////////////////////
-    CalcConcentration27(
-        para->getParD(level)->numberofthreads,
-        para->getParD(level)->concentration,
-        para->getParD(level)->typeOfGridNode,
-        para->getParD(level)->neighborX,
-        para->getParD(level)->neighborY,
-        para->getParD(level)->neighborZ,
-        para->getParD(level)->numberOfNodes,
-        para->getParD(level)->distributionsAD.f[0],
-        para->getParD(level)->isEvenTimestep);
-}
-
 ////////////////////////////////////////////////////////////////////////////////
 void ADKernelManager::setInitialNodeValuesAD(const int level, SPtr<CudaMemoryManager> cudaMemoryManager) const
 {
diff --git a/src/gpu/core/KernelManager/ADKernelManager.h b/src/gpu/core/KernelManager/ADKernelManager.h
index 46d818b04..86a59d583 100644
--- a/src/gpu/core/KernelManager/ADKernelManager.h
+++ b/src/gpu/core/KernelManager/ADKernelManager.h
@@ -54,9 +54,6 @@ public:
     //! \param parameter shared pointer to instance of class Parameter
     ADKernelManager(SPtr<Parameter> parameter, std::vector<SPtr<AdvectionDiffusionKernel>>& adkernels);
 
-    //! \brief initialize the advection diffusion distributions
-    void initAD(const int level) const;
-
     //! \brief set initial concentration values at all nodes
     //! \param cudaMemoryManager instance of class CudaMemoryManager
     void setInitialNodeValuesAD(const int level, SPtr<CudaMemoryManager> cudaMemoryManager) const;
diff --git a/src/gpu/core/LBM/Simulation.cpp b/src/gpu/core/LBM/Simulation.cpp
index 2a7787155..3566c26f8 100644
--- a/src/gpu/core/LBM/Simulation.cpp
+++ b/src/gpu/core/LBM/Simulation.cpp
@@ -156,6 +156,8 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     if (para->getDiffOn()) {
         VF_LOG_INFO("make AD Kernels");
         adKernels = kernelFactory->makeAdvDifKernels(para);
+        std::vector<PreProcessorType> preProADTypes = adKernels.at(0)->getPreProcessorTypes();
+        preProcessorAD = preProcessorFactory->makePreProcessor(preProADTypes, para);
     }
 
     //////////////////////////////////////////////////////////////////////////
@@ -270,7 +272,7 @@ void Simulation::init(GridProvider &gridProvider, BoundaryConditionFactory *bcFa
     // VF_LOG_INFO("done.");
 
     VF_LOG_INFO("init lattice...");
-    initLattice(para, preProcessor, cudaMemoryManager);
+    initLattice(para, preProcessor, preProcessorAD, cudaMemoryManager);
     VF_LOG_INFO("done");
 
     // VF_LOG_INFO("set geo for Q...\n");
diff --git a/src/gpu/core/LBM/Simulation.h b/src/gpu/core/LBM/Simulation.h
index e3bd7cf7b..4a9c89ba0 100644
--- a/src/gpu/core/LBM/Simulation.h
+++ b/src/gpu/core/LBM/Simulation.h
@@ -82,6 +82,7 @@ private:
     std::vector < SPtr< Kernel>> kernels;
     std::vector < SPtr< AdvectionDiffusionKernel>> adKernels;
     std::shared_ptr<PreProcessor> preProcessor;
+    std::shared_ptr<PreProcessor> preProcessorAD;
     SPtr<TurbulenceModelFactory> tmFactory;
 
     SPtr<RestartObject> restart_object;
-- 
GitLab