shiftField_lex_imp.cpp

Go to the documentation of this file.

#include "shiftField_lex.h"

#include <vector>

#include "bridgeIO.h"
using Bridge::vout;

//====================================================================
void ShiftField_lex::backward(Field& v, const Field& w,
                              const int mu)
{
  int bc = 1;

  if (mu == 0) {  // x-direction
    up_x(&v, &w, bc);
  } else if (mu == 1) {
    up_y(&v, &w, bc);
  } else if (mu == 2) {
    up_z(&v, &w, bc);
  } else if (mu == 3) {
    up_t(&v, &w, bc);
  } else {
    vout.crucial("Shift_Field_lex: parameter error, wrong parameter\n");
    exit(EXIT_FAILURE);
  }
}


//====================================================================
void ShiftField_lex::forward(Field& v, const Field& w,
                             const int mu)
{
  int bc = 1;

  if (mu == 0) {
    dn_x(&v, &w, bc);
  } else if (mu == 1) {
    dn_y(&v, &w, bc);
  } else if (mu == 2) {
    dn_z(&v, &w, bc);
  } else if (mu == 3) {
    dn_t(&v, &w, bc);
  } else {
    vout.crucial("Shift_Field_lex: parameter error, wrong parameter\n");
    exit(EXIT_FAILURE);
  }
}


//====================================================================
void ShiftField_lex::backward(Field& v, const Field& w,
                              const int bc, const int mu)
{
  if (mu == 0) {  // x-direction
    up_x(&v, &w, bc);
  } else if (mu == 1) {
    up_y(&v, &w, bc);
  } else if (mu == 2) {
    up_z(&v, &w, bc);
  } else if (mu == 3) {
    up_t(&v, &w, bc);
  } else {
    vout.crucial("Shift_Field_lex: parameter error, wrong parameter\n");
    exit(EXIT_FAILURE);
  }
}


//====================================================================
void ShiftField_lex::forward(Field& v, const Field& w,
                             const int bc, const int mu)
{
  if (mu == 0) {
    dn_x(&v, &w, bc);
  } else if (mu == 1) {
    dn_y(&v, &w, bc);
  } else if (mu == 2) {
    dn_z(&v, &w, bc);
  } else if (mu == 3) {
    dn_t(&v, &w, bc);
  } else {
    vout.crucial("Shift_Field_lex: parameter error, wrong parameter\n");
    exit(EXIT_FAILURE);
  }
}


//====================================================================
void ShiftField_lex::up_x(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(0) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Nx;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int s2 = 0; s2 < Ny * Nz * Nt; s2++) {
    // bulk
    for (int x = 0; x < Nx - 1; x++) {
      int ix = x + Nx * s2;
      int px = ix + 1;
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }

    // boundary (x=Nx-1)
    int px = Nx * s2;
    // int ix = Nx - 1 + px;
    for (int ex = 0; ex < Nex; ex++) {
      for (int in = 0; in < Nin; in++) {
        wt[in + Nin * (s2 + Nvol2 * ex)] = bc2 * wp[in + Nin * (px + Nvol * ex)];
      }
    }
  }

  int size = Nin * Nvol2 * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 0, 1, 0);

  for (int s2 = 0; s2 < Ny * Nz * Nt; s2++) {
    // boundary (x=Nx-1)
    int px = Nx * s2;
    int ix = Nx - 1 + px;
    for (int in = 0; in < Nin; in++) {
      for (int ex = 0; ex < Nex; ex++) {
        vp[in + Nin * (ix + Nvol * ex)] = vt[in + Nin * (s2 + Nvol2 * ex)];
      }
    }
  }
}


//====================================================================
void ShiftField_lex::dn_x(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(0) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Nx;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int s2 = 0; s2 < Ny * Nz * Nt; s2++) {
    // bulk
    for (int x = 1; x < Nx; x++) {
      int ix = x + Nx * s2;
      int px = ix - 1;
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }

    // boundary (x=0)
    int px = Nx - 1 + Nx * s2;
    for (int ex = 0; ex < Nex; ex++) {
      for (int in = 0; in < Nin; in++) {
        wt[in + Nin * (s2 + Nvol2 * ex)] = wp[in + Nin * (px + Nvol * ex)];
      }
    }
  }

  int size = Nin * Nvol2 * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 0, -1, 4);

  for (int s2 = 0; s2 < Ny * Nz * Nt; s2++) {
    // boundary (x=0)
    int ix = Nx * s2;
    for (int ex = 0; ex < Nex; ex++) {
      for (int in = 0; in < Nin; in++) {
        vp[in + Nin * (ix + Nvol * ex)] = bc2 * vt[in + Nin * (s2 + Nvol2 * ex)];
      }
    }
  }
}


