msm_longrng_sprec.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_sprec (NL_Msm *msm)

Variables
static const int	PolyDegree [NL_MSM_APPROX_END]
static const float	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 float	PhiStencil [NL_MSM_APPROX_END][MAX_NSTENCIL]

Macro Definition Documentation

D_STENCIL_1D

#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_sprec.c.

Referenced by interpolation().

STENCIL_1D

#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_sprec.c.

Referenced by anterpolation().

Enumeration Type Documentation

anonymous enum

anonymous enum

Approximation formulaes are up to degree 9 polynomials.

Enumerator
MAX_POLY_DEGREE

Definition at line 159 of file msm_longrng_sprec.c.

{ MAX_POLY_DEGREE = 9 };

anonymous enum

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_sprec.c.

{ MAX_NSTENCIL = 11 };

Function Documentation

anterpolation()

int anterpolation ( NL_Msm *  msm )

static

msm_longrng.c

Compute the long-range part of MSM forces, using single precision.

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_sprec.c.

References ASSERT, GRID_INDEX, GRID_INDEX_CHECK, GRID_ZERO, MAX_POLY_DEGREE, NL_MSM_ERROR_RANGE, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, PolyDegree, and STENCIL_1D.

Referenced by NL_msm_compute_long_range_sprec().

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

  float xphi[MAX_POLY_DEGREE+1];  /* Phi stencil along x-dimension */
  float yphi[MAX_POLY_DEGREE+1];  /* Phi stencil along y-dimension */
  float zphi[MAX_POLY_DEGREE+1];  /* Phi stencil along z-dimension */
  float rx_hx, ry_hy, rz_hz;      /* normalized distance from atom to origin */
#ifndef TEST_INLINING
  float delta;                    /* normalized distance to stencil point */
#else
  float t;                    /* normalized distance to stencil point */
