diff --git a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
index a10b5e43ed5b8c43ee6fd9c6b24fef3614761ff5..f71c08bb49f2319bf95001df8441500348a939c6 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
@@ -250,7 +250,346 @@ extern "C" __global__ void LBCalcMacCompSP27(
                     rho[k]-(velocityX[k] * velocityX[k] + velocityY[k] * velocityY[k] + velocityZ[k] * velocityZ[k]) * (c1o1+rho[k])) * (c1o1 / OxxPyyPzz - c1o2) +rho[k];
    }
 }
+////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+extern "C" __global__ void LBCalcMacADCompSP27( 
+	real* velocityX,
+	real* velocityY,
+	real* velocityZ,
+	real* rho,
+	real* pressure,
+	uint* typeOfGridNode,
+	uint* neighborX,
+	uint* neighborY,
+	uint* neighborZ,
+	uint size_Mat,
+	real* distributions,
+    real* distributionsAD,
+    real* forces,
+	bool isEvenTimestep)
+{
+   //////////////////////////////////////////////////////////////////////////
+   //! The velocity boundary condition is executed in the following steps
+   //!
+   ////////////////////////////////////////////////////////////////////////////////
+   //! - Get node index coordinates from thredIdx, 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 < size_Mat) && (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 dist;
+      if (isEvenTimestep)
+      {
+         dist.f[dirE   ] = &distributions[dirE   *size_Mat];
+         dist.f[dirW   ] = &distributions[dirW   *size_Mat];
+         dist.f[dirN   ] = &distributions[dirN   *size_Mat];
+         dist.f[dirS   ] = &distributions[dirS   *size_Mat];
+         dist.f[dirT   ] = &distributions[dirT   *size_Mat];
+         dist.f[dirB   ] = &distributions[dirB   *size_Mat];
+         dist.f[dirNE  ] = &distributions[dirNE  *size_Mat];
+         dist.f[dirSW  ] = &distributions[dirSW  *size_Mat];
+         dist.f[dirSE  ] = &distributions[dirSE  *size_Mat];
+         dist.f[dirNW  ] = &distributions[dirNW  *size_Mat];
+         dist.f[dirTE  ] = &distributions[dirTE  *size_Mat];
+         dist.f[dirBW  ] = &distributions[dirBW  *size_Mat];
+         dist.f[dirBE  ] = &distributions[dirBE  *size_Mat];
+         dist.f[dirTW  ] = &distributions[dirTW  *size_Mat];
+         dist.f[dirTN  ] = &distributions[dirTN  *size_Mat];
+         dist.f[dirBS  ] = &distributions[dirBS  *size_Mat];
+         dist.f[dirBN  ] = &distributions[dirBN  *size_Mat];
+         dist.f[dirTS  ] = &distributions[dirTS  *size_Mat];
+         dist.f[dirREST] = &distributions[dirREST*size_Mat];
+         dist.f[dirTNE ] = &distributions[dirTNE *size_Mat];
+         dist.f[dirTSW ] = &distributions[dirTSW *size_Mat];
+         dist.f[dirTSE ] = &distributions[dirTSE *size_Mat];
+         dist.f[dirTNW ] = &distributions[dirTNW *size_Mat];
+         dist.f[dirBNE ] = &distributions[dirBNE *size_Mat];
+         dist.f[dirBSW ] = &distributions[dirBSW *size_Mat];
+         dist.f[dirBSE ] = &distributions[dirBSE *size_Mat];
+         dist.f[dirBNW ] = &distributions[dirBNW *size_Mat];
+      } 
+      else
+      {
+         dist.f[dirW   ] = &distributions[dirE   *size_Mat];
+         dist.f[dirE   ] = &distributions[dirW   *size_Mat];
+         dist.f[dirS   ] = &distributions[dirN   *size_Mat];
+         dist.f[dirN   ] = &distributions[dirS   *size_Mat];
+         dist.f[dirB   ] = &distributions[dirT   *size_Mat];
+         dist.f[dirT   ] = &distributions[dirB   *size_Mat];
+         dist.f[dirSW  ] = &distributions[dirNE  *size_Mat];
+         dist.f[dirNE  ] = &distributions[dirSW  *size_Mat];
+         dist.f[dirNW  ] = &distributions[dirSE  *size_Mat];
+         dist.f[dirSE  ] = &distributions[dirNW  *size_Mat];
+         dist.f[dirBW  ] = &distributions[dirTE  *size_Mat];
+         dist.f[dirTE  ] = &distributions[dirBW  *size_Mat];
+         dist.f[dirTW  ] = &distributions[dirBE  *size_Mat];
+         dist.f[dirBE  ] = &distributions[dirTW  *size_Mat];
+         dist.f[dirBS  ] = &distributions[dirTN  *size_Mat];
+         dist.f[dirTN  ] = &distributions[dirBS  *size_Mat];
+         dist.f[dirTS  ] = &distributions[dirBN  *size_Mat];
+         dist.f[dirBN  ] = &distributions[dirTS  *size_Mat];
+         dist.f[dirREST] = &distributions[dirREST*size_Mat];
+         dist.f[dirTNE ] = &distributions[dirBSW *size_Mat];
+         dist.f[dirTSW ] = &distributions[dirBNE *size_Mat];
+         dist.f[dirTSE ] = &distributions[dirBNW *size_Mat];
+         dist.f[dirTNW ] = &distributions[dirBSE *size_Mat];
+         dist.f[dirBNE ] = &distributions[dirTSW *size_Mat];
+         dist.f[dirBSW ] = &distributions[dirTNE *size_Mat];
+         dist.f[dirBSE ] = &distributions[dirTNW *size_Mat];
+         dist.f[dirBNW ] = &distributions[dirTSE *size_Mat];
+      }
+      ////////////////////////////////////////////////////////////////////////////////
+       Distributions27 distAD;
+       if (isEvenTimestep)
+       {
+           distAD.f[dirE   ] = &distributionsAD[dirE   *size_Mat];
+           distAD.f[dirW   ] = &distributionsAD[dirW   *size_Mat];
+           distAD.f[dirN   ] = &distributionsAD[dirN   *size_Mat];
+           distAD.f[dirS   ] = &distributionsAD[dirS   *size_Mat];
+           distAD.f[dirT   ] = &distributionsAD[dirT   *size_Mat];
+           distAD.f[dirB   ] = &distributionsAD[dirB   *size_Mat];
+           distAD.f[dirNE  ] = &distributionsAD[dirNE  *size_Mat];
+           distAD.f[dirSW  ] = &distributionsAD[dirSW  *size_Mat];
+           distAD.f[dirSE  ] = &distributionsAD[dirSE  *size_Mat];
+           distAD.f[dirNW  ] = &distributionsAD[dirNW  *size_Mat];
+           distAD.f[dirTE  ] = &distributionsAD[dirTE  *size_Mat];
+           distAD.f[dirBW  ] = &distributionsAD[dirBW  *size_Mat];
+           distAD.f[dirBE  ] = &distributionsAD[dirBE  *size_Mat];
+           distAD.f[dirTW  ] = &distributionsAD[dirTW  *size_Mat];
+           distAD.f[dirTN  ] = &distributionsAD[dirTN  *size_Mat];
+           distAD.f[dirBS  ] = &distributionsAD[dirBS  *size_Mat];
+           distAD.f[dirBN  ] = &distributionsAD[dirBN  *size_Mat];
+           distAD.f[dirTS  ] = &distributionsAD[dirTS  *size_Mat];
+           distAD.f[dirREST] = &distributionsAD[dirREST*size_Mat];
+           distAD.f[dirTNE ] = &distributionsAD[dirTNE *size_Mat];
+           distAD.f[dirTSW ] = &distributionsAD[dirTSW *size_Mat];
+           distAD.f[dirTSE ] = &distributionsAD[dirTSE *size_Mat];
+           distAD.f[dirTNW ] = &distributionsAD[dirTNW *size_Mat];
+           distAD.f[dirBNE ] = &distributionsAD[dirBNE *size_Mat];
+           distAD.f[dirBSW ] = &distributionsAD[dirBSW *size_Mat];
+           distAD.f[dirBSE ] = &distributionsAD[dirBSE *size_Mat];
+           distAD.f[dirBNW ] = &distributionsAD[dirBNW *size_Mat];
+       }
+       else
+       {
+           distAD.f[dirW   ] = &distributionsAD[dirE   *size_Mat];
+           distAD.f[dirE   ] = &distributionsAD[dirW   *size_Mat];
+           distAD.f[dirS   ] = &distributionsAD[dirN   *size_Mat];
+           distAD.f[dirN   ] = &distributionsAD[dirS   *size_Mat];
+           distAD.f[dirB   ] = &distributionsAD[dirT   *size_Mat];
+           distAD.f[dirT   ] = &distributionsAD[dirB   *size_Mat];
+           distAD.f[dirSW  ] = &distributionsAD[dirNE  *size_Mat];
+           distAD.f[dirNE  ] = &distributionsAD[dirSW  *size_Mat];
+           distAD.f[dirNW  ] = &distributionsAD[dirSE  *size_Mat];
+           distAD.f[dirSE  ] = &distributionsAD[dirNW  *size_Mat];
+           distAD.f[dirBW  ] = &distributionsAD[dirTE  *size_Mat];
+           distAD.f[dirTE  ] = &distributionsAD[dirBW  *size_Mat];
+           distAD.f[dirTW  ] = &distributionsAD[dirBE  *size_Mat];
+           distAD.f[dirBE  ] = &distributionsAD[dirTW  *size_Mat];
+           distAD.f[dirBS  ] = &distributionsAD[dirTN  *size_Mat];
+           distAD.f[dirTN  ] = &distributionsAD[dirBS  *size_Mat];
+           distAD.f[dirTS  ] = &distributionsAD[dirBN  *size_Mat];
+           distAD.f[dirBN  ] = &distributionsAD[dirTS  *size_Mat];
+           distAD.f[dirREST] = &distributionsAD[dirREST*size_Mat];
+           distAD.f[dirBSW ] = &distributionsAD[dirTNE *size_Mat];
+           distAD.f[dirBNE ] = &distributionsAD[dirTSW *size_Mat];
+           distAD.f[dirBNW ] = &distributionsAD[dirTSE *size_Mat];
+           distAD.f[dirBSE ] = &distributionsAD[dirTNW *size_Mat];
+           distAD.f[dirTSW ] = &distributionsAD[dirBNE *size_Mat];
+           distAD.f[dirTNE ] = &distributionsAD[dirBSW *size_Mat];
+           distAD.f[dirTNW ] = &distributionsAD[dirBSE *size_Mat];
+           distAD.f[dirTSE ] = &distributionsAD[dirBNW *size_Mat];
+       }
+	  ////////////////////////////////////////////////////////////////////////////////
+	  //! - Set neighbor indices (necessary for indirect addressing)  
+	  //!
+	  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];
+	  ////////////////////////////////////////////////////////////////////////////////
+	  //! - Set local distributions  
+	  //!
+	  real mfcbb = (dist.f[dirE   ])[k   ];
+	  real mfabb = (dist.f[dirW   ])[kw  ];
+	  real mfbcb = (dist.f[dirN   ])[k   ];
+	  real mfbab = (dist.f[dirS   ])[ks  ];
+	  real mfbbc = (dist.f[dirT   ])[k   ];
+	  real mfbba = (dist.f[dirB   ])[kb  ];
+	  real mfccb = (dist.f[dirNE  ])[k   ];
+	  real mfaab = (dist.f[dirSW  ])[ksw ];
+	  real mfcab = (dist.f[dirSE  ])[ks  ];
+	  real mfacb = (dist.f[dirNW  ])[kw  ];
+	  real mfcbc = (dist.f[dirTE  ])[k   ];
+	  real mfaba = (dist.f[dirBW  ])[kbw ];
+	  real mfcba = (dist.f[dirBE  ])[kb  ];
+	  real mfabc = (dist.f[dirTW  ])[kw  ];
+	  real mfbcc = (dist.f[dirTN  ])[k   ];
+	  real mfbaa = (dist.f[dirBS  ])[kbs ];
+	  real mfbca = (dist.f[dirBN  ])[kb  ];
+	  real mfbac = (dist.f[dirTS  ])[ks  ];
+	  real mfbbb = (dist.f[dirREST])[k   ];
+	  real mfccc = (dist.f[dirTNE ])[k   ];
+	  real mfaac = (dist.f[dirTSW ])[ksw ];
+	  real mfcac = (dist.f[dirTSE ])[ks  ];
+	  real mfacc = (dist.f[dirTNW ])[kw  ];
+	  real mfcca = (dist.f[dirBNE ])[kb  ];
+	  real mfaaa = (dist.f[dirBSW ])[kbsw];
+	  real mfcaa = (dist.f[dirBSE ])[kbs ];
+	  real mfaca = (dist.f[dirBNW ])[kbw ];
+	  ////////////////////////////////////////////////////////////////////////////////////
+       //! - Set local distributions Advection Diffusion
+       //!
+       real fcbb = (distAD.f[dirE   ])[k];
+       real fabb = (distAD.f[dirW   ])[kw];
+       real fbcb = (distAD.f[dirN   ])[k];
+       real fbab = (distAD.f[dirS   ])[ks];
+       real fbbc = (distAD.f[dirT   ])[k];
+       real fbba = (distAD.f[dirB   ])[kb];
+       real fccb = (distAD.f[dirNE  ])[k];
+       real faab = (distAD.f[dirSW  ])[ksw];
+       real fcab = (distAD.f[dirSE  ])[ks];
+       real facb = (distAD.f[dirNW  ])[kw];
+       real fcbc = (distAD.f[dirTE  ])[k];
+       real faba = (distAD.f[dirBW  ])[kbw];
+       real fcba = (distAD.f[dirBE  ])[kb];
+       real fabc = (distAD.f[dirTW  ])[kw];
+       real fbcc = (distAD.f[dirTN  ])[k];
+       real fbaa = (distAD.f[dirBS  ])[kbs];
+       real fbca = (distAD.f[dirBN  ])[kb];
+       real fbac = (distAD.f[dirTS  ])[ks];
+       real fbbb = (distAD.f[dirREST])[k];
+       real fccc = (distAD.f[dirTNE ])[k];
+       real faac = (distAD.f[dirTSW ])[ksw];
+       real fcac = (distAD.f[dirTSE ])[ks];
+       real facc = (distAD.f[dirTNW ])[kw];
+       real fcca = (distAD.f[dirBNE ])[kb];
+       real faaa = (distAD.f[dirBSW ])[kbsw];
+       real fcaa = (distAD.f[dirBSE ])[kbs];
+       real faca = (distAD.f[dirBNW ])[kbw];
+	  ////////////////////////////////////////////////////////////////////////////////
+	  //! - Set pressure, density and velocity to \f$ 1.0 \f$  
+	  //!
+	  pressure[k]  = c0o1;
+	  rho[k]       = c0o1;
+	  velocityX[k] = c0o1;
+	  velocityY[k] = c0o1;
+	  velocityZ[k] = c0o1;
 
