msm_longrng.c File Reference

#include "msm_defn.h"

Go to the source code of this file.

Macros
#define	STENCIL_1D(_phi, _delta, _approx)
#define	D_STENCIL_1D(_dphi, _phi, _h_1, _delta, _approx)

Enumerations
enum	{ MAX_POLY_DEGREE = 9 }
enum	{ MAX_NSTENCIL = 11 }

Functions
static int	anterpolation (NL_Msm *)
static int	interpolation (NL_Msm *)
static int	restriction (NL_Msm *, int level)
static int	prolongation (NL_Msm *, int level)
static int	restriction_factored (NL_Msm *, int level)
static int	prolongation_factored (NL_Msm *, int level)
static int	gridcutoff (NL_Msm *, int level)
int	NL_msm_compute_long_range (NL_Msm *msm)

Variables
static const int	PolyDegree [NL_MSM_APPROX_END]
static const double	PhiStencilFactored [NL_MSM_APPROX_END][2 *MAX_POLY_DEGREE+1]
static const int	Nstencil [NL_MSM_APPROX_END]
static const int	IndexOffset [NL_MSM_APPROX_END][MAX_NSTENCIL]
static const double	PhiStencil [NL_MSM_APPROX_END][MAX_NSTENCIL]

Macro Definition Documentation

#define D_STENCIL_1D	(		_dphi,
			_phi,
			_h_1,
			_delta,
			_approx
	)

Calculate the stencil of basis function and derivatives of (1/h)Phi. dphi - stencil array (up to size MAX_POLY_DEGREE+1) phi - stencil array (up to size MAX_POLY_DEGREE+1) h_1 - 1/h, h is the grid spacing delta - normalized distance of atom from lowest grid point in stencil approx - APPROX enum value from msm.h

Definition at line 468 of file msm_longrng.c.

Referenced by interpolation().

#define STENCIL_1D	(		_phi,
			_delta,
			_approx
	)

Calculate the stencil of basis function values of Phi. phi - stencil array (up to size MAX_POLY_DEGREE+1) delta - normalized distance of atom from lowest grid point in stencil approx - APPROX enum value from msm.h

Definition at line 280 of file msm_longrng.c.

Referenced by anterpolation().

Enumeration Type Documentation

anonymous enum

Approximation formulaes are up to degree 9 polynomials.

Enumerator
MAX_POLY_DEGREE

Definition at line 159 of file msm_longrng.c.

{ MAX_POLY_DEGREE = 9 };

anonymous enum

Max stencil length is basically PolyDegree+2 for those approximations that interpolate. (We skip over zero in the complete stencils above.)

Enumerator
MAX_NSTENCIL

Definition at line 205 of file msm_longrng.c.

{ MAX_NSTENCIL = 11 };

Function Documentation

int anterpolation ( NL_Msm *  msm )

static

msm_longrng.c

Compute the long-range part of MSM forces.

The entire procedure can be depicted as an inverse V-cycle with bars across each level to represent the short-range gridded calculations. The algorithm first performs anterpolation of the charge to the finest level grid. Each intermediate grid level, below the coarsest level at the top, calculates a more slowly varying long-range part of the interaction by restricting the charges to a grid with twice the spacing. Each grid level also has a short-range cutoff part between charges that contributes to the potentials. Also summed to the potentials at each intermediate grid level are the long-range contributions prolongated from a grid with twice the spacing. Finally, the forces are interpolated (or approximated if the basis functions do not interpolate) from the finest level grid of potentials. The steps of the algorithm are detailed in section 2.3 of the thesis.

The restriction and prolongation operations use the factored approach having a work constant which is linear in the degree of the polynomial, i.e. O(pM), where p is the degree and M is the number of grid points.

The algorithm below uses stencil widths expressed as a function of the degree of the basis function polynomial Phi. The constants for restriction and prolongation stencils are stored in an array indexed by the APPROX enum defined in msm.h header. Calculations of the grid stencils of the 1D Phi and its derivative, needed for interpolation and anterpolation, are expressed as macros that use the APPROX enum to select the polynomial pieces from a switch statement. Periodic boundaries are handled by wrapping around the edges of the grid, as discussed in section 6.2 of the thesis.

XXX The short-range gridded calculations contribute to the virial, but it is not calculated.

XXX The factored restriction and prolongation procedures are not suitable for parallel and GPU computation due to the sequential dependence along each dimension of the grid.

Definition at line 779 of file msm_longrng.c.

References NL_Msm_t::approx, ASSERT, NL_Msm_t::atom, GRID_INDEX, GRID_INDEX_CHECK, GRID_ZERO, NL_Msm_t::gx, NL_Msm_t::gy, NL_Msm_t::gz, NL_Msm_t::hx, NL_Msm_t::hy, NL_Msm_t::hz, MAX_POLY_DEGREE, NL_Msm_t::msmflags, NL_MSM_ERROR_RANGE, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::numatoms, PolyDegree, NL_Msm_t::qh, and STENCIL_1D.

Referenced by NL_msm_compute_long_range().