#endif
  float ck, cjk;
  const float hx_1 = 1/msm->hx_f;
  const float hy_1 = 1/msm->hy_f;
  const float hz_1 = 1/msm->hz_f;
  const float xm0 = msm->gx_f;      /* grid origin x */
  const float ym0 = msm->gy_f;      /* grid origin y */
  const float zm0 = msm->gz_f;      /* grid origin z */
  float q;

  NL_Msmgrid_float *qhgrid = &(msm->qh_f[0]);
  float *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) floorf(rx_hx) - s_edge;
    jlo = (int) floorf(ry_hy) - s_edge;
    klo = (int) floorf(rz_hz) - s_edge;

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

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

    t = rz_hz - (double) klo;
        zphi[0] = 0.5f * (1 - t) * (2 - t) * (2 - t); \
        t--; \
        zphi[1] = (1 - t) * (1 + t - 1.5f * t * t); \
        t--; \
        zphi[2] = (1 + t) * (1 - t - 1.5f * t * t); \
        t--; \
        zphi[3] = 0.5f * (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 */

gridcutoff()

int gridcutoff ( NL_Msm *  msm, int  level  )

static

Definition at line 1690 of file msm_longrng_sprec.c.

References ASSERT, GRID_INDEX, GRID_INDEX_CHECK, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, and NL_MSM_SUCCESS.

Referenced by NL_msm_compute_long_range_sprec().

{
  float eh_sum;

  /* grids of charge and potential */
  const NL_Msmgrid_float *qhgrid = &(msm->qh_f[level]);
  const float *qh = qhgrid->data;
  NL_Msmgrid_float *ehgrid = &(msm->eh_f[level]);
  float *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_float *gcgrid = &(msm->gc_f[level]);
  const float  *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;

#ifdef NL_DEBUG
  if (level==0) {
    const float *gcbuf = gcgrid->buffer;
    int index = 0;
    printf("+++ qh[0,0,0] = %g\n", qh[0]);
#if 0
    for (k = gka;  k <= gkb;  k++) {
      for (j = gja;  j <= gjb;  j++) {
        for (i = gia;  i <= gib;  i++, index++) {
          printf("=== gc[%d,%d,%d] = %g\n", i, j, k, gcbuf[index]);
        }
      }
    }
#endif
  }
#endif

  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;
}

interpolation()

int interpolation ( NL_Msm *  msm )

static

Definition at line 960 of file msm_longrng_sprec.c.

References ASSERT, D_STENCIL_1D, GRID_INDEX, GRID_INDEX_CHECK, MAX_POLY_DEGREE, NL_MSM_ERROR_RANGE, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, and PolyDegree.

Referenced by NL_msm_compute_long_range_sprec().

{
  float *f = msm->felec_f;
  const float *atom = msm->atom_f;
  const int natoms = msm->numatoms;

  float xphi[MAX_POLY_DEGREE+1];  /* Phi stencil along x-dimension */
  float yphi[MAX_POLY_DEGREE+1];  /* Phi stencil along y-dimension */
  float zphi[MAX_POLY_DEGREE+1];  /* Phi stencil along z-dimension */
  float dxphi[MAX_POLY_DEGREE+1]; /* derivative of Phi along x-dimension */
  float dyphi[MAX_POLY_DEGREE+1]; /* derivative of Phi along y-dimension */
  float dzphi[MAX_POLY_DEGREE+1]; /* derivative of Phi along z-dimension */
  float rx_hx, ry_hy, rz_hz;      /* normalized distance from atom to origin */
#ifndef TEST_INLINING
  float delta;                    /* normalized distance to stencil point */
#else
  float t;                    /* normalized distance to stencil point */
#endif
  const float hx_1 = 1/msm->hx_f;
  const float hy_1 = 1/msm->hy_f;
  const float hz_1 = 1/msm->hz_f;
  const float xm0 = msm->gx_f;      /* grid origin x */
  const float ym0 = msm->gy_f;      /* grid origin y */
  const float zm0 = msm->gz_f;      /* grid origin z */
  float q;
  double u = 0;  /* need double precision for accumulating potential */

  const NL_Msmgrid_float *ehgrid = &(msm->eh_f[0]);
  const float *eh = ehgrid->data;
  const float *ebuf = ehgrid->buffer;
  const NL_Msmgrid_float *qhgrid = &(msm->qh_f[0]);
  const float *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;

  float 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 */

#ifdef NL_DEBUG
#if 1
  printf("+++ eh[0,0,0] = %g\n", eh[0]);
#else
  n = 0;
  for (k = ka;  k <= kb;  k++) {
    for (j = ja;  j <= jb;  j++) {
      for (i = ia;  i <= ib;  i++, n++) {
        printf("=== eh[%d,%d,%d] = %g\n", i, j, k, ebuf[n]);
      }
    }
  }
#endif
#endif

  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) floorf(rx_hx) - s_edge;
    jlo = (int) floorf(ry_hy) - s_edge;
    klo = (int) floorf(rz_hz) - s_edge;

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

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

    t = rz_hz - (double) klo;
        zphi[0] = 0.5f * (1 - t) * (2 - t) * (2 - t); \
        dzphi[0] = (1.5f * t - 2) * (2 - t) * hz_1; \
        t--; \
        zphi[1] = (1 - t) * (1 + t - 1.5f * t * t); \
        dzphi[1] = (-5 + 4.5f * t) * t * hz_1; \
        t--; \
        zphi[2] = (1 + t) * (1 - t - 1.5f * t * t); \
        dzphi[2] = (-5 - 4.5f * t) * t * hz_1; \
        t--; \
        zphi[3] = 0.5f * (1 + t) * (2 + t) * (2 + t); \
        dzphi[3] = (1.5f * 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;

  } /* end loop over atoms */

  /* compute potential, subtract self potential */
  u = 0;
  if (1) {
    for (n = 0;  n < natoms;  n++) {
      float q = atom[4*n + 3];
      u += q * q;
    }
    u *= -msm->gzero_f;
  }
  for (index = 0;  index < nn;  index++) {
    u += qbuf[index] * ebuf[index];
  }
  msm->uelec += 0.5 * u;

  return NL_MSM_SUCCESS;
} /* interpolation() */