+      //////////////////////////////////////////////////////////////////////////
+	  //! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq. (J1)-(J3) \ref
+	  //! <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>
+	  //!
+	  real drho = ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
+      	(((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
+      	((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
+      
+      real rhoDev = c1o1 + drho;
+      ////////////////////////////////////////////////////////////////////////////////////
+      real vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+      	(((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+      	(mfcbb - mfabb)) / rhoDev;
+      real vvy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+      	(((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+      	(mfbcb - mfbab)) / rhoDev;
+      real vvz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+      	(((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+      	(mfbbc - mfbba)) / rhoDev;
+       ////////////////////////////////////////////////////////////////////////////////////
+       // scalar component (ADE)
+       real rhoScalar =
+               ((((fccc + faaa) + (faca + fcac)) + ((facc + fcaa) + (faac + fcca))) +
+                (((fbac + fbca) + (fbaa + fbcc)) + ((fabc + fcba) + (faba + fcbc)) + ((facb + fcab) + (faab + fccb))) +
+                ((fabb + fcbb) + (fbab + fbcb) + (fbba + fbbc))) + fbbb;
+       ////////////////////////////////////////////////////////////////////////////////////
+       //! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref
+       //! <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>
+       //!
+       real fx = forces[0];
+       real fy = forces[1];
+       real fz = -rhoScalar*forces[2];
+      
+      rho[k]        = drho;
+      velocityX[k]  = vvx+fx*c1o2;
+      velocityY[k]  = vvy+fy*c1o2;
+      velocityZ[k]  = vvz+fz*c1o2;
+      
+	  //////////////////////////////////////////////////////////////////////////
+	  //! - Calculate pressure
+	  //!
+	  real OxxPyyPzz = c1o1;
+      pressure[k] =((dist.f[dirE   ])[ke  ]+ (dist.f[dirW   ])[kw  ]+ 
+                    (dist.f[dirN   ])[kn  ]+ (dist.f[dirS   ])[ks  ]+
+                    (dist.f[dirT   ])[kt  ]+ (dist.f[dirB   ])[kb  ]+
+                    c2o1*(
+                    (dist.f[dirNE  ])[kne ]+ (dist.f[dirSW  ])[ksw ]+
+                    (dist.f[dirSE  ])[kse ]+ (dist.f[dirNW  ])[knw ]+
+                    (dist.f[dirTE  ])[kte ]+ (dist.f[dirBW  ])[kbw ]+
+                    (dist.f[dirBE  ])[kbe ]+ (dist.f[dirTW  ])[ktw ]+
+                    (dist.f[dirTN  ])[ktn ]+ (dist.f[dirBS  ])[kbs ]+
+                    (dist.f[dirBN  ])[kbn ]+ (dist.f[dirTS  ])[kts ])+
+                    c3o1*(
+                    (dist.f[dirTNE ])[ktne]+ (dist.f[dirTSW ])[ktsw]+ 
+                    (dist.f[dirTSE ])[ktse]+ (dist.f[dirTNW ])[ktnw]+ 
+                    (dist.f[dirBNE ])[kbne]+ (dist.f[dirBSW ])[kbsw]+ 
+                    (dist.f[dirBSE ])[kbse]+ (dist.f[dirBNW ])[kbnw])-
+                    rho[k]-(velocityX[k] * velocityX[k] + velocityY[k] * velocityY[k] + velocityZ[k] * velocityZ[k]) * (c1o1+rho[k])) * (c1o1 / OxxPyyPzz - c1o2) +rho[k];
+   }
+}
 
 
 
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CudaKernelManager.cpp b/src/gpu/VirtualFluids_GPU/GPU/CudaKernelManager.cpp
index 6705edffc1ef7e777125e71f94b07c92aeb3ac8a..451174fa70031aca85c806f72f79e378d667220d 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CudaKernelManager.cpp
+++ b/src/gpu/VirtualFluids_GPU/GPU/CudaKernelManager.cpp
@@ -91,20 +91,39 @@ void CudaKernelManager::runVelocityBCKernel(SPtr<Parameter> para)
 
 void CudaKernelManager::calculateMacroscopicValues(SPtr<Parameter> para)
 {
-	CalcMacCompSP27(
-		para->getParD()->velocityX,
-		para->getParD()->velocityY,
-		para->getParD()->velocityZ,
-		para->getParD()->rho,
-		para->getParD()->pressure,
-		para->getParD()->typeOfGridNode,
-		para->getParD()->neighborX,
-		para->getParD()->neighborY,
-		para->getParD()->neighborZ,
-		para->getParD()->numberOfNodes,
-		para->getParD()->numberofthreads,
-		para->getParD()->distributions.f[0],
-		para->getParD()->isEvenTimestep);
+    if (para->getIsADcalculationOn()) {
+		CalcMacADCompSP27(
+			para->getParD()->velocityX,
+			para->getParD()->velocityY,
+			para->getParD()->velocityZ,
+			para->getParD()->rho,
+			para->getParD()->pressure,
+			para->getParD()->typeOfGridNode,
+			para->getParD()->neighborX,
+			para->getParD()->neighborY,
+			para->getParD()->neighborZ,
+			para->getParD()->numberOfNodes,
+			para->getParD()->numberofthreads,
+			para->getParD()->distributions.f[0],
+			para->getParD()->distributionsAD.f[0],
+            para->getParD()->forcing,
+			para->getParD()->isEvenTimestep);
+    } else {
+		CalcMacCompSP27(
+			para->getParD()->velocityX,
+			para->getParD()->velocityY,
+			para->getParD()->velocityZ,
+			para->getParD()->rho,
+			para->getParD()->pressure,
+			para->getParD()->typeOfGridNode,
+			para->getParD()->neighborX,
+			para->getParD()->neighborY,
+			para->getParD()->neighborZ,
+			para->getParD()->numberOfNodes,
+			para->getParD()->numberofthreads,
+			para->getParD()->distributions.f[0],
+			para->getParD()->isEvenTimestep);
+	}
 }
 
 
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h b/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h
index bbe25173a87369a1dbf082f77005aa8076e7ebd3..ac4d43ac785e7ae8134d367a3c30b806e8e556c5 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h
+++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h
@@ -94,6 +94,25 @@ extern "C" void CalcMacCompSP27(
 	real* distributions,
 	bool isEvenTimestep);
 
+//////////////////////////////////////////////////////////////////////////
+//! \brief calculates the macroscopic values from distribution functions
+extern "C" void CalcMacADCompSP27(
+	real* velocityX,
+	real* velocityY,
+	real* velocityZ,
+	real* rho,
+	real* pressure,
+	uint* typeOfGridNode,
+	uint* neighborX,
+	uint* neighborY,
+	uint* neighborZ,
+	uint size_Mat,
+	uint numberOfThreads, 
+	real* distributions,
+    real* distributionsAD,
+    real* forces,
+	bool isEvenTimestep);
+
 //////////////////////////////////////////////////////////////////////////
 //! \brief defines the behavior of a velocity boundary condition (plain bounce-back plus velocity)
 extern "C" void QVelDevicePlainBB27(
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
index fde257f57fb30615829042f2eccfe0ca6dcc62cf..32cca83ed57dcad805202e0f9e5c0649b5c6e531 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
+++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
@@ -94,6 +94,24 @@ extern "C" __global__ void LBCalcMacCompSP27(
 	real* distributions,
 	bool isEvenTimestep);
 
+//////////////////////////////////////////////////////////////////////////
+//! \brief \ref LBCalcMacCompSP27 : device function to calculate macroscopic values from distributions
+extern "C" __global__ void LBCalcMacADCompSP27(
+	real* velocityX,
+	real* velocityY,
+	real* velocityZ,
+	real* rho,
+	real* pressure,
+	uint* typeOfGridNode,
+	uint* neighborX,
+	uint* neighborY,
+	uint* neighborZ,
+	uint size_Mat,
+	real* distributions,
+    real* distributionsAD,
+    real* forces,
+	bool isEvenTimestep);
+
 //////////////////////////////////////////////////////////////////////////
 //! \brief \ref QVelDevPlainBB27 : device function for the velocity boundary condition (plain bounce-back + velocity)
 extern "C" __global__ void QVelDevPlainBB27(
diff --git a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
index 3791556fb7e9d957f0771ab746a0d002e7a92c94..784793b54018f322fc48cef5660a618fad1ecd99 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -160,7 +160,46 @@ extern "C" void CalcMacCompSP27(
 		size_Mat,
 		distributions,
 		isEvenTimestep);
-	getLastCudaError("LBCalcMacSP27 execution failed");
+	getLastCudaError("LBCalcMacCompSP27 execution failed");
+}
+//////////////////////////////////////////////////////////////////////////
+extern "C" void CalcMacADCompSP27(
+	real* velocityX,
+	real* velocityY,
+	real* velocityZ,
+	real* rho,
+	real* pressure,
+	uint* typeOfGridNode,
+	uint* neighborX,
+	uint* neighborY,
+	uint* neighborZ,
+	uint size_Mat,
+	uint numberOfThreads,
+	real* distributions,
+    real* distributionsAD,
+    real* forces,
+	bool isEvenTimestep)
+{
+	int Grid = (size_Mat / numberOfThreads) + 1;
+	dim3 grid(Grid, 1, 1);
+	dim3 threads(numberOfThreads, 1, 1);
+
+	LBCalcMacADCompSP27<<<grid, threads>>>(
+		velocityX,
+		velocityY,
+		velocityZ,
+		rho,
+		pressure,
+		typeOfGridNode,
+		neighborX,
+		neighborY,
+		neighborZ,
+		size_Mat,
+		distributions,
+        distributionsAD,
+        forces,
+		isEvenTimestep);
+	getLastCudaError("LBCalcMacADCompSP27 execution failed");
 }
 //////////////////////////////////////////////////////////////////////////
 extern "C" void QVelDevicePlainBB27(