//====================================================================
void ShiftField_lex::up_y(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(1) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Ny;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int zt = 0; zt < Nz * Nt; zt++) {
    for (int x = 0; x < Nx; x++) {
      int s1 = x + Nx * Ny * zt;
      int s2 = x + Nx * zt;

      // bulk
      for (int y = 0; y < Ny - 1; y++) {
        int ix = s1 + Nx * y;
        int px = ix + Nx;
        for (int ex = 0; ex < Nex; ex++) {
          for (int in = 0; in < Nin; in++) {
            vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
          }
        }
      }

      // boundary (y=Ny-1)
      int px = s1;
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          wt[in + Nin * (s2 + Nvol2 * ex)] = bc2 * wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }
  }

  int size = Nin * Nvol2 * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 1, 1, 1);

  for (int zt = 0; zt < Nz * Nt; zt++) {
    for (int x = 0; x < Nx; x++) {
      int s1 = x + Nx * Ny * zt;
      int s2 = x + Nx * zt;

      // boundary (y=Ny-1)
      int ix = s1 + Nx * (Ny - 1);
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          vp[in + Nin * (ix + Nvol * ex)] = vt[in + Nin * (s2 + Nvol2 * ex)];
        }
      }
    }
  }
}


//====================================================================
void ShiftField_lex::dn_y(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(1) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Ny;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int zt = 0; zt < Nz * Nt; zt++) {
    for (int x = 0; x < Nx; x++) {
      int s1 = x + Nx * Ny * zt;
      int s2 = x + Nx * zt;

      // bulk
      for (int y = 1; y < Ny; y++) {
        int ix = s1 + Nx * y;
        int px = ix - Nx;
        for (int ex = 0; ex < Nex; ex++) {
          for (int in = 0; in < Nin; in++) {
            vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
          }
        }
      }

      // boundary (y=0)
      int px = s1 + Nx * (Ny - 1);
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          wt[in + Nin * (s2 + Nvol2 * ex)] = wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }
  }

  int size = Nin * Nvol2 * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 1, -1, 5);

  for (int zt = 0; zt < Nz * Nt; zt++) {
    for (int x = 0; x < Nx; x++) {
      int s1 = x + Nx * Ny * zt;
      int s2 = x + Nx * zt;

      // boundary (y=0)
      int ix = s1;
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          vp[in + Nin * (ix + Nvol * ex)] = bc2 * vt[in + Nin * (s2 + Nvol2 * ex)];
        }
      }
    }
  }
}


//====================================================================
void ShiftField_lex::up_z(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(2) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Nz;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int t = 0; t < Nt; t++) {
    for (int xy = 0; xy < Nx * Ny; xy++) {
      int s1 = xy + Nx * Ny * Nz * t;
      int s2 = xy + Nx * Ny * t;

      // bulk
      for (int z = 0; z < Nz - 1; z++) {
        int ix = s1 + Nx * Ny * z;
        int px = ix + Nx * Ny;
        for (int ex = 0; ex < Nex; ex++) {
          for (int in = 0; in < Nin; in++) {
            vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
          }
        }
      }

      // boundary (z=Nz-1)
      int px = s1;
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          wt[in + Nin * (s2 + Nvol2 * ex)] = bc2 * wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }
  }

  int size = Nin * (Nvol / Nz) * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 2, 1, 2);

  for (int t = 0; t < Nt; t++) {
    for (int xy = 0; xy < Nx * Ny; xy++) {
      int s1 = xy + Nx * Ny * Nz * t;
      int s2 = xy + Nx * Ny * t;

      // boundary (z=Nz-1)
      int ix = s1 + Nx * Ny * (Nz - 1);
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          vp[in + Nin * (ix + Nvol * ex)] = vt[in + Nin * (s2 + Nvol2 * ex)];
        }
      }
    }
  }
}