{
  const double *atom = msm->atom;
  const int natoms = msm->numatoms;

  double xphi[MAX_POLY_DEGREE+1];  /* Phi stencil along x-dimension */
  double yphi[MAX_POLY_DEGREE+1];  /* Phi stencil along y-dimension */
  double zphi[MAX_POLY_DEGREE+1];  /* Phi stencil along z-dimension */
  double rx_hx, ry_hy, rz_hz;      /* normalized distance from atom to origin */
#ifndef TEST_INLINING
  double delta;                    /* normalized distance to stencil point */
#else
  double t;                    /* normalized distance to stencil point */
#endif
  double ck, cjk;
  const double hx_1 = 1/msm->hx;
  const double hy_1 = 1/msm->hy;
  const double hz_1 = 1/msm->hz;
  const double xm0 = msm->gx;      /* grid origin x */
  const double ym0 = msm->gy;      /* grid origin y */
  const double zm0 = msm->gz;      /* grid origin z */
  double q;

  NL_Msmgrid_double *qhgrid = &(msm->qh[0]);
  double *qh = qhgrid->data;
  const int ni = qhgrid->ni;
  const int nj = qhgrid->nj;
  const int nk = qhgrid->nk;
  const int ia = qhgrid->i0;
  const int ib = ia + ni - 1;
  const int ja = qhgrid->j0;
  const int jb = ja + nj - 1;
  const int ka = qhgrid->k0;
  const int kb = ka + nk - 1;

  const int ispx = (msm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (msm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (msm->msmflags & NL_MSM_PERIODIC_VEC3);
  int iswithin;

  int n, i, j, k, ilo, jlo, klo, koff;
  int jkoff, index;

  const int approx = msm->approx;
  const int s_size = PolyDegree[approx] + 1;         /* stencil size */
  const int s_edge = (PolyDegree[approx] - 1) >> 1;  /* stencil "edge" size */

  GRID_ZERO(qhgrid);

  for (n = 0;  n < natoms;  n++) {

    /* atomic charge */
    q = atom[4*n + 3];
    if (0==q) continue;

    /* distance between atom and origin measured in grid points */
    rx_hx = (atom[4*n    ] - xm0) * hx_1;
    ry_hy = (atom[4*n + 1] - ym0) * hy_1;
    rz_hz = (atom[4*n + 2] - zm0) * hz_1;

    /* find smallest numbered grid point in stencil */
    ilo = (int) floor(rx_hx) - s_edge;
    jlo = (int) floor(ry_hy) - s_edge;
    klo = (int) floor(rz_hz) - s_edge;

    /* calculate Phi stencils along each dimension */
#ifndef TEST_INLINING
    delta = rx_hx - (double) ilo;
    STENCIL_1D(xphi, delta, approx);
    delta = ry_hy - (double) jlo;
    STENCIL_1D(yphi, delta, approx);
    delta = rz_hz - (double) klo;
    STENCIL_1D(zphi, delta, approx);
#else
    t = rx_hx - (double) ilo;
        xphi[0] = 0.5 * (1 - t) * (2 - t) * (2 - t); \
        t--; \
        xphi[1] = (1 - t) * (1 + t - 1.5 * t * t); \
        t--; \
        xphi[2] = (1 + t) * (1 - t - 1.5 * t * t); \
        t--; \
        xphi[3] = 0.5 * (1 + t) * (2 + t) * (2 + t); \

    t = ry_hy - (double) jlo;
        yphi[0] = 0.5 * (1 - t) * (2 - t) * (2 - t); \
        t--; \
        yphi[1] = (1 - t) * (1 + t - 1.5 * t * t); \
        t--; \
        yphi[2] = (1 + t) * (1 - t - 1.5 * t * t); \
        t--; \
        yphi[3] = 0.5 * (1 + t) * (2 + t) * (2 + t); \

    t = rz_hz - (double) klo;
        zphi[0] = 0.5 * (1 - t) * (2 - t) * (2 - t); \
        t--; \
        zphi[1] = (1 - t) * (1 + t - 1.5 * t * t); \
        t--; \
        zphi[2] = (1 + t) * (1 - t - 1.5 * t * t); \
        t--; \
        zphi[3] = 0.5 * (1 + t) * (2 + t) * (2 + t); \

#endif

    /* test to see if stencil is within edges of grid */
    iswithin = ( ia <= ilo && (ilo+(s_size-1)) <= ib &&
                 ja <= jlo && (jlo+(s_size-1)) <= jb &&
                 ka <= klo && (klo+(s_size-1)) <= kb );

    if ( iswithin ) {  /* no wrapping needed */

      /* determine charge on cube of grid points around atom */
      for (k = 0;  k < s_size;  k++) {
        koff = (k + klo) * nj;
        ck = zphi[k] * q;
        for (j = 0;  j < s_size;  j++) {
          jkoff = (koff + (j + jlo)) * ni;
          cjk = yphi[j] * ck;
          for (i = 0;  i < s_size;  i++) {
            index = jkoff + (i + ilo);
            GRID_INDEX_CHECK(qhgrid, i+ilo, j+jlo, k+klo);
            ASSERT(GRID_INDEX(qhgrid, i+ilo, j+jlo, k+klo) == index);
            qh[index] += xphi[i] * cjk;
          }
        }
      }
    } /* if */

    else {  /* requires wrapping around grid */
      int ip, jp, kp;

      /* adjust ilo, jlo, klo so they are within grid indexing */
      if (ispx) {
        if      (ilo < ia) do { ilo += ni; } while (ilo < ia);
        else if (ilo > ib) do { ilo -= ni; } while (ilo > ib);
      }
      else if (ilo < ia || (ilo+(s_size-1)) > ib) {
        return NL_MSM_ERROR_RANGE;
      }

      if (ispy) {
        if      (jlo < ja) do { jlo += nj; } while (jlo < ja);
        else if (jlo > jb) do { jlo -= nj; } while (jlo > jb);
      }
      else if (jlo < ja || (jlo+(s_size-1)) > jb) {
        return NL_MSM_ERROR_RANGE;
      }

      if (ispz) {
        if      (klo < ka) do { klo += nk; } while (klo < ka);
        else if (klo > kb) do { klo -= nk; } while (klo > kb);
      }
      else if (klo < ka || (klo+(s_size-1)) > kb) {
        return NL_MSM_ERROR_RANGE;
      }

      /* determine charge on cube of grid points around atom, with wrapping */
      for (k = 0, kp = klo;  k < s_size;  k++, kp++) {
        if (kp > kb) kp = ka;  /* wrap stencil around grid */
        koff = kp * nj;
        ck = zphi[k] * q;
        for (j = 0, jp = jlo;  j < s_size;  j++, jp++) {
          if (jp > jb) jp = ja;  /* wrap stencil around grid */
          jkoff = (koff + jp) * ni;
          cjk = yphi[j] * ck;
          for (i = 0, ip = ilo;  i < s_size;  i++, ip++) {
            if (ip > ib) ip = ia;  /* wrap stencil around grid */
            index = jkoff + ip;
            GRID_INDEX_CHECK(qhgrid, ip, jp, kp);
            ASSERT(GRID_INDEX(qhgrid, ip, jp, kp) == index);
            qh[index] += xphi[i] * cjk;
          }
        }
      }
    } /* else */

  } /* end loop over atoms */

  return NL_MSM_SUCCESS;
} /* anterpolation */

int gridcutoff ( NL_Msm *  msm, int  level  )

static

Definition at line 1677 of file msm_longrng.c.

References ASSERT, NL_Msm_t::eh, NL_Msm_t::gc, GRID_INDEX, GRID_INDEX_CHECK, NL_Msm_t::msmflags, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, and NL_Msm_t::qh.

Referenced by NL_msm_compute_long_range().

{
  double eh_sum;

  /* grids of charge and potential */
  const NL_Msmgrid_double *qhgrid = &(msm->qh[level]);
  const double *qh = qhgrid->data;
  NL_Msmgrid_double *ehgrid = &(msm->eh[level]);
  double *eh = ehgrid->data;
  const int ia = qhgrid->i0;            /* lowest x-index */
  const int ib = ia + qhgrid->ni - 1;   /* highest x-index */
  const int ja = qhgrid->j0;            /* lowest y-index */
  const int jb = ja + qhgrid->nj - 1;   /* highest y-index */
  const int ka = qhgrid->k0;            /* lowest z-index */
  const int kb = ka + qhgrid->nk - 1;   /* highest z-index */
  const int ni = qhgrid->ni;            /* length along x-dim */
  const int nj = qhgrid->nj;            /* length along y-dim */
  const int nk = qhgrid->nk;            /* length along z-dim */

  /* grid of weights for pairwise grid point interactions within cutoff */
  const NL_Msmgrid_double *gcgrid = &(msm->gc[level]);
  const double *gc = gcgrid->data;
  const int gia = gcgrid->i0;            /* lowest x-index */
  const int gib = gia + gcgrid->ni - 1;  /* highest x-index */
  const int gja = gcgrid->j0;            /* lowest y-index */
  const int gjb = gja + gcgrid->nj - 1;  /* highest y-index */
  const int gka = gcgrid->k0;            /* lowest z-index */
  const int gkb = gka + gcgrid->nk - 1;  /* highest z-index */
  const int gni = gcgrid->ni;            /* length along x-dim */
  const int gnj = gcgrid->nj;            /* length along y-dim */

  const int ispx = (msm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (msm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (msm->msmflags & NL_MSM_PERIODIC_VEC3);

  const int ispnone = !(ispx || ispy || ispz);

  int i, j, k;
  int gia_clip, gib_clip;
  int gja_clip, gjb_clip;
  int gka_clip, gkb_clip;
  int koff;
  long jkoff, index;
  int id, jd, kd;
  int knoff;
  long jknoff, nindex;
  int kgoff, jkgoff, ngindex;

  if ( ispnone ) {  /* non-periodic boundaries */

    /* loop over all grid points */
    for (k = ka;  k <= kb;  k++) {

      /* clip gc ranges to keep offset for k index within grid */
      gka_clip = (k + gka < ka ? ka - k : gka);
      gkb_clip = (k + gkb > kb ? kb - k : gkb);

      koff = k * nj;  /* find eh flat index */

      for (j = ja;  j <= jb;  j++) {

        /* clip gc ranges to keep offset for j index within grid */
        gja_clip = (j + gja < ja ? ja - j : gja);
        gjb_clip = (j + gjb > jb ? jb - j : gjb);

        jkoff = (koff + j) * ni;  /* find eh flat index */

        for (i = ia;  i <= ib;  i++) {

          /* clip gc ranges to keep offset for i index within grid */
          gia_clip = (i + gia < ia ? ia - i : gia);
          gib_clip = (i + gib > ib ? ib - i : gib);

          index = jkoff + i;  /* eh flat index */

          /* sum over "sphere" of weighted charge */
          eh_sum = 0;
          for (kd = gka_clip;  kd <= gkb_clip;  kd++) {
            knoff = (k + kd) * nj;  /* find qh flat index */
            kgoff = kd * gnj;       /* find gc flat index */

            for (jd = gja_clip;  jd <= gjb_clip;  jd++) {
              jknoff = (knoff + (j + jd)) * ni;  /* find qh flat index */
              jkgoff = (kgoff + jd) * gni;       /* find gc flat index */

              for (id = gia_clip;  id <= gib_clip;  id++) {
                nindex = jknoff + (i + id);  /* qh flat index */
                ngindex = jkgoff + id;       /* gc flat index */

                GRID_INDEX_CHECK(qhgrid, i+id, j+jd, k+kd);
                ASSERT(GRID_INDEX(qhgrid, i+id, j+jd, k+kd) == nindex);

                GRID_INDEX_CHECK(gcgrid, id, jd, kd);
                ASSERT(GRID_INDEX(gcgrid, id, jd, kd) == ngindex);

                eh_sum += qh[nindex] * gc[ngindex];  /* sum weighted charge */
              }
            }
          } /* end loop over "sphere" of charge */

          GRID_INDEX_CHECK(ehgrid, i, j, k);
          ASSERT(GRID_INDEX(ehgrid, i, j, k) == index);
          eh[index] = eh_sum;  /* store potential */
        }
      }
    } /* end loop over all grid points */

  } /* if nonperiodic boundaries */
  else {
    /* some boundary is periodic */
    int ilo, jlo, klo;
    int ip, jp, kp;

    /* loop over all grid points */
    for (k = ka;  k <= kb;  k++) {
      klo = k + gka;
      if ( ! ispz ) {  /* nonperiodic z */
        /* clip gc ranges to keep offset for k index within grid */
        gka_clip = (k + gka < ka ? ka - k : gka);
        gkb_clip = (k + gkb > kb ? kb - k : gkb);
        if (klo < ka) klo = ka;  /* keep lowest qh index within grid */
      }
      else {  /* periodic z */
        gka_clip = gka;
        gkb_clip = gkb;
        if (klo < ka) do { klo += nk; } while (klo < ka);
      }
      ASSERT(klo <= kb);

      koff = k * nj;  /* find eh flat index */

      for (j = ja;  j <= jb;  j++) {
        jlo = j + gja;
        if ( ! ispy ) {  /* nonperiodic y */
          /* clip gc ranges to keep offset for j index within grid */
          gja_clip = (j + gja < ja ? ja - j : gja);
          gjb_clip = (j + gjb > jb ? jb - j : gjb);
          if (jlo < ja) jlo = ja;  /* keep lowest qh index within grid */
        }
        else {  /* periodic y */
          gja_clip = gja;
          gjb_clip = gjb;
          if (jlo < ja) do { jlo += nj; } while (jlo < ja);
        }
        ASSERT(jlo <= jb);

        jkoff = (koff + j) * ni;  /* find eh flat index */

        for (i = ia;  i <= ib;  i++) {
          ilo = i + gia;
          if ( ! ispx ) {  /* nonperiodic x */
            /* clip gc ranges to keep offset for i index within grid */
            gia_clip = (i + gia < ia ? ia - i : gia);
            gib_clip = (i + gib > ib ? ib - i : gib);
            if (ilo < ia) ilo = ia;  /* keep lowest qh index within grid */
          }
          else {  /* periodic x */
            gia_clip = gia;
            gib_clip = gib;
            if (ilo < ia) do { ilo += ni; } while (ilo < ia);
          }
          ASSERT(ilo <= ib);

          index = jkoff + i;  /* eh flat index */

          /* sum over "sphere" of weighted charge */
          eh_sum = 0;
          for (kd = gka_clip, kp = klo;  kd <= gkb_clip;  kd++, kp++) {
            /* clipping makes conditional always fail for nonperiodic */
            if (kp > kb) kp = ka;  /* wrap z direction */
            knoff = kp * nj;       /* find qh flat index */
            kgoff = kd * gnj;      /* find gc flat index */

            for (jd = gja_clip, jp = jlo;  jd <= gjb_clip;  jd++, jp++) {
              /* clipping makes conditional always fail for nonperiodic */
              if (jp > jb) jp = ja;              /* wrap y direction */
              jknoff = (knoff + jp) * ni;  /* find qh flat index */
              jkgoff = (kgoff + jd) * gni;       /* find gc flat index */

              for (id = gia_clip, ip = ilo;  id <= gib_clip;  id++, ip++) {
                /* clipping makes conditional always fail for nonperiodic */
                if (ip > ib) ip = ia;   /* wrap x direction */
                nindex = jknoff +  ip;  /* qh flat index */
                ngindex = jkgoff + id;  /* gc flat index */

                GRID_INDEX_CHECK(qhgrid, ip, jp, kp);
                ASSERT(GRID_INDEX(qhgrid, ip, jp, kp) == nindex);

                GRID_INDEX_CHECK(gcgrid, id, jd, kd);
                ASSERT(GRID_INDEX(gcgrid, id, jd, kd) == ngindex);

                eh_sum += qh[nindex] * gc[ngindex];  /* sum weighted charge */

              }
            }
          } /* end loop over "sphere" of charge */

          GRID_INDEX_CHECK(ehgrid, i, j, k);
          ASSERT(GRID_INDEX(ehgrid, i, j, k) == index);
          eh[index] = eh_sum;  /* store potential */
        }
      }
    } /* end loop over all grid points */

  } /* else some boundary is periodic */

  return NL_MSM_SUCCESS;
}

int interpolation ( NL_Msm *  msm )

static

Definition at line 960 of file msm_longrng.c.

References NL_Msm_t::approx, ASSERT, NL_Msm_t::atom, D_STENCIL_1D, NL_Msm_t::eh, NL_Msm_t::felec, GRID_INDEX, GRID_INDEX_CHECK, NL_Msm_t::gx, NL_Msm_t::gy, NL_Msm_t::gz, NL_Msm_t::gzero, NL_Msm_t::hx, NL_Msm_t::hy, NL_Msm_t::hz, MAX_POLY_DEGREE, NL_Msm_t::msmflags, NL_MSM_ERROR_RANGE, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::numatoms, PolyDegree, NL_Msm_t::qh, and NL_Msm_t::uelec.

Referenced by NL_msm_compute_long_range().

{
  double *f = msm->felec;
  const double *atom = msm->atom;
  const int natoms = msm->numatoms;

  double xphi[MAX_POLY_DEGREE+1];  /* Phi stencil along x-dimension */
  double yphi[MAX_POLY_DEGREE+1];  /* Phi stencil along y-dimension */
  double zphi[MAX_POLY_DEGREE+1];  /* Phi stencil along z-dimension */
  double dxphi[MAX_POLY_DEGREE+1]; /* derivative of Phi along x-dimension */
  double dyphi[MAX_POLY_DEGREE+1]; /* derivative of Phi along y-dimension */
  double dzphi[MAX_POLY_DEGREE+1]; /* derivative of Phi along z-dimension */
  double rx_hx, ry_hy, rz_hz;      /* normalized distance from atom to origin */
#ifndef TEST_INLINING
  double delta;                    /* normalized distance to stencil point */
#else
  double t;                    /* normalized distance to stencil point */
#endif
  const double hx_1 = 1/msm->hx;
  const double hy_1 = 1/msm->hy;
  const double hz_1 = 1/msm->hz;
  const double xm0 = msm->gx;      /* grid origin x */
  const double ym0 = msm->gy;      /* grid origin y */
  const double zm0 = msm->gz;      /* grid origin z */
  double q;
  double u = 0;

  const NL_Msmgrid_double *ehgrid = &(msm->eh[0]);
  const double *eh = ehgrid->data;
  const double *ebuf = ehgrid->buffer;
  const NL_Msmgrid_double *qhgrid = &(msm->qh[0]);
  const double *qbuf = qhgrid->buffer;
  const int ni = ehgrid->ni;
  const int nj = ehgrid->nj;
  const int nk = ehgrid->nk;
  const int ia = ehgrid->i0;
  const int ib = ia + ni - 1;
  const int ja = ehgrid->j0;
  const int jb = ja + nj - 1;
  const int ka = ehgrid->k0;
  const int kb = ka + nk - 1;

  const int ispx = (msm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (msm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (msm->msmflags & NL_MSM_PERIODIC_VEC3);
  int iswithin;

  double fx, fy, fz, cx, cy, cz;

  int n, i, j, k, ilo, jlo, klo, koff;
  long jkoff, index;
  const int nn = (ni*nj) * nk;

  const int approx = msm->approx;
  const int s_size = PolyDegree[approx] + 1;         /* stencil size */
  const int s_edge = (PolyDegree[approx] - 1) >> 1;  /* stencil "edge" size */

  for (n = 0;  n < natoms;  n++) {

    /* atomic charge */
    q = atom[4*n + 3];
    if (0==q) continue;

    /* distance between atom and origin measured in grid points */
    rx_hx = (atom[4*n    ] - xm0) * hx_1;
    ry_hy = (atom[4*n + 1] - ym0) * hy_1;
    rz_hz = (atom[4*n + 2] - zm0) * hz_1;

    /* find smallest numbered grid point in stencil */
    ilo = (int) floor(rx_hx) - s_edge;
    jlo = (int) floor(ry_hy) - s_edge;
    klo = (int) floor(rz_hz) - s_edge;

    /* calculate Phi stencils and derivatives along each dimension */
#ifndef TEST_INLINING
    delta = rx_hx - (double) ilo;
    //STENCIL_1D(xphi, delta, approx);
    D_STENCIL_1D(dxphi, xphi, hx_1, delta, approx);
    delta = ry_hy - (double) jlo;
    //STENCIL_1D(yphi, delta, approx);
    D_STENCIL_1D(dyphi, yphi, hy_1, delta, approx);
    delta = rz_hz - (double) klo;
    //STENCIL_1D(zphi, delta, approx);
    D_STENCIL_1D(dzphi, zphi, hz_1, delta, approx);
#else
    t = rx_hx - (double) ilo;
        xphi[0] = 0.5 * (1 - t) * (2 - t) * (2 - t); \
        dxphi[0] = (1.5 * t - 2) * (2 - t) * hx_1; \
        t--; \
        xphi[1] = (1 - t) * (1 + t - 1.5 * t * t); \
        dxphi[1] = (-5 + 4.5 * t) * t * hx_1; \
        t--; \
        xphi[2] = (1 + t) * (1 - t - 1.5 * t * t); \
        dxphi[2] = (-5 - 4.5 * t) * t * hx_1; \
        t--; \
        xphi[3] = 0.5 * (1 + t) * (2 + t) * (2 + t); \
        dxphi[3] = (1.5 * t + 2) * (2 + t) * hx_1; \

    t = ry_hy - (double) jlo;
        yphi[0] = 0.5 * (1 - t) * (2 - t) * (2 - t); \
        dyphi[0] = (1.5 * t - 2) * (2 - t) * hy_1; \
        t--; \
        yphi[1] = (1 - t) * (1 + t - 1.5 * t * t); \
        dyphi[1] = (-5 + 4.5 * t) * t * hy_1; \
        t--; \
        yphi[2] = (1 + t) * (1 - t - 1.5 * t * t); \
        dyphi[2] = (-5 - 4.5 * t) * t * hy_1; \
        t--; \
        yphi[3] = 0.5 * (1 + t) * (2 + t) * (2 + t); \
        dyphi[3] = (1.5 * t + 2) * (2 + t) * hy_1; \

    t = rz_hz - (double) klo;
        zphi[0] = 0.5 * (1 - t) * (2 - t) * (2 - t); \
        dzphi[0] = (1.5 * t - 2) * (2 - t) * hz_1; \
        t--; \
        zphi[1] = (1 - t) * (1 + t - 1.5 * t * t); \
        dzphi[1] = (-5 + 4.5 * t) * t * hz_1; \
        t--; \
        zphi[2] = (1 + t) * (1 - t - 1.5 * t * t); \
        dzphi[2] = (-5 - 4.5 * t) * t * hz_1; \
        t--; \
        zphi[3] = 0.5 * (1 + t) * (2 + t) * (2 + t); \
        dzphi[3] = (1.5 * t + 2) * (2 + t) * hz_1; \

#endif

    /* test to see if stencil is within edges of grid */
    iswithin = ( ia <= ilo && (ilo+(s_size-1)) <= ib &&
                 ja <= jlo && (jlo+(s_size-1)) <= jb &&
                 ka <= klo && (klo+(s_size-1)) <= kb );

    if ( iswithin ) {  /* no wrapping needed */

      /* determine force on atom from cube of grid point potentials */
      fx = fy = fz = 0;
      for (k = 0;  k < s_size;  k++) {
        koff = (k + klo) * nj;
        for (j = 0;  j < s_size;  j++) {
          jkoff = (koff + (j + jlo)) * ni;
          cx = yphi[j] * zphi[k];
          cy = dyphi[j] * zphi[k];
          cz = yphi[j] * dzphi[k];
          for (i = 0;  i < s_size;  i++) {
            index = jkoff + (i + ilo);
            GRID_INDEX_CHECK(ehgrid, i+ilo, j+jlo, k+klo);
            ASSERT(GRID_INDEX(ehgrid, i+ilo, j+jlo, k+klo) == index);
            fx += eh[index] * dxphi[i] * cx;
            fy += eh[index] * xphi[i] * cy;
            fz += eh[index] * xphi[i] * cz;
          }
        }
      }
    } /* if */

    else {  /* requires wrapping around grid */
      int ip, jp, kp;

      /* adjust ilo, jlo, klo so they are within grid indexing */
      if (ispx) {
        if      (ilo < ia) do { ilo += ni; } while (ilo < ia);
        else if (ilo > ib) do { ilo -= ni; } while (ilo > ib);
      }
      else if (ilo < ia || (ilo+(s_size-1)) > ib) {
        return NL_MSM_ERROR_RANGE;
      }

      if (ispy) {
        if      (jlo < ja) do { jlo += nj; } while (jlo < ja);
        else if (jlo > jb) do { jlo -= nj; } while (jlo > jb);
      }
      else if (jlo < ja || (jlo+(s_size-1)) > jb) {
        return NL_MSM_ERROR_RANGE;
      }

      if (ispz) {
        if      (klo < ka) do { klo += nk; } while (klo < ka);
        else if (klo > kb) do { klo -= nk; } while (klo > kb);
      }
      else if (klo < ka || (klo+(s_size-1)) > kb) {
        return NL_MSM_ERROR_RANGE;
      }

      /* determine force on atom from cube of grid point potentials, wrapping */
      fx = fy = fz = 0;
      for (k = 0, kp = klo;  k < s_size;  k++, kp++) {
        if (kp > kb) kp = ka;  /* wrap stencil around grid */
        koff = kp * nj;
        for (j = 0, jp = jlo;  j < s_size;  j++, jp++) {
          if (jp > jb) jp = ja;  /* wrap stencil around grid */
          jkoff = (koff + jp) * ni;
          cx = yphi[j] * zphi[k];
          cy = dyphi[j] * zphi[k];
          cz = yphi[j] * dzphi[k];
          for (i = 0, ip = ilo;  i < s_size;  i++, ip++) {
            if (ip > ib) ip = ia;  /* wrap stencil around grid */
            index = jkoff + ip;
            GRID_INDEX_CHECK(ehgrid, ip, jp, kp);
            ASSERT(GRID_INDEX(ehgrid, ip, jp, kp) == index);
            fx += eh[index] * dxphi[i] * cx;
            fy += eh[index] * xphi[i] * cy;
            fz += eh[index] * xphi[i] * cz;
          }
        }
      }
    } /* else */

    /* update force */
    f[3*n  ] -= q * fx;
    f[3*n+1] -= q * fy;
    f[3*n+2] -= q * fz;
//    printf("force[%d] = %g %g %g\n", n, -q*fx, -q*fy, -q*fz);

  } /* end loop over atoms */

  /* compute potential, subtract self potential */
  u = 0;
  if (1) {
    for (n = 0;  n < natoms;  n++) {
      double q = atom[4*n + 3];
      u += q * q;
    }
    u *= -msm->gzero;
  }
  for (index = 0;  index < nn;  index++) {
    u += qbuf[index] * ebuf[index];
  }
  msm->uelec += 0.5 * u;
//  printf("MSM long-range potential:  %g\n", 0.5 * u);

  return NL_MSM_SUCCESS;
} /* interpolation() */

int NL_msm_compute_long_range

(

NL_Msm *

msm

)

Definition at line 52 of file msm_longrng.c.

References anterpolation(), gridcutoff(), interpolation(), NL_Msm_t::msmflags, NL_MSM_COMPUTE_CUDA_FALL_BACK, NL_MSM_COMPUTE_CUDA_GRID_CUTOFF, NL_MSM_COMPUTE_NONFACTORED, NL_MSM_ERROR_SUPPORT, NL_MSM_SUCCESS, NL_Msm_t::nlevels, prolongation(), prolongation_factored(), NL_Msm_t::report_timings, restriction(), restriction_factored(), NL_Msm_t::timer_longrng, wkf_timer_start(), wkf_timer_stop(), and wkf_timer_time().

Referenced by NL_msm_compute_force().

                                           {
  double time_delta;
  int rc = 0;
  int level;
  int nlevels = msm->nlevels;
  int use_cuda_gridcutoff = (msm->msmflags & NL_MSM_COMPUTE_CUDA_GRID_CUTOFF);
  int fallback_cpu = (msm->msmflags & NL_MSM_COMPUTE_CUDA_FALL_BACK);
  int use_nonfactored = (msm->msmflags & NL_MSM_COMPUTE_NONFACTORED);

  wkf_timer_start(msm->timer_longrng);
  rc = anterpolation(msm);
  if (rc) return rc;
  wkf_timer_stop(msm->timer_longrng);
  time_delta = wkf_timer_time(msm->timer_longrng);
  if (msm->report_timings) {
    printf("MSM anterpolation:  %6.3f sec\n", time_delta);
  }

  wkf_timer_start(msm->timer_longrng);
  if (use_nonfactored) {
    for (level = 0;  level < nlevels - 1;  level++) {
      rc = restriction(msm, level);
      if (rc) return rc;
    }
  }
  else {
    for (level = 0;  level < nlevels - 1;  level++) {
      rc = restriction_factored(msm, level);
      if (rc) return rc;
    }
  }
  wkf_timer_stop(msm->timer_longrng);
  time_delta = wkf_timer_time(msm->timer_longrng);
  if (msm->report_timings) {
    printf("MSM restriction:    %6.3f sec\n", time_delta);
  }

  wkf_timer_start(msm->timer_longrng);
  if (use_cuda_gridcutoff && nlevels > 1) {
#ifdef NL_USE_CUDA
    if ((rc = NL_msm_cuda_condense_qgrids(msm)) != NL_MSM_SUCCESS ||
        (rc = NL_msm_cuda_compute_gridcutoff(msm)) != NL_MSM_SUCCESS ||
        (rc = NL_msm_cuda_expand_egrids(msm)) != NL_MSM_SUCCESS) {
      if (fallback_cpu) {
        printf("Falling back on CPU for grid cutoff computation\n");
        use_cuda_gridcutoff = 0;
      }
      else return rc;
    }
#else
    if (fallback_cpu) {
      printf("Falling back on CPU for grid cutoff computation\n");
      use_cuda_gridcutoff = 0;
    }
    else return NL_MSM_ERROR_SUPPORT;
#endif
  }

  if ( ! use_cuda_gridcutoff ) {
    for (level = 0;  level < nlevels - 1;  level++) {
      rc = gridcutoff(msm, level);
      if (rc) return rc;
    }
  }

  rc = gridcutoff(msm, level);  /* top level */
  if (rc) return rc;

  wkf_timer_stop(msm->timer_longrng);
  time_delta = wkf_timer_time(msm->timer_longrng);
  if (msm->report_timings) {
    printf("MSM grid cutoff:    %6.3f sec\n", time_delta);
  }

  wkf_timer_start(msm->timer_longrng);
  if (use_nonfactored) {
    for (level--;  level >= 0;  level--) {
      rc = prolongation(msm, level);
      if (rc) return rc;
    }
  }
  else {
    for (level--;  level >= 0;  level--) {
      rc = prolongation_factored(msm, level);
      if (rc) return rc;
    }
  }
  wkf_timer_stop(msm->timer_longrng);
  time_delta = wkf_timer_time(msm->timer_longrng);
  if (msm->report_timings) {
    printf("MSM prolongation:   %6.3f sec\n", time_delta);
  }

  wkf_timer_start(msm->timer_longrng);
  rc = interpolation(msm);
  if (rc) return rc;
  wkf_timer_stop(msm->timer_longrng);
  time_delta = wkf_timer_time(msm->timer_longrng);
  if (msm->report_timings) {
    printf("MSM interpolation:  %6.3f sec\n", time_delta);
  }

  return 0;
}

int prolongation ( NL_Msm *  msm, int  level  )

static

Definition at line 1582 of file msm_longrng.c.

References NL_Msm_t::approx, NL_Msm_t::eh, IndexOffset, NL_Msm_t::msmflags, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Nstencil, and PhiStencil.

Referenced by MsmBlock::addPotential(), MsmC1HermiteBlock::addPotential(), and NL_msm_compute_long_range().

                                         {
  /* grids of charge, finer grid and coarser grid */
  NL_Msmgrid_double *ehgrid = &(msm->eh[level]);
  double *eh = ehgrid->data;        /* index the offset data buffer */
  const NL_Msmgrid_double *e2hgrid = &(msm->eh[level+1]);
  const double *e2h_buffer = e2hgrid->buffer;   /* index the raw buffer */

  /* finer grid index ranges and dimensions */
  const int ia1 = ehgrid->i0;             /* lowest x-index */
  const int ib1 = ia1 + ehgrid->ni - 1;   /* highest x-index */
  const int ja1 = ehgrid->j0;             /* lowest y-index */
  const int jb1 = ja1 + ehgrid->nj - 1;   /* highest y-index */
  const int ka1 = ehgrid->k0;             /* lowest z-index */
  const int kb1 = ka1 + ehgrid->nk - 1;   /* highest z-index */
  const int ni1 = ehgrid->ni;             /* length along x-dim */
  const int nj1 = ehgrid->nj;             /* length along y-dim */
  const int nk1 = ehgrid->nk;             /* length along z-dim */

  /* coarser grid index ranges and dimensions */
  const int ia2 = e2hgrid->i0;            /* lowest x-index */
  const int ib2 = ia2 + e2hgrid->ni - 1;  /* highest x-index */
  const int ja2 = e2hgrid->j0;            /* lowest y-index */
  const int jb2 = ja2 + e2hgrid->nj - 1;  /* highest y-index */
  const int ka2 = e2hgrid->k0;            /* lowest z-index */
  const int kb2 = ka2 + e2hgrid->nk - 1;  /* highest z-index */

  const int ispx = (msm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (msm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (msm->msmflags & NL_MSM_PERIODIC_VEC3);

  const int nstencil = Nstencil[msm->approx];
  const int *offset = IndexOffset[msm->approx];
  const double *phi = PhiStencil[msm->approx];

  double cjk;

  int i1, j1, k1, index1, jk1off, k1off;
  int i2, j2, k2;
  int index2;
  int i, j, k;

  for (index2 = 0, k2 = ka2;  k2 <= kb2;  k2++) {
    k1 = k2 * 2;
    for (j2 = ja2;  j2 <= jb2;  j2++) {
      j1 = j2 * 2;
      for (i2 = ia2;  i2 <= ib2;  i2++, index2++) {
        i1 = i2 * 2;

        for (k = 0;  k < nstencil;  k++) {
          k1off = k1 + offset[k];
          if (k1off < ka1) {
            if (ispz) do { k1off += nk1; } while (k1off < ka1);
            else continue;
          }
          else if (k1off > kb1) {
            if (ispz) do { k1off -= nk1; } while (k1off > kb1);
            else break;
          }
          k1off *= nj1;
          for (j = 0;  j < nstencil;  j++) {
            jk1off = j1 + offset[j];
            if (jk1off < ja1) {
              if (ispy) do { jk1off += nj1; } while (jk1off < ja1);
              else continue;
            }
            else if (jk1off > jb1) {
              if (ispy) do { jk1off -= nj1; } while (jk1off > jb1);
              else break;
            }
            jk1off = (k1off + jk1off) * ni1;
            cjk = phi[j] * phi[k];
            for (i = 0;  i < nstencil;  i++) {
              index1 = i1 + offset[i];
              if (index1 < ia1) {
                if (ispx) do { index1 += ni1; } while (index1 < ia1);
                else continue;
              }
              else if (index1 > ib1) {
                if (ispx) do { index1 -= ni1; } while (index1 > ib1);
                else break;
              }
              index1 += jk1off;
              eh[index1] += e2h_buffer[index2] * phi[i] * cjk;
            }
          }
        } /* end loop over finer grid stencil */

      }
    }
  } /* end loop over each coarser grid point */

  return NL_MSM_SUCCESS;
}

int prolongation_factored ( NL_Msm *  msm, int  level  )

static

Definition at line 1340 of file msm_longrng.c.

References NL_Msm_t::approx, NL_Msm_t::eh, NL_Msm_t::lyzd, NL_Msm_t::lzd, NL_Msm_t::msmflags, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, PhiStencilFactored, and PolyDegree.

Referenced by NL_msm_compute_long_range().

                                                  {
  /* grids of potential, finer grid and coarser grid */
  NL_Msmgrid_double *ehgrid = &(msm->eh[level]);
  double *eh = ehgrid->data;
  const NL_Msmgrid_double *e2hgrid = &(msm->eh[level+1]);
  const double *e2h = e2hgrid->data;

  /* finer grid index ranges and dimensions */
  const int ia1 = ehgrid->i0;             /* lowest x-index */
  const int ib1 = ia1 + ehgrid->ni - 1;   /* highest x-index */
  const int ja1 = ehgrid->j0;             /* lowest y-index */
  const int jb1 = ja1 + ehgrid->nj - 1;   /* highest y-index */
  const int ka1 = ehgrid->k0;             /* lowest z-index */
  const int kb1 = ka1 + ehgrid->nk - 1;   /* highest z-index */
  const int ni1 = ehgrid->ni;             /* length along x-dim */
  const int nj1 = ehgrid->nj;             /* length along y-dim */
  const int nk1 = ehgrid->nk;             /* length along z-dim */

  /* coarser grid index ranges and dimensions */
  const int ia2 = e2hgrid->i0;            /* lowest x-index */
  const int ib2 = ia2 + e2hgrid->ni - 1;  /* highest x-index */
  const int ja2 = e2hgrid->j0;            /* lowest y-index */
  const int jb2 = ja2 + e2hgrid->nj - 1;  /* highest y-index */
  const int ka2 = e2hgrid->k0;            /* lowest z-index */
  const int kb2 = ka2 + e2hgrid->nk - 1;  /* highest z-index */
  const int nrow_e2 = e2hgrid->ni;
  const int nstride_e2 = nrow_e2 * e2hgrid->nj;

  const int ispx = (msm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (msm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (msm->msmflags & NL_MSM_PERIODIC_VEC3);

  /* set buffer using indexing offset, so that indexing matches eh grid */
  double *ezd = msm->lzd + (-ka1);
  double *eyzd = msm->lyzd + (-ka1*nj1 + -ja1);

  const size_t size_lzd = nk1 * sizeof(double);
  const size_t size_lyzd = nj1 * nk1 * sizeof(double);

  const double *phi = NULL;

  int i2, j2, k2;
  int im, jm, km;
  int i, j, k;
  int index_plane_e2, index_e2;
  int index_jk, offset_k;
  int offset;

  const double *phi_factored = PhiStencilFactored[msm->approx];
  const int r_stencil = PolyDegree[msm->approx];  /* "radius" of stencil */
  const int diam_stencil = 2*r_stencil + 1;       /* "diameter" of stencil */

  for (i2 = ia2;  i2 <= ib2;  i2++) {
    memset(msm->lyzd, 0, size_lyzd);

    for (j2 = ja2;  j2 <= jb2;  j2++) {
      memset(msm->lzd, 0, size_lzd);
      index_plane_e2 = j2 * nrow_e2 + i2;

      for (k2 = ka2;  k2 <= kb2;  k2++) {
        index_e2 = k2 * nstride_e2 + index_plane_e2;
        km = (k2 << 1);  /* = 2*k2 */
        if ( ! ispz ) {  /* nonperiodic */
          int lower = km - r_stencil;
          int upper = km + r_stencil;
          if (lower < ka1) lower = ka1;
          if (upper > kb1) upper = kb1;  /* clip edges */
          phi = phi_factored + r_stencil;  /* center of stencil */
          for (k = lower;  k <= upper;  k++) {
            ezd[k] += phi[k-km] * e2h[index_e2];
          }
        }
        else {  /* periodic */
          int kp = km - r_stencil;  /* index at left end of stencil */
          if (kp < ka1) do { kp += nk1; } while (kp < ka1);  /* start inside */
          phi = phi_factored;  /* left end of stencil */
          for (k = 0;  k < diam_stencil;  k++, kp++) {
            if (kp > kb1) kp = ka1;  /* wrap around edge of grid */
            ezd[kp] += phi[k] * e2h[index_e2];
          }
        }
      } /* for k2 */

      for (k = ka1;  k <= kb1;  k++) {
        offset = k * nj1;
        jm = (j2 << 1);  /* = 2*j2 */
        if ( ! ispy ) {  /* nonperiodic */
          int lower = jm - r_stencil;
          int upper = jm + r_stencil;
          if (lower < ja1) lower = ja1;
          if (upper > jb1) upper = jb1;  /* clip edges */
          phi = phi_factored + r_stencil;  /* center of stencil */
          for (j = lower;  j <= upper;  j++) {
            eyzd[offset + j] += phi[j-jm] * ezd[k];
          }
        }
        else {  /* periodic */
          int jp = jm - r_stencil;  /* index at left end of stencil */
          if (jp < ja1) do { jp += nj1; } while (jp < ja1);  /* start inside */
          phi = phi_factored;  /* left end of stencil */
          for (j = 0;  j < diam_stencil;  j++, jp++) {
            if (jp > jb1) jp = ja1;  /* wrap around edge of grid */
            eyzd[offset + jp] += phi[j] * ezd[k];
          }
        }
      } /* for k */

    } /* for j2 */

    for (k = ka1;  k <= kb1;  k++) {
      offset_k = k * nj1;

      for (j = ja1;  j <= jb1;  j++) {
        index_jk = offset_k + j;
        offset = index_jk * ni1;
        im = (i2 << 1);  /* = 2*i2 */
        if ( ! ispx ) {  /* nonperiodic */
          int lower = im - r_stencil;
          int upper = im + r_stencil;
          if (lower < ia1) lower = ia1;
          if (upper > ib1) upper = ib1;  /* clip edges */
          phi = phi_factored + r_stencil;  /* center of stencil */
          for (i = lower;  i <= upper;  i++) {
            eh[offset + i] += phi[i-im] * eyzd[index_jk];
          }
        }
        else {  /* periodic */
          int ip = im - r_stencil;  /* index at left end of stencil */
          if (ip < ia1) do { ip += ni1; } while (ip < ia1);  /* start inside */
          phi = phi_factored;  /* left end of stencil */
          for (i = 0;  i < diam_stencil;  i++, ip++) {
            if (ip > ib1) ip = ia1;  /* wrap around edge of grid */
            eh[offset + ip] += phi[i] * eyzd[index_jk];
          }
        }
      } /* for j */

    } /* for k */

  } /* for i2 */

  return NL_MSM_SUCCESS;
} /* prolongation_factored */

int restriction ( NL_Msm *  msm, int  level  )

static

Definition at line 1485 of file msm_longrng.c.

References NL_Msm_t::approx, IndexOffset, NL_Msm_t::msmflags, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Nstencil, PhiStencil, and NL_Msm_t::qh.

Referenced by MsmBlock::addCharge(), MsmC1HermiteBlock::addCharge(), and NL_msm_compute_long_range().

                                        {
  /* grids of charge, finer grid and coarser grid */
  const NL_Msmgrid_double *qhgrid = &(msm->qh[level]);
  const double *qh = qhgrid->data;        /* index the offset data buffer */
  NL_Msmgrid_double *q2hgrid = &(msm->qh[level+1]);
  double *q2h_buffer = q2hgrid->buffer;   /* index the raw buffer */

  /* finer grid index ranges and dimensions */
  const int ia1 = qhgrid->i0;             /* lowest x-index */
  const int ib1 = ia1 + qhgrid->ni - 1;   /* highest x-index */
  const int ja1 = qhgrid->j0;             /* lowest y-index */
  const int jb1 = ja1 + qhgrid->nj - 1;   /* highest y-index */
  const int ka1 = qhgrid->k0;             /* lowest z-index */
  const int kb1 = ka1 + qhgrid->nk - 1;   /* highest z-index */
  const int ni1 = qhgrid->ni;             /* length along x-dim */
  const int nj1 = qhgrid->nj;             /* length along y-dim */
  const int nk1 = qhgrid->nk;             /* length along z-dim */

  /* coarser grid index ranges and dimensions */
  const int ia2 = q2hgrid->i0;            /* lowest x-index */
  const int ib2 = ia2 + q2hgrid->ni - 1;  /* highest x-index */
  const int ja2 = q2hgrid->j0;            /* lowest y-index */
  const int jb2 = ja2 + q2hgrid->nj - 1;  /* highest y-index */
  const int ka2 = q2hgrid->k0;            /* lowest z-index */
  const int kb2 = ka2 + q2hgrid->nk - 1;  /* highest z-index */

  const int ispx = (msm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (msm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (msm->msmflags & NL_MSM_PERIODIC_VEC3);

  const int nstencil = Nstencil[msm->approx];
  const int *offset = IndexOffset[msm->approx];
  const double *phi = PhiStencil[msm->approx];

  double q2h_sum, cjk;

  int i1, j1, k1, index1, jk1off, k1off;
  int i2, j2, k2;
  int index2;
  int i, j, k;

  for (index2 = 0, k2 = ka2;  k2 <= kb2;  k2++) {
    k1 = k2 * 2;
    for (j2 = ja2;  j2 <= jb2;  j2++) {
      j1 = j2 * 2;
      for (i2 = ia2;  i2 <= ib2;  i2++, index2++) {
        i1 = i2 * 2;

        q2h_sum = 0;
        for (k = 0;  k < nstencil;  k++) {
          k1off = k1 + offset[k];
          if (k1off < ka1) {
            if (ispz) do { k1off += nk1; } while (k1off < ka1);
            else continue;
          }
          else if (k1off > kb1) {
            if (ispz) do { k1off -= nk1; } while (k1off > kb1);
            else break;
          }
          k1off *= nj1;
          for (j = 0;  j < nstencil;  j++) {
            jk1off = j1 + offset[j];
            if (jk1off < ja1) {
              if (ispy) do { jk1off += nj1; } while (jk1off < ja1);
              else continue;
            }
            else if (jk1off > jb1) {
              if (ispy) do { jk1off -= nj1; } while (jk1off > jb1);
              else break;
            }
            jk1off = (k1off + jk1off) * ni1;
            cjk = phi[j] * phi[k];
            for (i = 0;  i < nstencil;  i++) {
              index1 = i1 + offset[i];
              if (index1 < ia1) {
                if (ispx) do { index1 += ni1; } while (index1 < ia1);
                else continue;
              }
              else if (index1 > ib1) {
                if (ispx) do { index1 -= ni1; } while (index1 > ib1);
                else break;
              }
              index1 += jk1off;
              q2h_sum += qh[index1] * phi[i] * cjk;
            }
          }
        } /* end loop over finer grid stencil */

        q2h_buffer[index2] = q2h_sum;  /* store charge to coarser grid */
      }
    }
  } /* end loop over each coarser grid point */

  return NL_MSM_SUCCESS;
}

int restriction_factored ( NL_Msm *  msm, int  level  )

static

Definition at line 1193 of file msm_longrng.c.

References NL_Msm_t::approx, NL_Msm_t::lyzd, NL_Msm_t::lzd, NL_Msm_t::msmflags, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, PhiStencilFactored, PolyDegree, and NL_Msm_t::qh.

Referenced by NL_msm_compute_long_range().

                                                 {
  /* grids of charge, finer grid and coarser grid */
  const NL_Msmgrid_double *qhgrid = &(msm->qh[level]);
  const double *qh = qhgrid->data;
  NL_Msmgrid_double *q2hgrid = &(msm->qh[level+1]);
  double *q2h = q2hgrid->data;

  /* finer grid index ranges and dimensions */
  const int ia1 = qhgrid->i0;             /* lowest x-index */
  const int ib1 = ia1 + qhgrid->ni - 1;   /* highest x-index */
  const int ja1 = qhgrid->j0;             /* lowest y-index */
  const int jb1 = ja1 + qhgrid->nj - 1;   /* highest y-index */
  const int ka1 = qhgrid->k0;             /* lowest z-index */
  const int kb1 = ka1 + qhgrid->nk - 1;   /* highest z-index */
  const int ni1 = qhgrid->ni;             /* length along x-dim */
  const int nj1 = qhgrid->nj;             /* length along y-dim */
  const int nk1 = qhgrid->nk;             /* length along z-dim */

  /* coarser grid index ranges and dimensions */
  const int ia2 = q2hgrid->i0;            /* lowest x-index */
  const int ib2 = ia2 + q2hgrid->ni - 1;  /* highest x-index */
  const int ja2 = q2hgrid->j0;            /* lowest y-index */
  const int jb2 = ja2 + q2hgrid->nj - 1;  /* highest y-index */
  const int ka2 = q2hgrid->k0;            /* lowest z-index */
  const int kb2 = ka2 + q2hgrid->nk - 1;  /* highest z-index */
  const int nrow_q2 = q2hgrid->ni;
  const int nstride_q2 = nrow_q2 * q2hgrid->nj;

  const int ispx = (msm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (msm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (msm->msmflags & NL_MSM_PERIODIC_VEC3);

  /* set buffer using indexing offset, so that indexing matches qh grid */
  double *qzd = msm->lzd + (-ka1);
  double *qyzd = msm->lyzd + (-ka1*nj1 + -ja1);
  double qsum;

  const double *phi = NULL;

  int i2, j2, k2;
  int im, jm, km;
  int i, j, k;
  int index_plane_q2, index_q2;
  int index_jk, offset_k;
  int offset;

  const double *phi_factored = PhiStencilFactored[msm->approx];
  const int r_stencil = PolyDegree[msm->approx];  /* "radius" of stencil */
  const int diam_stencil = 2*r_stencil + 1;       /* "diameter" of stencil */

  for (i2 = ia2;  i2 <= ib2;  i2++) {

    for (k = ka1;  k <= kb1;  k++) {
      offset_k = k * nj1;

      for (j = ja1;  j <= jb1;  j++) {
        index_jk = offset_k + j;
        offset = index_jk * ni1;
        im = (i2 << 1);  /* = 2*i2 */
        qsum = 0;
        if ( ! ispx ) {  /* nonperiodic */
          int lower = im - r_stencil;
          int upper = im + r_stencil;
          if (lower < ia1) lower = ia1;
          if (upper > ib1) upper = ib1;  /* clip edges */
          phi = phi_factored + r_stencil;  /* center of stencil */
          for (i = lower;  i <= upper;  i++) {
            qsum += phi[i-im] * qh[offset + i];
          }
        }
        else {  /* periodic */
          int ip = im - r_stencil;  /* index at left end of stencil */
          if (ip < ia1) do { ip += ni1; } while (ip < ia1);  /* start inside */
          phi = phi_factored;  /* left end of stencil */
          for (i = 0;  i < diam_stencil;  i++, ip++) {
            if (ip > ib1) ip = ia1;  /* wrap around edge of grid */
            qsum += phi[i] * qh[offset + ip];
          }
        }
        qyzd[index_jk] = qsum;
      } /* for j */

    } /* for k */

    for (j2 = ja2;  j2 <= jb2;  j2++) {
      index_plane_q2 = j2 * nrow_q2 + i2;

      for (k = ka1;  k <= kb1;  k++) {
        offset = k * nj1;
        jm = (j2 << 1);  /* = 2*j2 */
        qsum = 0;
        if ( ! ispy ) {  /* nonperiodic */
          int lower = jm - r_stencil;
          int upper = jm + r_stencil;
          if (lower < ja1) lower = ja1;
          if (upper > jb1) upper = jb1;  /* clip edges */
          phi = phi_factored + r_stencil;  /* center of stencil */
          for (j = lower;  j <= upper;  j++) {
            qsum += phi[j-jm] * qyzd[offset + j];
          }
        }
        else {  /* periodic */
          int jp = jm - r_stencil;  /* index at left end of stencil */
          if (jp < ja1) do { jp += nj1; } while (jp < ja1);  /* start inside */
          phi = phi_factored;  /* left end of stencil */
          for (j = 0;  j < diam_stencil;  j++, jp++) {
            if (jp > jb1) jp = ja1;  /* wrap around edge of grid */
            qsum += phi[j] * qyzd[offset + jp];
          }
        }
        qzd[k] = qsum;
      } /* for k */

      for (k2 = ka2;  k2 <= kb2;  k2++) {
        index_q2 = k2 * nstride_q2 + index_plane_q2;
        km = (k2 << 1);  /* = 2*k2 */
        qsum = 0;
        if ( ! ispz ) {  /* nonperiodic */
          int lower = km - r_stencil;
          int upper = km + r_stencil;
          if (lower < ka1) lower = ka1;
          if (upper > kb1) upper = kb1;  /* clip edges */
          phi = phi_factored + r_stencil;  /* center of stencil */
          for (k = lower;  k <= upper;  k++) {
            qsum += phi[k-km] * qzd[k];
          }
        }
        else {  /* periodic */
          int kp = km - r_stencil;  /* index at left end of stencil */
          if (kp < ka1) do { kp += nk1; } while (kp < ka1);  /* start inside */
          phi = phi_factored;  /* left end of stencil */
          for (k = 0;  k < diam_stencil;  k++, kp++) {
            if (kp > kb1) kp = ka1;  /* wrap around edge of grid */
            qsum += phi[k] * qzd[kp];
          }
        }
        q2h[index_q2] = qsum;
      } /* for k2 */

    } /* for j2 */

  } /* for i2 */

  return NL_MSM_SUCCESS;
} /* restriction_factored */

Variable Documentation

const int IndexOffset[NL_MSM_APPROX_END][MAX_NSTENCIL]

static

Initial value:

= {
  
  {-3, -1, 0, 1, 3},

  
  {-5, -3, -1, 0, 1, 3, 5},

  
  {-5, -3, -1, 0, 1, 3, 5},

  
  {-7, -5, -3, -1, 0, 1, 3, 5, 7},

  
  {-7, -5, -3, -1, 0, 1, 3, 5, 7},

  
  {-9, -7, -5, -3, -1, 0, 1, 3, 5, 7, 9},

  
  {-9, -7, -5, -3, -1, 0, 1, 3, 5, 7, 9},

  
  {-3, -2, -1, 0, 1, 2, 3},
}

Index offsets from the stencil-centered grid element, to get to the correct contributing grid element.

Definition at line 214 of file msm_longrng.c.

Referenced by prolongation(), and restriction().

const int Nstencil[NL_MSM_APPROX_END]

static

Initial value:

= {
  5, 7, 7, 9, 9, 11, 11, 7
}

The stencil array lengths below.

Definition at line 208 of file msm_longrng.c.

Referenced by prolongation(), and restriction().

const double PhiStencil[NL_MSM_APPROX_END][MAX_NSTENCIL]

static

The grid transfer stencils for the non-factored restriction and prolongation procedures.

Definition at line 242 of file msm_longrng.c.

Referenced by prolongation(), and restriction().

const double PhiStencilFactored[NL_MSM_APPROX_END][2 *MAX_POLY_DEGREE+1]

static

The grid transfer stencils for factored restriction and prolongation. Must be listed in same order as APPROX enum from msm.h

Definition at line 170 of file msm_longrng.c.

Referenced by prolongation_factored(), and restriction_factored().

const int PolyDegree[NL_MSM_APPROX_END]

static

Initial value:

= {
  3, 5, 5, 7, 7, 9, 9, 3,
}

Degree of polynomial basis function Phi. Must be listed in same order as APPROX enum from msm.h

Definition at line 163 of file msm_longrng.c.

Referenced by anterpolation(), interpolation(), prolongation_factored(), and restriction_factored().