src_openmp4.5/Advection.hpp
// Perform the interpolation of the 4d advection
Perform the interpolation of the 4d advection
Perform the interpolation of the 4d advection
static inline float64 interp_4D(RealView4D &tmp_fn, float64 xmin[DIMENSION], 
                                  float64 inv_dx[DIMENSION], float64 xstar[DIMENSION]) {
    int ipos[DIMENSION];
    for(int j = 0; j < DIMENSION; j++) {
      xstar[j] = inv_dx[j] * (xstar[j] - xmin[j]);
    }
    #ifdef LAG_ORDER
      float64 coefx0[LAG_PTS];
      float64 coefx1[LAG_PTS];
      float64 coefx2[LAG_PTS];
      float64 coefx3[LAG_PTS];
      ipos[0] = floor(xstar[0]) - LAG_OFFSET;
      ipos[1] = floor(xstar[1]) - LAG_OFFSET;
      ipos[2] = floor(xstar[2]) - LAG_OFFSET;
      ipos[3] = floor(xstar[3]) - LAG_OFFSET;
      lag_basis((xstar[0] - ipos[0]), coefx0);
      lag_basis((xstar[1] - ipos[1]), coefx1);
      lag_basis((xstar[2] - ipos[2]), coefx2);
      lag_basis((xstar[3] - ipos[3]), coefx3);
      float64 ftmp1 = 0.;
      for(int k3 = 0; k3 <= LAG_ORDER; k3++) {
        for(int k2 = 0; k2 <= LAG_ORDER; k2++) {
          float64 ftmp2 = 0.;
          for(int k1 = 0; k1 <= LAG_ORDER; k1++) {
            float64 ftmp3 = 0.;
            for(int k0 = 0; k0 <= LAG_ORDER; k0++) {
              ftmp3 += coefx0[k0] * tmp_fn(ipos[0] + k0, ipos[1] + k1, ipos[2] + k2, ipos[3] + k3);
            }
            ftmp2 += ftmp3 * coefx1[k1];
          }
          ftmp1 += ftmp2 * coefx2[k2] * coefx3[k3];
        }
      }
      return ftmp1;
    #else
      float64 eta[DIMENSION][DIMENSION];
      for(int j = 0; j < DIMENSION; j++) {
        ipos[j] = floor(xstar[j]);
        const float64 wx = xstar[j] - ipos[j];
        const float64 etax3 = (1./6.) * wx * wx * wx;
        const float64 etax0 = (1./6.) + 0.5 * wx * (wx - 1.) - etax3;
        const float64 etax2 = wx + etax0 - 2. * etax3;
        const float64 etax1 = 1. - etax0 - etax2 - etax3;
        eta[j][0] = etax0;
        eta[j][1] = etax1;
        eta[j][2] = etax2;
        eta[j][3] = etax3;
      }
      float64 ftmp1 = 0.;
      for(int k3 = 0; k3 <= 3; k3++) {
        for(int k2 = 0; k2 <= 3; k2++) {
          float64 ftmp2 = 0.;
          for(int k1 = 0; k1 <= 3; k1++) {
            float64 ftmp3 = 0.;
            for(int k0 = 0; k0 <= 3; k0++) {
              ftmp3 += eta[0][k0]
                     * tmp_fn(ipos[0]+k0-1, ipos[1]+k1-1, ipos[2]+k2-1, ipos[3]+k3-1);
            }
            ftmp2 += ftmp3 * eta[1][k1];
          }
          ftmp1 += ftmp2 * eta[2][k2] * eta[3][k3];
        }
      }
      return ftmp1;
    #endif
  }