NL_msm_compute_long_range_sprec()

int NL_msm_compute_long_range_sprec

(

NL_Msm *

msm

)

Definition at line 52 of file msm_longrng_sprec.c.

References anterpolation(), gridcutoff(), interpolation(), NL_MSM_COMPUTE_CUDA_FALL_BACK, NL_MSM_COMPUTE_CUDA_GRID_CUTOFF, NL_MSM_COMPUTE_NONFACTORED, NL_MSM_ERROR_SUPPORT, NL_MSM_SUCCESS, prolongation(), prolongation_factored(), restriction(), restriction_factored(), wkf_timer_start(), wkf_timer_stop(), and wkf_timer_time().

Referenced by NL_msm_compute_force_sprec().

                                                 {
  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;
}

prolongation()

int prolongation ( NL_Msm *  msm, int  level  )

static

Definition at line 1595 of file msm_longrng_sprec.c.

References IndexOffset, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Nstencil, and PhiStencil.

Referenced by NL_msm_compute_long_range_sprec().

                                         {
  /* grids of charge, finer grid and coarser grid */
  NL_Msmgrid_float *ehgrid = &(msm->eh_f[level]);
  float *eh = ehgrid->data;        /* index the offset data buffer */
  const NL_Msmgrid_float *e2hgrid = &(msm->eh_f[level+1]);
  const float *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 float *phi = PhiStencil[msm->approx];

  float 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;
}

prolongation_factored()

int prolongation_factored ( NL_Msm *  msm, int  level  )

static

Definition at line 1353 of file msm_longrng_sprec.c.

References 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_sprec().

                                                  {
  /* grids of potential, finer grid and coarser grid */
  NL_Msmgrid_float *ehgrid = &(msm->eh_f[level]);
  float *eh = ehgrid->data;
  const NL_Msmgrid_float *e2hgrid = &(msm->eh_f[level+1]);
  const float *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 */
  float *ezd = msm->lzd_f + (-ka1);
  float *eyzd = msm->lyzd_f + (-ka1*nj1 + -ja1);

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

  const float *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 float *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_f, 0, size_lyzd);

    for (j2 = ja2;  j2 <= jb2;  j2++) {
      memset(msm->lzd_f, 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 */

restriction()

int restriction ( NL_Msm *  msm, int  level  )

static

Definition at line 1498 of file msm_longrng_sprec.c.

References IndexOffset, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, Nstencil, and PhiStencil.

Referenced by NL_msm_compute_long_range_sprec().

                                        {
  /* grids of charge, finer grid and coarser grid */
  const NL_Msmgrid_float *qhgrid = &(msm->qh_f[level]);
  const float *qh = qhgrid->data;        /* index the offset data buffer */
  NL_Msmgrid_float *q2hgrid = &(msm->qh_f[level+1]);
  float *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 float *phi = PhiStencil[msm->approx];

  float 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;
}

restriction_factored()

int restriction_factored ( NL_Msm *  msm, int  level  )

static

Definition at line 1206 of file msm_longrng_sprec.c.

References 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_sprec().

                                                 {
  /* grids of charge, finer grid and coarser grid */
  const NL_Msmgrid_float *qhgrid = &(msm->qh_f[level]);
  const float *qh = qhgrid->data;
  NL_Msmgrid_float *q2hgrid = &(msm->qh_f[level+1]);
  float *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 */
  float *qzd = msm->lzd_f + (-ka1);
  float *qyzd = msm->lyzd_f + (-ka1*nj1 + -ja1);
  float qsum;

  const float *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 float *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

IndexOffset

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_sprec.c.

Referenced by prolongation(), and restriction().

Nstencil

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_sprec.c.

Referenced by prolongation(), and restriction().

PhiStencil

const float 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_sprec.c.

Referenced by prolongation(), and restriction().

PhiStencilFactored

const float 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_sprec.c.

Referenced by prolongation_factored(), and restriction_factored().

PolyDegree

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_sprec.c.

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