//====================================================================
void ShiftField_lex::dn_z(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(2) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Nz;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int t = 0; t < Nt; t++) {
    for (int xy = 0; xy < Nx * Ny; xy++) {
      int s1 = xy + Nx * Ny * Nz * t;
      int s2 = xy + Nx * Ny * t;

      // bulk
      for (int z = 1; z < Nz; z++) {
        int ix = s1 + Nx * Ny * z;
        int px = ix - Nx * Ny;
        for (int ex = 0; ex < Nex; ex++) {
          for (int in = 0; in < Nin; in++) {
            vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
          }
        }
      }
      // boundary (z=0)
      int px = s1 + Nx * Ny * (Nz - 1);
      for (int in = 0; in < Nin; in++) {
        for (int ex = 0; ex < Nex; ex++) {
          wt[in + Nin * (s2 + Nvol2 * ex)] = wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }
  }

  int size = Nin * Nvol2 * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 2, -1, 6);

  for (int t = 0; t < Nt; t++) {
    for (int xy = 0; xy < Nx * Ny; xy++) {
      int s1 = xy + Nx * Ny * Nz * t;
      int s2 = xy + Nx * Ny * t;

      // boundary (z=0)
      int ix = s1;
      for (int in = 0; in < Nin; in++) {
        for (int ex = 0; ex < Nex; ex++) {
          vp[in + Nin * (ix + Nvol * ex)] = bc2 * vt[in + Nin * (s2 + Nvol2 * ex)];
        }
      }
    }
  }
}


//====================================================================
void ShiftField_lex::up_t(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(3) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Nt;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int s2 = 0; s2 < Nx * Ny * Nz; s2++) {
    // bulk
    for (int t = 0; t < Nt - 1; t++) {
      int ix = s2 + Nx * Ny * Nz * t;
      int px = ix + Nx * Ny * Nz;
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }

    // boundary (t=Nt-1)
    int px = s2;
    for (int ex = 0; ex < Nex; ex++) {
      for (int in = 0; in < Nin; in++) {
        wt[in + Nin * (s2 + Nvol2 * ex)] = bc2 * wp[in + Nin * (px + Nvol * ex)];
      }
    }
  }

  int size = Nin * Nvol2 * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 3, 1, 3);

  for (int s2 = 0; s2 < Nx * Ny * Nz; s2++) {
    // boundary (t=Nt-1)
    int ix = s2 + Nx * Ny * Nz * (Nt - 1);
    for (int ex = 0; ex < Nex; ex++) {
      for (int in = 0; in < Nin; in++) {
        vp[in + Nin * (ix + Nvol * ex)] = vt[in + Nin * (s2 + Nvol2 * ex)];
      }
    }
  }
}


//====================================================================
void ShiftField_lex::dn_t(Field *v, const Field *w, const int bc)
{
  double bc2;

  if (Communicator::ipe(3) == 0) {
    bc2 = bc;
  } else {
    bc2 = 1.0;
  }

  int Nin   = w->nin();
  int Nex   = w->nex();
  int Nvol  = w->nvol();
  int Nvol2 = Nvol / Nt;

  double       *vp = v->ptr(0);
  const double *wp = w->ptr(0);

  std::vector<double> wt(Nin * Nvol2 * Nex);
  std::vector<double> vt(Nin * Nvol2 * Nex);

  for (int s2 = 0; s2 < Nx * Ny * Nz; s2++) {
    // bulk
    for (int t = 1; t < Nt; t++) {
      int ix = s2 + Nx * Ny * Nz * t;
      int px = ix - Nx * Ny * Nz;
      for (int ex = 0; ex < Nex; ex++) {
        for (int in = 0; in < Nin; in++) {
          vp[in + Nin * (ix + Nvol * ex)] = wp[in + Nin * (px + Nvol * ex)];
        }
      }
    }
    // boundary (t=0)
    int px = s2 + Nx * Ny * Nz * (Nt - 1);
    for (int ex = 0; ex < Nex; ex++) {
      for (int in = 0; in < Nin; in++) {
        wt[in + Nin * (s2 + Nvol2 * ex)] = wp[in + Nin * (px + Nvol * ex)];
      }
    }
  }

  int size = Nin * Nvol2 * Nex;
  Communicator::exchange(size, &vt[0], &wt[0], 3, -1, 7);

  for (int s2 = 0; s2 < Nx * Ny * Nz; s2++) {
    // boundary (t=0)
    int ix = s2;
    for (int ex = 0; ex < Nex; ex++) {
      for (int in = 0; in < Nin; in++) {
        vp[in + Nin * (ix + Nvol * ex)] = bc2 * vt[in + Nin * (s2 + Nvol2 * ex)];
      }
    }
  }
}


//====================================================================
//============================================================END=====