static inline float64 interp_4D ( RealView4D & tmp_fn , float64 xmin [ DIMENSION ] , float64 inv_dx [ DIMENSION ] , float64 xstar [ DIMENSION ] ) { int ipos [ DIMENSION ] ; for ( int j = 0 ; j < DIMENSION ; j ++ ) { xstar [ j ] = inv_dx [ j ] * ( xstar [ j ] - xmin [ j ] ) ; } #ifdef LAG_ORDER float64 coefx0 [ LAG_PTS ] ; float64 coefx1 [ LAG_PTS ] ; float64 coefx2 [ LAG_PTS ] ; float64 coefx3 [ LAG_PTS ] ; ipos [ 0 ] = floor ( xstar [ 0 ] ) - LAG_OFFSET ; ipos [ 1 ] = floor ( xstar [ 1 ] ) - LAG_OFFSET ; ipos [ 2 ] = floor ( xstar [ 2 ] ) - LAG_OFFSET ; ipos [ 3 ] = floor ( xstar [ 3 ] ) - LAG_OFFSET ; lag_basis ( ( xstar [ 0 ] - ipos [ 0 ] ) , coefx0 ) ; lag_basis ( ( xstar [ 1 ] - ipos [ 1 ] ) , coefx1 ) ; lag_basis ( ( xstar [ 2 ] - ipos [ 2 ] ) , coefx2 ) ; lag_basis ( ( xstar [ 3 ] - ipos [ 3 ] ) , coefx3 ) ; float64 ftmp1 = 0. ; for ( int k3 = 0 ; k3 <= LAG_ORDER ; k3 ++ ) { for ( int k2 = 0 ; k2 <= LAG_ORDER ; k2 ++ ) { float64 ftmp2 = 0. ; for ( int k1 = 0 ; k1 <= LAG_ORDER ; k1 ++ ) { float64 ftmp3 = 0. ; for ( int k0 = 0 ; k0 <= LAG_ORDER ; k0 ++ ) { ftmp3 += coefx0 [ k0 ] * tmp_fn ( ipos [ 0 ] + k0 , ipos [ 1 ] + k1 , ipos [ 2 ] + k2 , ipos [ 3 ] + k3 ) ; } ftmp2 += ftmp3 * coefx1 [ k1 ] ; } ftmp1 += ftmp2 * coefx2 [ k2 ] * coefx3 [ k3 ] ; } } return ftmp1 ; #else float64 eta [ DIMENSION ] [ DIMENSION ] ; for ( int j = 0 ; j < DIMENSION ; j ++ ) { ipos [ j ] = floor ( xstar [ j ] ) ; const float64 wx = xstar [ j ] - ipos [ j ] ; const float64 etax3 = ( 1. / 6. ) * wx * wx * wx ; const float64 etax0 = ( 1. / 6. ) + 0.5 * wx * ( wx - 1. ) - etax3 ; const float64 etax2 = wx + etax0 - 2. * etax3 ; const float64 etax1 = 1. - etax0 - etax2 - etax3 ; eta [ j ] [ 0 ] = etax0 ; eta [ j ] [ 1 ] = etax1 ; eta [ j ] [ 2 ] = etax2 ; eta [ j ] [ 3 ] = etax3 ; } float64 ftmp1 = 0. ; for ( int k3 = 0 ; k3 <= 3 ; k3 ++ ) { for ( int k2 = 0 ; k2 <= 3 ; k2 ++ ) { float64 ftmp2 = 0. ; for ( int k1 = 0 ; k1 <= 3 ; k1 ++ ) { float64 ftmp3 = 0. ; for ( int k0 = 0 ; k0 <= 3 ; k0 ++ ) { ftmp3 += eta [ 0 ] [ k0 ] * tmp_fn ( ipos [ 0 ] + k0 - 1 , ipos [ 1 ] + k1 - 1 , ipos [ 2 ] + k2 - 1 , ipos [ 3 ] + k3 - 1 ) ; } ftmp2 += ftmp3 * eta [ 1 ] [ k1 ] ; } ftmp1 += ftmp2 * eta [ 2 ] [ k2 ] * eta [ 3 ] [ k3 ] ; } } return ftmp1 ; #endif }
Perform the interpolation of the 4d advection
Perform the interpolation of the 4d advection
//LAG_ORDER