Initialize at openlb-1-3

1 files changed, 1394 insertions, 0 deletions

diff --git a/src/functors/lattice/turbulentF3D.hh b/src/functors/lattice/turbulentF3D.hh
new file mode 100644
index 0000000..66983ab
--- /dev/null
+++ b/src/functors/lattice/turbulentF3D.hh
@@ -0,0 +1,1394 @@
+/*  This file is part of the OpenLB library
+ *
+ *  Copyright (C) 2013 Patrick Nathen, Mathias J. Krause
+ *  E-mail contact: info@openlb.net
+ *  The most recent release of OpenLB can be downloaded at
+ *  <http://www.openlb.net/>
+ *
+ *  This program 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 2
+ *  of the License, or (at your option) any later version.
+ *
+ *  This program 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 this program; if not, write to the Free
+ *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ *  Boston, MA  02110-1301, USA.
+*/
+
+#ifndef TURBULENT_F_3D_HH
+#define TURBULENT_F_3D_HH
+
+#include<vector>
+#include<cmath>
+
+#include "turbulentF3D.h"
+#include "blockLatticeLocalF3D.h"
+#include "dynamics/smagorinskyBGKdynamics.h"
+#include "core/superLattice3D.h"
+#include "core/finiteDifference.h"
+#include "geometry/superGeometry3D.h"
+#include "utilities/vectorHelpers.h"
+#include "dynamics/lbHelpers.h"  // for computation of lattice rho and velocity
+
+
+namespace olb {
+
+
+///////////////////////////// SuperLatticeYplus3D //////////////////////////////
+template <typename T, typename DESCRIPTOR>
+SuperLatticeYplus3D<T,DESCRIPTOR>::SuperLatticeYplus3D(SuperLattice3D<T,DESCRIPTOR>& sLattice,
+    const UnitConverter<T,DESCRIPTOR>& converter, SuperGeometry3D<T>& superGeometry,
+    IndicatorF3D<T>& indicator, const int material )
+  : SuperLatticePhysF3D<T,DESCRIPTOR>(sLattice,converter,1),
+    _superGeometry(superGeometry), _indicator(indicator), _material(material)
+{
+  this->getName() = "yPlus";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool SuperLatticeYplus3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  int globIC = input[0];
+  int locix = input[1];
+  int lociy = input[2];
+  int lociz = input[3];
+
+  output[0]=T();
+  if ( this->_sLattice.getLoadBalancer().rank(globIC) == singleton::mpi().getRank() ) {
+    std::vector<T> normalTemp(3,T());
+    std::vector<T> normal(3,T());       // normalized
+    T counter = T();
+    T distance = T();
+    if (_superGeometry.get(input) == 1) {
+      for (int iPop = 1; iPop < DESCRIPTOR::q; ++iPop) {
+        if (_superGeometry.get(input[0],
+                               input[1] + descriptors::c<DESCRIPTOR>(iPop,0),
+                               input[2] + descriptors::c<DESCRIPTOR>(iPop,1),
+                               input[3] + descriptors::c<DESCRIPTOR>(iPop,2)) == _material) {
+          counter++;
+          normalTemp[0] += descriptors::c<DESCRIPTOR>(iPop,0);
+          normalTemp[1] += descriptors::c<DESCRIPTOR>(iPop,1);
+          normalTemp[2] += descriptors::c<DESCRIPTOR>(iPop,2);
+        }
+      }
+      if ( !util::nearZero(counter) ) {
+        // get physical Coordinates at intersection
+
+        std::vector<T> physR(3, T());
+        _superGeometry.getCuboidGeometry().getPhysR(&(physR[0]), &(input[0]));
+
+        T voxelSize = _superGeometry.getCuboidGeometry().get(globIC).getDeltaR();
+
+        normal = util::normalize(normalTemp);
+
+        std::vector<T> direction(normal);
+        direction[0] = voxelSize*normal[0]*2.;
+        direction[1] = voxelSize*normal[1]*2.;
+        direction[2] = voxelSize*normal[2]*2.;
+
+        // calculate distance to STL file
+        if ( _indicator.distance(distance, physR, direction) ) {
+          // call stress at this point
+          T rho;
+          T u[3];
+          T pi[6];
+          this->_sLattice.getBlockLattice(this->_sLattice.getLoadBalancer().loc(globIC)).get(locix, lociy, lociz).computeRhoU(rho, u);
+          this->_sLattice.getBlockLattice(this->_sLattice.getLoadBalancer().loc(globIC)).computeStress(locix, lociy, lociz, pi);
+
+          // Compute phys stress tau = mu*du/dx
+          T omega = 1. / this->_converter.getLatticeRelaxationTime();
+          T dt = this->_converter.getConversionFactorTime();
+          T physFactor = -omega*descriptors::invCs2<T,DESCRIPTOR>()/rho/2./dt*this->_converter.getPhysDensity(rho)*this->_converter.getPhysViscosity();
+
+          //  Totel Stress projected from cell in normal direction on obstacle
+          T Rx = pi[0]*physFactor*normal[0] + pi[1]*physFactor*normal[1] + pi[2]*physFactor*normal[2];
+          T Ry = pi[1]*physFactor*normal[0] + pi[3]*physFactor*normal[1] + pi[4]*physFactor*normal[2];
+          T Rz = pi[2]*physFactor*normal[0] + pi[4]*physFactor*normal[1] + pi[5]*physFactor*normal[2];
+
+          // Stress appearing as pressure in corresponding direction is calculated and substracted
+          T R_res_pressure = normal[0]*pi[0]*physFactor*normal[0] + normal[0]*pi[1]*physFactor*normal[1] + normal[0]*pi[2]*physFactor*normal[2]
+                             +normal[1]*pi[1]*physFactor*normal[0] + normal[1]*pi[3]*physFactor*normal[1] + normal[1]*pi[4]*physFactor*normal[2]
+                             +normal[2]*pi[2]*physFactor*normal[0] + normal[2]*pi[4]*physFactor*normal[1] + normal[2]*pi[5]*physFactor*normal[2];
+
+          Rx -= R_res_pressure *normal[0];
+          Ry -= R_res_pressure *normal[1];
+          Rz -= R_res_pressure *normal[2];
+
+          T tau_wall = sqrt(Rx*Rx+Ry*Ry+Rz*Rz);
+          T u_tau = sqrt(tau_wall/this->_converter.getPhysDensity(rho));
+          //y_plus
+          output[0] = distance*u_tau / this->_converter.getPhysViscosity();
+        } // if 4
+      }
+    }
+  }
+  return true;
+}
+
+/*template <typename T, typename DESCRIPTOR>
+BlockLatticeADM3D<T,DESCRIPTOR>::BlockLatticeADM3D
+(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, T sigma, int order, bool adaptive, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice,5), _sigma(sigma), _order(order), _adaptive(adaptive), _converter(converter)
+{
+  this->getName() = "ADMfilter";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeADM3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  // Declaration of all variables needed for filtering
+
+  int globX = input[0];
+  int globY = input[1];
+  int globZ = input[2];
+
+  T output000[4] = {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY, globZ).computeRhoU(output000[0],output000+1);
+
+  T outputp00[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX+1, globY, globZ).computeRhoU(outputp00[0],outputp00+1);
+
+  T output2p00[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX+2, globY, globZ).computeRhoU(output2p00[0],output2p00+1);
+
+  T outputn00[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX-1, globY, globZ).computeRhoU(outputn00[0],outputn00+1);
+
+  T output2n00[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX-2, globY, globZ).computeRhoU(output2n00[0],output2n00+1);
+
+
+  T output0p0[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY+1, globZ).computeRhoU(output0p0[0],output0p0+1);
+
+  T output02p0[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY+2, globZ).computeRhoU(output02p0[0],output02p0+1);
+
+  T output0n0[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY-1, globZ).computeRhoU(output0n0[0],output0n0+1);
+
+  T output02n0[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY-2, globZ).computeRhoU(output02n0[0],output02n0+1);
+
+
+  T output00p[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY, globZ+1).computeRhoU(output00p[0],output00p+1);
+
+  T output002p[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY, globZ+2).computeRhoU(output002p[0],output002p+1);
+
+  T output00n[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY, globZ-1).computeRhoU(output00n[0],output00n+1);
+
+  T output002n[4]= {1.,0.,0.,0.};
+  this-> _blockLattice.get(globX, globY, globZ-2).computeRhoU(output002n[0],output002n+1);
+
+ T relax=_sigma;
+
+   if (_adaptive==true ){
+    ///////////////////////////////////////////////DISS VERSION///////////////////////////////////////////////////
+
+  //  T diss = dissipation(u_000)[0];
+
+    // std::cout << "diss: "<< diss << std::endl;
+
+  //  T* avDiss = this-> _blockLattice.get(globX, globY , globZ)[localAvDissBeginsAt];
+
+    // // std::cout <<"avDiss:" << *avDiss << std::endl;
+
+  //   *avDiss = (*avDiss * this->_blockLattice.getStatistics().getTime() + diss) / (this->_blockLattice.getStatistics().getTime() + (int) 1 );
+
+    // // std::cout <<"avDiss nach mittelung:" << *avDiss << std::endl;
+
+   //  T TKE = 0.5*(velocity(u_000)[0]*velocity(u_000)[0]+velocity(u_000)[1]*velocity(u_000)[1]+velocity(u_000)[2]+velocity(u_000)[2]);
+
+   //  T* avTKE = this-> _blockLattice.get(globX, globY , globZ)[localAvTKEBeginsAt];
+
+   //  *avTKE = (*avTKE * this->_blockLattice.getStatistics().getTime() + TKE) / (this->_blockLattice.getStatistics().getTime() + (int) 1 );
+
+    // std::cout << "TKE: "<< *avTKE << std::endl;
+
+
+   //  relax = sqrt((diss - *avDiss)*(diss - *avDiss));// / (*avTKE * converter.getDeltaT());
+    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+    ////////////////////////////////////////////// Stress Version ////////////////////////////////////////////////
+    T stress[9] = {0.,0.,0.,0.,0.,0.,0.,0.,0.};
+    this-> _blockLattice.computeStress(globX, globY, globZ, stress);
+
+    T ux = stress(u_000)[0];
+    T uy = stress(u_000)[1];
+    T uz = stress(u_000)[2];
+
+    T vx = stress(u_000)[3];
+    T vy = stress(u_000)[4];
+    T vz = stress(u_000)[5];
+
+    T wx = stress(u_000)[6];
+    T wy = stress(u_000)[7];
+    T wz = stress(u_000)[8];
+
+
+    T ux = stress[0];
+    T uy = stress[1];
+    T uz = stress[2];
+
+    T vx = stress[3];
+    T vy = stress[4];
+    T vz = stress[5];
+
+    T wx = stress[6];
+    T wy = stress[7];
+    T wz = stress[8];
+
+    T norm = sqrt(  (wx*uz + vx*uy + ux*ux) + (wy*vz + vy*vy + uy*vx) + (wz*wz + vz*wy + uz*wx) ) ;
+
+
+    T* avNorm = this-> _blockLattice.get(globX, globY , globZ)[_localAvDissBeginsAt];
+
+    *avNorm = (*avNorm * this->_blockLattice.getStatistics().getTime() + norm) / (this->_blockLattice.getStatistics().getTime() + (int) 1 );
+
+    // relax = sigma;
+    // / (*avTKE * converter.getDeltaT());
+
+
+    // if (this->_blockLattice.getStatistics().getTime() >= 30000){
+
+
+     relax = sqrt((norm - *avNorm)*(norm - *avNorm)) ;
+
+      // std::cout << "adaptive relaxation: " << relax <<  " time: "<<this->_blockLattice.getStatistics().getTime()<< endl;
+    // }
+
+   }
+
+  if (_order==2) { // second order
+    T d_0 = 6./16.;
+    T d_1 = -4./16.;
+    T d_2 = 1./16.;
+
+    output[0] = output000[0] - relax*(d_2*(output2n00[0]+output02n0[0]+output002n[0]) +
+                                       d_1*(outputn00[0]+output0n0[0]+output00n[0])+
+                                       d_0*(output000[0]+output000[0]+output000[0])+
+                                       d_1*(outputp00[0]+output0p0[0]+output00p[0])+
+                                       d_2*(output2p00[0]+output02p0[0]+output002p[0]) );
+    for (int i = 1; i < 4; ++i ) {
+      output[i] = output000[i] - relax*(d_2*(output2n00[i] + output02n0[i] + output002n[i]) +
+                                         d_1*(outputn00[i] + output0n0[i] + output00n[i])+
+                                         d_0*(output000[i] + output000[i] + output000[i])+
+                                         d_1*(outputp00[i] + output0p0[i] + output00p[i])+
+                                         d_2*(output2p00[i] + output02p0[i] + output002p[i]) );
+    }
+  } else { // third order
+
+    T output3p00[4]= {1.,0.,0.,0.};
+    this-> _blockLattice.get(globX+3, globY, globZ).computeRhoU(output3p00[0],output3p00+1);
+
+    T output3n00[4]= {1.,0.,0.,0.};
+    this-> _blockLattice.get(globX-3, globY, globZ).computeRhoU(output3n00[0],output3n00+1);
+
+    T output03p0[4]= {1.,0.,0.,0.};
+    this-> _blockLattice.get(globX, globY+3, globZ).computeRhoU(output03p0[0],output03p0+1);
+
+    T output03n0[4]= {1.,0.,0.,0.};
+    this-> _blockLattice.get(globX, globY-3, globZ).computeRhoU(output03n0[0],output03n0+1);
+
+    T output003p[4]= {1.,0.,0.,0.};
+    this-> _blockLattice.get(globX, globY, globZ+3).computeRhoU(output003p[0],output003p+1);
+
+    T output003n[4] = {1.,0.,0.,0.};
+    this-> _blockLattice.get(globX, globY, globZ-3).computeRhoU(output003n[0],output003n+1);
+
+    T d_0 = 5./16.;
+    T d_1 = -15./64.;
+    T d_2 = 3./32.;
+    T d_3 = -1./64.;
+    output[0] = output000[0] - _sigma*(d_3*(output3n00[0]+output03n0[0]+output003n[0])+
+                                       d_2*(output2n00[0]+output02n0[0]+output002n[0]) +
+                                       d_1*(outputn00[0]+output0n0[0]+output00n[0])+
+                                       d_0*(output000[0]+output000[0]+output000[0])+
+                                       d_1*(outputp00[0]+output0p0[0]+output00p[0])+
+                                       d_2*(output2p00[0]+output02p0[0]+output002p[0])+
+                                       d_3*(output3p00[0]+output03p0[0]+output003p[0]) );
+    for (int i = 1; i < 4; ++i ) {
+      output[i] = output000[i] - _sigma*(d_3*(output3n00[i]+output03n0[i]+output003n[i])+
+                                         d_2*(output2n00[i]+output02n0[i]+output002n[i]) +
+                                         d_1*(outputn00[i]+output0n0[i]+output00n[i])+
+                                         d_0*(output000[i]+output000[i]+output000[i])+
+                                         d_1*(outputp00[i]+output0p0[i]+output00p[i])+
+                                         d_2*(output2p00[i]+output02p0[i]+output002p[i])+
+                                         d_3*(output3p00[i]+output03p0[i]+output003p[i]) );
+    }
+  }
+  output[4]=relax;
+  return true;
+}
+
+template <typename T, typename DESCRIPTOR>
+void BlockLatticeADM3D<T,DESCRIPTOR>::execute(const int input[])
+{
+  T output[5] = {T(),T(),T(),T(),T()};
+  this->operator()(output,input);
+  this-> _blockLattice.get(input[0],input[1],input[2]).defineField<descriptors::FIL_RHO>( &output[0] );
+  this-> _blockLattice.get(input[0],input[1],input[2]).defineField<descriptors::LOCAL_FIL_VEL_X>( &output[1]);
+  this-> _blockLattice.get(input[0],input[1],input[2]).defineField<descriptors::LOCAL_FIL_VEL_Y>( &output[2]);
+  this-> _blockLattice.get(input[0],input[1],input[2]).defineField<descriptors::LOCAL_FIL_VEL_Z>( &output[3]);
+  this-> _blockLattice.get(input[0],input[1],input[2]).defineField<descriptors::LOCAL_SIGMA_ADM>( &output[4]);
+}
+
+template <typename T, typename DESCRIPTOR>
+void BlockLatticeADM3D<T,DESCRIPTOR>::execute()
+{
+  int nX = this-> _blockLattice.getNx();
+  int nY = this-> _blockLattice.getNy();
+  int nZ = this-> _blockLattice.getNz();
+  int i[3];
+  for (i[0]=0; i[0]<nX; ++i[0]) {
+    for (i[1]=0; i[1]<nY; ++i[1]) {
+      for (i[2]=0; i[2]<nZ; ++i[2]) {
+        execute(i);
+      }
+    }
+  }
+}
+
+
+template <typename T, typename DESCRIPTOR>
+SuperLatticeADM3D<T,DESCRIPTOR>::SuperLatticeADM3D
+(SuperLattice3D<T,DESCRIPTOR>& sLattice, T sigma, int order, bool adaptive, const UnitConverter<T,DESCRIPTOR>& converter)
+  : SuperLatticeF3D<T,DESCRIPTOR>(sLattice,5), _sigma(sigma), _order(order), _adaptive(adaptive), _converter(converter)
+{
+  this->getName() = "ADMfilter";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool SuperLatticeADM3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  int globIC = input[0];
+  if ( this->_sLattice.getLoadBalancer().rank(globIC) == singleton::mpi().getRank() ) {
+    int inputLocal[3]= {};
+    T overlap = this->_sLattice.getOverlap();
+    inputLocal[0] = input[1] + overlap;
+    inputLocal[1] = input[2] + overlap;
+    inputLocal[2] = input[3] + overlap;
+
+    BlockLatticeADM3D<T,DESCRIPTOR> blockLatticeF( this->_sLattice.getExtendedBlockLattice(this->_sLattice.getLoadBalancer().loc(globIC)), _sigma, _order,_adaptive, _converter );
+
+    blockLatticeF(output,inputLocal);
+    return true;
+  } else {
+    return false;
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+void SuperLatticeADM3D<T,DESCRIPTOR>::execute(SuperGeometry3D<T>& superGeometry, const int material)
+{
+  this->_sLattice.communicate();
+  CuboidGeometry3D<T>& cGeometry =  this->_sLattice.getCuboidGeometry();
+  LoadBalancer<T>& load = this->_sLattice.getLoadBalancer();
+  int overlap = this->_sLattice.getOverlap();
+
+  for (int iC = 0; iC < load.size(); ++iC) {
+    BlockLatticeADM3D<T,DESCRIPTOR> blockLatticeF( this->_sLattice.getExtendedBlockLattice(iC), _sigma, _order, _adaptive, _converter );
+    int nX = cGeometry.get(load.glob(iC)).getNx();
+    int nY = cGeometry.get(load.glob(iC)).getNy();
+    int nZ = cGeometry.get(load.glob(iC)).getNz();
+
+    int i[3];
+    for (i[0]=overlap; i[0]<nX+overlap; ++i[0]) {
+      for (i[1]=overlap; i[1]<nY+overlap; ++i[1]) {
+        for (i[2]=overlap; i[2]<nZ+overlap; ++i[2]) {
+          if (superGeometry.getExtendedBlockGeometry(iC).get(i[0],i[1],i[2]) == material) {
+            blockLatticeF.execute(i);
+          }
+        }
+      }
+    }
+  }
+}
+
+*/
+
+////////////////////////BlockFiniteDifference3D//////////////////////////////////
+template <typename T>
+BlockFiniteDifference3D<T>::BlockFiniteDifference3D
+(BlockGeometryStructure3D<T>& blockGeometry, BlockF3D<T>& blockFunctor, std::list<int>& matNumber)
+  : BlockF3D<T>(blockFunctor.getBlockStructure(), 3*blockFunctor.getTargetDim()), _blockGeometry(blockGeometry), _blockFunctor(blockFunctor), _matNumber(matNumber)
+{
+  this->getName() = "FiniteDifference";
+  _targetDim = _blockFunctor.getTargetDim();
+  _n[0] = this-> _blockGeometry.getNx()-1;
+  _n[1] = this-> _blockGeometry.getNy()-1;
+  _n[2] = this-> _blockGeometry.getNz()-1;
+
+}
+
+template <typename T>
+bool BlockFiniteDifference3D<T>::operator() (T output[], const int input[])
+{
+//  // derivation tensor
+  std::vector<std::vector<T>> fdGrad;
+
+  fdGrad.resize(_targetDim);
+  for (int i = 0; i < _targetDim; i++) {
+    fdGrad[i].resize(3);
+  }
+
+  for (int i = 0; i < 3; i++) {
+    int fInput_p[3];
+    fInput_p[0] = input[0];
+    fInput_p[1] = input[1];
+    fInput_p[2] = input[2];
+    fInput_p[i]+=1;
+
+    int fInput_2p[3];
+    fInput_2p[0] = input[0];
+    fInput_2p[1] = input[1];
+    fInput_2p[2] = input[2];
+    fInput_2p[i]+=2;
+
+    int fInput_3p[3];
+    fInput_3p[0] = input[0];
+    fInput_3p[1] = input[1];
+    fInput_3p[2] = input[2];
+    fInput_3p[i]+=3;
+
+    int fInput_4p[3];
+    fInput_4p[0] = input[0];
+    fInput_4p[1] = input[1];
+    fInput_4p[2] = input[2];
+    fInput_4p[i]+=4;
+
+    int fInput_n[3];
+    fInput_n[0] = input[0];
+    fInput_n[1] = input[1];
+    fInput_n[2] = input[2];
+    fInput_n[i]-=1;
+
+    int fInput_2n[3];
+    fInput_2n[0] = input[0];
+    fInput_2n[1] = input[1];
+    fInput_2n[2] = input[2];
+    fInput_2n[i]-=2;
+
+    int fInput_3n[3];
+    fInput_3n[0] = input[0];
+    fInput_3n[1] = input[1];
+    fInput_3n[2] = input[2];
+    fInput_3n[i]-=3;
+
+    int fInput_4n[3];
+    fInput_4n[0] = input[0];
+    fInput_4n[1] = input[1];
+    fInput_4n[2] = input[2];
+    fInput_4n[i]-=4;
+
+    T fOutput[_targetDim];
+    _blockFunctor(fOutput,input);
+
+    if (input[i] < 3) {
+      if (std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_2p[0], fInput_2p[1], fInput_2p[2])) == _matNumber.end()) {
+        T fOutput_p[_targetDim];
+        _blockFunctor(fOutput_p,fInput_p);
+        for (int j=0; j < _targetDim; j++) {
+          fdGrad[j][i]= -fOutput[j] + fOutput_p[j];
+        }
+      } else {
+        T fOutput_p[_targetDim];
+        _blockFunctor(fOutput_p,fInput_p);
+        T fOutput_2p[_targetDim];
+        _blockFunctor(fOutput_2p,fInput_2p);
+        for (int j=0; j < _targetDim; j++) {
+          fdGrad[j][i]=fd::boundaryGradient(fOutput[j], fOutput_p[j], fOutput_2p[j]);
+        }
+      }
+    } else if (input[i] > _n[i]-3) {
+      if (std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_2n[0], fInput_2n[1], fInput_2n[2])) == _matNumber.end()) {
+        T fOutput_n[_targetDim];
+        _blockFunctor(fOutput_n,fInput_n);
+        for (int j=0; j < _targetDim; j++) {
+          fdGrad[j][i]= -fOutput_n[j] + fOutput[j];
+        }
+      } else {
+        T fOutput_n[_targetDim];
+        _blockFunctor(fOutput_n,fInput_n);
+        T fOutput_2n[_targetDim];
+        _blockFunctor(fOutput_2n,fInput_2n);
+        for (int j=0; j < _targetDim; j++) {
+          fdGrad[j][i]=fd::boundaryGradient(-fOutput[j], -fOutput_n[j], -fOutput_2n[j]);
+        }
+      }
+    } else {
+      if ( std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_n[0], fInput_n[1], fInput_n[2])) == _matNumber.end()  &&
+           std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_p[0], fInput_p[1], fInput_p[2])) == _matNumber.end() ) {
+        for (int j=0; j < _targetDim; j++) {
+          fdGrad[j][i]=0.;
+        }
+        // boundary treatment with Second-order asymmetric gradient
+      } else if (std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_n[0], fInput_n[1], fInput_n[2])) == _matNumber.end()) {
+        if (std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_2p[0], fInput_2p[1], fInput_2p[2])) == _matNumber.end()) {
+          T fOutput_p[_targetDim];
+          _blockFunctor(fOutput_p,fInput_p);
+          for (int j=0; j < _targetDim; j++) {
+            fdGrad[j][i]= -fOutput[j] + fOutput_p[j];
+          }
+        } else {
+          T fOutput_p[_targetDim];
+          _blockFunctor(fOutput_p,fInput_p);
+          T fOutput_2p[_targetDim];
+          _blockFunctor(fOutput_2p,fInput_2p);
+          for (int j=0; j < _targetDim; j++) {
+            fdGrad[j][i]=fd::boundaryGradient(fOutput[j], fOutput_p[j], fOutput_2p[j]);
+          }
+        }
+      } else if (std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_p[0], fInput_p[1], fInput_p[2])) == _matNumber.end() ) {
+        if (std::find(_matNumber.begin(), _matNumber.end(), _blockGeometry.get(fInput_2n[0], fInput_2n[1], fInput_2n[2])) == _matNumber.end()) {
+          T fOutput_n[_targetDim];
+          _blockFunctor(fOutput_n,fInput_n);
+          for (int j=0; j < _targetDim; j++) {
+            fdGrad[j][i]= -fOutput_n[j] + fOutput[j];
+          }
+        } else {
+          T fOutput_n[_targetDim];
+          _blockFunctor(fOutput_n,fInput_n);
+          T fOutput_2n[_targetDim];
+          _blockFunctor(fOutput_2n,fInput_2n);
+          for (int j=0; j < _targetDim; j++) {
+            fdGrad[j][i]=fd::boundaryGradient(-fOutput[j], -fOutput_n[j], -fOutput_2n[j]);
+          }
+        }
+      } else {
+        //inner domain 8th order central difference
+        T fOutput_n[_targetDim];
+        _blockFunctor(fOutput_n,fInput_n);
+
+        T fOutput_2n[_targetDim];
+        _blockFunctor(fOutput_2n,fInput_2n);
+
+        T fOutput_3n[_targetDim];
+        _blockFunctor(fOutput_3n,fInput_3n);
+
+        T fOutput_4n[_targetDim];
+        _blockFunctor(fOutput_4n,fInput_4n);
+
+        T fOutput_p[_targetDim];
+        _blockFunctor(fOutput_p,fInput_p);
+
+        T fOutput_2p[_targetDim];
+        _blockFunctor(fOutput_2p,fInput_2p);
+
+        T fOutput_3p[_targetDim];
+        _blockFunctor(fOutput_3p,fInput_3p);
+
+        T fOutput_4p[_targetDim];
+        _blockFunctor(fOutput_4p,fInput_4p);
+        for (int j=0; j < _targetDim; j++) {
+          //fdGrad[j][i]=fd::centralGradient(fOutput_p[j], fOutput_n[j]);
+          fdGrad[j][i]=((T)672*(fOutput_p[j]-fOutput_n[j])+(T)168*(fOutput_2n[j]-fOutput_2p[j])
+                        +(T)32*(fOutput_3p[j]-fOutput_3n[j])+(T)3*(fOutput_4n[j]-fOutput_4p[j])) / 840.;
+        }
+      }
+    }
+    for (int i=0; i < 3; i++) {
+      for (int j=0; j < _targetDim; j++) {
+        output[i*3+j] = fdGrad[i][j];
+      }
+    }
+  }
+  return true;
+}
+
+
+
+////////////////////////SuperFiniteDifference3D//////////////////////////////////
+template <typename T>
+SuperFiniteDifference3D<T>::SuperFiniteDifference3D
+(SuperGeometry3D<T>& sGeometry, SuperF3D<T>& sFunctor, std::list<int>& matNumber) : SuperF3D<T>(sFunctor.getSuperStructure(),3*sFunctor.getTargetDim()),
+  _sGeometry(sGeometry) ,_sFunctor(sFunctor), _matNumber(matNumber)
+{
+  this->getName() = "FiniteDifference";
+  int maxC = this->_superStructure.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int iC = 0; iC < maxC; iC++) {
+    this->_blockF.emplace_back(new BlockFiniteDifference3D<T> ( _sGeometry.getBlockGeometry(iC), _sFunctor.getBlockF(iC), _matNumber ));
+  }
+}
+
+template <typename T>
+bool SuperFiniteDifference3D<T>::operator() (T output[], const int input[])
+{
+  if (this->_superStructure.getLoadBalancer().rank(input[0]) == singleton::mpi().getRank()) {
+    return this->getBlockF(this->_superStructure.getLoadBalancer().loc(input[0]) )(output,&input[1]);
+  } else {
+    return false;
+  }
+}
+
+////////////////////////BlockPhysFiniteDifference3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+BlockPhysFiniteDifference3D<T,DESCRIPTOR>::BlockPhysFiniteDifference3D
+(BlockF3D<T>& blockFinDiff, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockF3D<T>(blockFinDiff.getBlockStructure(), 3*blockFinDiff.getTargetDim()), _blockFinDiff(blockFinDiff), _converter(converter)
+{
+  this->getName() = "PhysFiniteDifference";
+  _targetDim = _blockFinDiff.getTargetDim();
+
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockPhysFiniteDifference3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  _blockFinDiff(output,input);
+  for (int i = 0; i < _targetDim; i++) {
+    output[i] /= _converter.getConversionFactorLength();
+  }
+  return true;
+}
+
+
+
+////////////////////////SuperPhysFiniteDifference3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+SuperPhysFiniteDifference3D<T,DESCRIPTOR>::SuperPhysFiniteDifference3D
+(SuperGeometry3D<T>& sGeometry, SuperF3D<T>& sFunctor, std::list<int>& matNumber, const UnitConverter<T,DESCRIPTOR>& converter) : SuperF3D<T>(sFunctor.getSuperStructure(),3*sFunctor.getTargetDim()),
+  _sFinDiff(sGeometry,sFunctor,matNumber),_converter(converter)
+{
+  this->getName() = "PhysFiniteDifference";
+  int maxC = this->_superStructure.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int iC = 0; iC < maxC; iC++) {
+    this->_blockF.emplace_back(new BlockPhysFiniteDifference3D<T,DESCRIPTOR> (_sFinDiff.getBlockF(iC), _converter ));
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+bool SuperPhysFiniteDifference3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  if (this->_superStructure.getLoadBalancer().rank(input[0]) == singleton::mpi().getRank()) {
+    return this->getBlockF(this->_superStructure.getLoadBalancer().loc(input[0]) )(output,&input[1]);
+  } else {
+    return false;
+  }
+}
+////////////////////////BlockLatticeVelocityGradientFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+BlockLatticeVelocityGradientFD3D<T,DESCRIPTOR>::BlockLatticeVelocityGradientFD3D
+(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, BlockF3D<T>& blockFinDiff)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 9), _blockFinDiff(blockFinDiff)
+{
+  this->getName() = "VelocityGradientFD";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeVelocityGradientFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  //1 dudx 2 dudy 3 dudz
+  //4 dydx 5 dydy 6 dydz
+  //7 dwdx 8 dwdy 9 dwdz
+  _blockFinDiff(output,input);
+  return true;
+}
+
+////////////////////////BlockLatticeExternalVelocityGradientFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+BlockLatticeExternalVelocityGradientFD3D<T,DESCRIPTOR>::BlockLatticeExternalVelocityGradientFD3D
+(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, BlockF3D<T>& blockFinDiff)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 9), _blockFinDiff(blockFinDiff)
+{
+  this->getName() = "externalVelocityGradientFD";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeExternalVelocityGradientFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  //1 dudx 2 dudy 3 dudz
+  //4 dydx 5 dydy 6 dydz
+  //7 dwdx 8 dwdy 9 dwdz
+  _blockFinDiff(output,input);
+  return true;
+}
+
+////////////////////////SuperLatticeVelocityGradientFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+SuperLatticeVelocityGradientFD3D<T,DESCRIPTOR>::SuperLatticeVelocityGradientFD3D
+(SuperGeometry3D<T>& sGeometry, SuperLattice3D<T,DESCRIPTOR>& sLattice, std::list<int>& matNumber) : SuperLatticeF3D<T,DESCRIPTOR>(sLattice,9),
+  _sVelocity(sLattice), _sFinDiff(sGeometry, _sVelocity, matNumber)
+{
+  this->getName() = "VelocityGradientFD";
+  int maxC = this->_superStructure.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int iC = 0; iC < maxC; iC++) {
+    this->_blockF.emplace_back(new BlockLatticeVelocityGradientFD3D<T,DESCRIPTOR> (this->_sLattice.getBlockLattice(iC), this->_sFinDiff.getBlockF(iC)));
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+bool SuperLatticeVelocityGradientFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  if (this->_sLattice.getLoadBalancer().rank(input[0]) == singleton::mpi().getRank()) {
+    return this->getBlockF(this->_sLattice.getLoadBalancer().loc(input[0]) )(output,&input[1]);
+  } else {
+    return false;
+  }
+}
+////////////////////////SuperLatticeExternalVelocityGradientFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+SuperLatticeExternalVelocityGradientFD3D<T,DESCRIPTOR>::SuperLatticeExternalVelocityGradientFD3D
+(SuperGeometry3D<T>& sGeometry, SuperLattice3D<T,DESCRIPTOR>& sLattice, std::list<int>& matNumber) : SuperLatticeF3D<T,DESCRIPTOR>(sLattice,9),
+  _sVelocity(sLattice), _sFinDiff(sGeometry, _sVelocity, matNumber)
+{
+  this->getName() = "externalVelocityGradientFD";
+  int maxC = this->_superStructure.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int iC = 0; iC < maxC; iC++) {
+    this->_blockF.emplace_back(new BlockLatticeExternalVelocityGradientFD3D<T,DESCRIPTOR> (this->_sLattice.getBlockLattice(iC), this->_sFinDiff.getBlockF(iC)));
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+bool SuperLatticeExternalVelocityGradientFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  if (this->_sLattice.getLoadBalancer().rank(input[0]) == singleton::mpi().getRank()) {
+    return this->getBlockF(this->_sLattice.getLoadBalancer().loc(input[0]) )(output,&input[1]);
+  } else {
+    return false;
+  }
+}
+////////////////////////BlockLatticePhysVelocityGradientFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysVelocityGradientFD3D<T,DESCRIPTOR>::BlockLatticePhysVelocityGradientFD3D
+(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, BlockF3D<T>& blockFinDiff, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 9), _blockFinDiff(blockFinDiff), _converter(converter)
+{
+  this->getName() = "PhysVelocityGradientFD";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysVelocityGradientFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  _blockFinDiff(output,input);
+  return true;
+}
+
+
+
+////////////////////////SuperLatticePhysVelocityGradientFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+SuperLatticePhysVelocityGradientFD3D<T,DESCRIPTOR>::SuperLatticePhysVelocityGradientFD3D
+(SuperGeometry3D<T>& sGeometry, SuperLattice3D<T,DESCRIPTOR>& sLattice, std::list<int>& matNumber, const UnitConverter<T,DESCRIPTOR>& converter) : SuperLatticeF3D<T,DESCRIPTOR>(sLattice,9),
+  _sVelocity(sLattice, converter), _sFinDiff(sGeometry, _sVelocity, matNumber, converter), _converter(converter)
+{
+  this->getName() = "PhysVelocityGradientFD";
+  int maxC = this->_superStructure.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int iC = 0; iC < maxC; iC++) {
+    this->_blockF.emplace_back(new BlockLatticePhysVelocityGradientFD3D<T,DESCRIPTOR> (this->_sLattice.getBlockLattice(iC), this->_sFinDiff.getBlockF(iC), this->_converter));
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+bool SuperLatticePhysVelocityGradientFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  if (this->_sLattice.getLoadBalancer().rank(input[0]) == singleton::mpi().getRank()) {
+    return this->getBlockF(this->_sLattice.getLoadBalancer().loc(input[0]) )(output,&input[1]);
+  } else {
+    return false;
+  }
+}
+////////////////////////BlockLatticeStrainRateFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+BlockLatticeStrainRateFD3D<T,DESCRIPTOR>::BlockLatticeStrainRateFD3D
+(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, BlockF3D<T>& blockVeloGrad)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 9), _blockVeloGrad(blockVeloGrad)
+{
+  this->getName() = "StrainRateFD";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticeStrainRateFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  T velograd[9];
+  _blockVeloGrad(velograd,input);
+  output[0] = velograd[0];
+  output[1] = 0.5 * (velograd[1] + velograd[3]);
+  output[2] = 0.5 * (velograd[2] + velograd[6]);
+  output[3] = output[1];
+  output[4] = velograd[4];
+  output[5] = 0.5 * (velograd[5] + velograd[7]);
+  output[6] = output[2];
+  output[7] = output[5];
+  output[8] = velograd[8];
+  return true;
+}
+
+
+
+////////////////////////SuperLatticeStrainRateFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+SuperLatticeStrainRateFD3D<T,DESCRIPTOR>::SuperLatticeStrainRateFD3D
+(SuperGeometry3D<T>& sGeometry, SuperLattice3D<T,DESCRIPTOR>& sLattice, std::list<int>& matNumber) : SuperLatticeF3D<T,DESCRIPTOR>(sLattice,9),
+  _sVeloGrad(sGeometry, sLattice, matNumber)
+{
+  this->getName() = "StrainRateFD";
+  int maxC = this->_superStructure.getLoadBalancer().size();
+  this->_blockF.reserve(maxC);
+  for (int iC = 0; iC < maxC; iC++) {
+    this->_blockF.emplace_back(new BlockLatticeStrainRateFD3D<T,DESCRIPTOR> (this->_sLattice.getBlockLattice(iC), this->_sVeloGrad.getBlockF(iC), this->_converter));
+  }
+}
+
+template <typename T, typename DESCRIPTOR>
+bool SuperLatticeStrainRateFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  if (this->_sLattice.getLoadBalancer().rank(input[0]) == singleton::mpi().getRank()) {
+    return this->getBlockF(this->_sLattice.getLoadBalancer().loc(input[0]) )(output,&input[1]);
+  } else {
+    return false;
+  }
+}
+////////////////////////BlockLatticePhysStrainRateFD3D//////////////////////////////////
+template <typename T, typename DESCRIPTOR>
+BlockLatticePhysStrainRateFD3D<T,DESCRIPTOR>::BlockLatticePhysStrainRateFD3D
+(BlockLatticeStructure3D<T,DESCRIPTOR>& blockLattice, BlockF3D<T>& blockVeloGrad, const UnitConverter<T,DESCRIPTOR>& converter)
+  : BlockLatticeF3D<T,DESCRIPTOR>(blockLattice, 9), _blockVeloGrad(blockVeloGrad), _converter(converter)
+{
+  this->getName() = "PhysStrainRateFD";
+}
+
+template <typename T, typename DESCRIPTOR>
+bool BlockLatticePhysStrainRateFD3D<T,DESCRIPTOR>::operator() (T output[], const int input[])
+{
+  T velograd[9];
+  _blockVeloGrad(velograd,input);
+  output[0] = velograd[0];
+  output[1] = 0.5 * (velograd[1] + velograd[3]);
+  output[2] = 0.5 * (velograd[2] + velograd[6]);
+  output[3] = output[1];
+  output[4] = velograd[4];
+  output[5] = 0.5 * (velograd[5] + velograd[7]);
+  output[6] = output[2];
+  output[7] = output[5];
+  output[8] = velograd[8];
+