msm_setup.c File Reference

#include "msm_defn.h"

Go to the source code of this file.

Classes
struct	InterpParams_t

Typedefs
typedef struct InterpParams_t	InterpParams

Functions
void	NL_msm_cleanup (NL_Msm *pm)
static int	setup_bins_1away (NL_Msm *pm)
static int	setup_bins_k_away (NL_Msm *pm)
static int	setup_cell_vectors (NL_Msm *pm)
static int	setup_grids (NL_Msm *pm)
static int	setup_hgrid_1d (NL_Msm *pm, double len, double *hh, int *nn, int *aindex, int *bindex, int isperiodic)
static int	print_status (NL_Msm *msm)
int	NL_msm_setup (NL_Msm *pm, double cutoff, double cellvec1[3], double cellvec2[3], double cellvec3[3], double cellcenter[3], int msmflags)

Variables
static InterpParams	INTERP_PARAMS []

Typedef Documentation

InterpParams

typedef struct InterpParams_t InterpParams

Function Documentation

NL_msm_cleanup()

void NL_msm_cleanup

(

NL_Msm *

pm

)

Definition at line 7 of file msm_setup.c.

References NL_Msm_t::eh, NL_Msm_t::eh_f, NL_Msm_t::gc, NL_Msm_t::gc_f, GRID_DONE, NL_Msm_t::lyzd, NL_Msm_t::lyzd_f, NL_Msm_t::lzd, NL_Msm_t::lzd_f, NL_Msm_t::maxlevels, NL_Msm_t::msmflags, NL_MSM_COMPUTE_CUDA_GRID_CUTOFF, NL_MSM_COMPUTE_SPREC, NL_Msm_t::qh, and NL_Msm_t::qh_f.

Referenced by NL_msm_destroy().

                                {
  int i;
#ifdef NL_USE_CUDA
  if (pm->msmflags & NL_MSM_COMPUTE_CUDA_GRID_CUTOFF) {
    NL_msm_cuda_cleanup_gridcutoff(pm);
  }
#endif /* NL_USE_CUDA */
  if (pm->msmflags & NL_MSM_COMPUTE_SPREC) {
    for (i = 0;  i < pm->maxlevels;  i++) {
      GRID_DONE( &(pm->qh_f[i]) );
      GRID_DONE( &(pm->eh_f[i]) );
      GRID_DONE( &(pm->gc_f[i]) );
    }
  }
  else {
    for (i = 0;  i < pm->maxlevels;  i++) {
      GRID_DONE( &(pm->qh[i]) );
      GRID_DONE( &(pm->eh[i]) );
      GRID_DONE( &(pm->gc[i]) );
    }
  }
  free(pm->lzd);
  free(pm->lyzd);
  free(pm->lzd_f);
  free(pm->lyzd_f);
  free(pm->qh);
  free(pm->eh);
  free(pm->gc);
  free(pm->qh_f);
  free(pm->eh_f);
  free(pm->gc_f);
}

NL_msm_setup()

int NL_msm_setup	(	NL_Msm *	msm,
		double	cutoff,
		double	cellvec1[3],
		double	cellvec2[3],
		double	cellvec3[3],
		double	cellcenter[3],
		int	msmflags
	)

Setup the MSM solver for the molecular system. These parameters come directly from the simulation. Here the cutoff provides some control over the accuracy. The region of space containing the atoms is defined as a parallelepiped, and the "isperiodic" flag determines whether or not each individual cell dimension has periodic boundary conditions.
The defined cell is expected to contain the atoms regardless of the choice of boundary conditions. The cell volume must contain the cutoff-sphere.
XXX For now, the MSM solver permits only cell vectors that align with the x-, y-, and z-axis, respectively.

Parameters

pm	the MSM solver object
cutoff	cutoff distance for short-range part
cellvec1	cell basis vector 1
cellvec2	cell basis vector 2
cellvec3	cell basis vector 3
cellcenter	center of cell
msmflags	flags for periodicity and calculation

Definition at line 62 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::a_f, NL_Msm_t::cellcenter, NL_Msm_t::cellcenter_f, NL_Msm_t::cellvec1, NL_Msm_t::cellvec1_f, NL_Msm_t::cellvec2, NL_Msm_t::cellvec2_f, NL_Msm_t::cellvec3, NL_Msm_t::cellvec3_f, NL_Msm_t::gx, NL_Msm_t::gx_f, NL_Msm_t::gy, NL_Msm_t::gy_f, NL_Msm_t::gz, NL_Msm_t::gz_f, NL_Msm_t::hx, NL_Msm_t::hx_f, NL_Msm_t::hy, NL_Msm_t::hy_f, NL_Msm_t::hz, NL_Msm_t::hz_f, NL_Msm_t::msmflags, NL_MSM_ALL_FLAGS, NL_MSM_COMPUTE_1AWAY, NL_MSM_COMPUTE_ALL, NL_MSM_COMPUTE_CUDA_FALL_BACK, NL_MSM_COMPUTE_CUDA_GRID_CUTOFF, NL_MSM_COMPUTE_LONG_RANGE, NL_MSM_COMPUTE_SHORT_RANGE, NL_MSM_COMPUTE_SPREC, NL_MSM_ERROR_PARAM, NL_MSM_ERROR_SUPPORT, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_REPORT_TIMINGS, NL_MSM_SUCCESS, print_status(), NL_Msm_t::recipvec1, NL_Msm_t::recipvec1_f, NL_Msm_t::recipvec2, NL_Msm_t::recipvec2_f, NL_Msm_t::recipvec3, NL_Msm_t::recipvec3_f, NL_Msm_t::report_timings, setup_bins_1away(), setup_bins_k_away(), setup_cell_vectors(), and setup_grids().

Referenced by ComputeMsmSerialMgr::recvCoord().

      {

  int rc = 0;

  /* check sanity of msmflags */
  if ((~NL_MSM_ALL_FLAGS & msmflags) != 0 ||
      (NL_MSM_COMPUTE_ALL & msmflags) == 0) {
    return NL_MSM_ERROR_PARAM;
  }

  /* report timings? */
  pm->report_timings = ((msmflags & NL_MSM_REPORT_TIMINGS) != 0);

  /* for now, permit only orthogonal cell aligned with coordinate axes */
  if (cellvec1[1] != 0 || cellvec1[2] != 0 ||
      cellvec2[0] != 0 || cellvec2[2] != 0 ||
      cellvec3[0] != 0 || cellvec3[1] != 0 ||
      cellvec1[0] <= 0 || cellvec2[1] <= 0 || cellvec3[2] <= 0) {
    return NL_MSM_ERROR_SUPPORT;
  }

  /* check sanity of cutoff wrt expected cell;
   * XXX cell widths must be at least the cutoff distance length */
  if (cutoff <= 0 ||
      (cutoff > cellvec1[0] && (msmflags & NL_MSM_PERIODIC_VEC1)) ||
      (cutoff > cellvec2[1] && (msmflags & NL_MSM_PERIODIC_VEC2)) ||
      (cutoff > cellvec3[2] && (msmflags & NL_MSM_PERIODIC_VEC3))) {
    return NL_MSM_ERROR_PARAM;
  }

  pm->msmflags = msmflags;
  pm->cellvec1[0] = cellvec1[0];
  pm->cellvec1[1] = cellvec1[1];
  pm->cellvec1[2] = cellvec1[2];
  pm->cellvec2[0] = cellvec2[0];
  pm->cellvec2[1] = cellvec2[1];
  pm->cellvec2[2] = cellvec2[2];
  pm->cellvec3[0] = cellvec3[0];
  pm->cellvec3[1] = cellvec3[1];
  pm->cellvec3[2] = cellvec3[2];
  pm->cellcenter[0] = cellcenter[0];
  pm->cellcenter[1] = cellcenter[1];
  pm->cellcenter[2] = cellcenter[2];

  pm->a = cutoff;

  rc = setup_cell_vectors(pm);
  if (rc) return rc;

  if (msmflags & NL_MSM_COMPUTE_SHORT_RANGE) {
    /* set up bins for short-range part */
    if (msmflags & NL_MSM_COMPUTE_1AWAY) {
      rc = setup_bins_1away(pm);
    }
    else {
      rc = setup_bins_k_away(pm);
    }
    if (rc) return rc;
  }

  if (msmflags & NL_MSM_COMPUTE_LONG_RANGE) {
    /* set up grid hierarchy for long-range part */
    rc = setup_grids(pm);
    if (rc) return rc;
  }

  if (msmflags & NL_MSM_COMPUTE_SPREC) {
    /* fill out single precision data if needed */
    pm->cellvec1_f[0] = (float) pm->cellvec1[0];
    pm->cellvec1_f[1] = (float) pm->cellvec1[1];
    pm->cellvec1_f[2] = (float) pm->cellvec1[2];
    pm->cellvec2_f[0] = (float) pm->cellvec2[0];
    pm->cellvec2_f[1] = (float) pm->cellvec2[1];
    pm->cellvec2_f[2] = (float) pm->cellvec2[2];
    pm->cellvec3_f[0] = (float) pm->cellvec3[0];
    pm->cellvec3_f[1] = (float) pm->cellvec3[1];
    pm->cellvec3_f[2] = (float) pm->cellvec3[2];
    pm->cellcenter_f[0] = (float) pm->cellcenter[0];
    pm->cellcenter_f[1] = (float) pm->cellcenter[1];
    pm->cellcenter_f[2] = (float) pm->cellcenter[2];
    pm->recipvec1_f[0] = (float) pm->recipvec1[0];
    pm->recipvec1_f[1] = (float) pm->recipvec1[1];
    pm->recipvec1_f[2] = (float) pm->recipvec1[2];
    pm->recipvec2_f[0] = (float) pm->recipvec2[0];
    pm->recipvec2_f[1] = (float) pm->recipvec2[1];
    pm->recipvec2_f[2] = (float) pm->recipvec2[2];
    pm->recipvec3_f[0] = (float) pm->recipvec3[0];
    pm->recipvec3_f[1] = (float) pm->recipvec3[1];
    pm->recipvec3_f[2] = (float) pm->recipvec3[2];
    pm->hx_f = pm->hx;
    pm->hy_f = pm->hy;
    pm->hz_f = pm->hz;
    pm->a_f = pm->a;
    pm->gx_f = pm->gx;
    pm->gy_f = pm->gy;
    pm->gz_f = pm->gz;
  }

  print_status(pm);

#ifdef NL_USE_CUDA
  if (msmflags & NL_MSM_COMPUTE_CUDA_GRID_CUTOFF) {
    rc = NL_msm_cuda_setup_gridcutoff(pm);
    if (rc == NL_MSM_SUCCESS) {
      printf("Using CUDA for grid cutoff computation\n");
    }
    else {
      printf("Unable to set up CUDA for grid cutoff computation\n");
      if (msmflags & NL_MSM_COMPUTE_CUDA_FALL_BACK) {
        NL_msm_cuda_cleanup_gridcutoff(pm);
        printf("Falling back on CPU\n");
        pm->msmflags &= ~NL_MSM_COMPUTE_CUDA_GRID_CUTOFF;
      }
      else return rc;
    }
  }
#else
  if (msmflags & NL_MSM_COMPUTE_CUDA_GRID_CUTOFF) {
    if (msmflags & NL_MSM_COMPUTE_CUDA_FALL_BACK) {
      printf("Falling back on CPU\n");
      pm->msmflags &= ~NL_MSM_COMPUTE_CUDA_GRID_CUTOFF;
    }
    else return NL_MSM_ERROR_SUPPORT;
  }
#endif /* NL_USE_CUDA */

  return NL_MSM_SUCCESS;
}

print_status()

int print_status ( NL_Msm *  msm )

static

Definition at line 218 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::approx, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::gridspacing, NL_Msm_t::gx, NL_Msm_t::gy, NL_Msm_t::gz, NL_Msm_t::hx, NL_Msm_t::hy, NL_Msm_t::hz, NL_Msm_t::msmflags, NL_msm_approx_name(), NL_MSM_COMPUTE_SPREC, NL_MSM_PERIODIC_ALL, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_msm_split_name(), NL_MSM_SUCCESS, NL_Msm_t::nlevels, NL_Msm_t::nx, NL_Msm_t::ny, NL_Msm_t::nz, NL_Msm_t::qh, NL_Msm_t::qh_f, and NL_Msm_t::split.

Referenced by NL_msm_setup().

                             {
  int k;
  int ispx = (pm->msmflags & NL_MSM_PERIODIC_VEC1);
  int ispy = (pm->msmflags & NL_MSM_PERIODIC_VEC2);
  int ispz = (pm->msmflags & NL_MSM_PERIODIC_VEC3);
  int ispany = (pm->msmflags & NL_MSM_PERIODIC_ALL);
  int ispall = (ispany == NL_MSM_PERIODIC_ALL);

  const double xlen = pm->cellvec1[0];  /* XXX */
  const double ylen = pm->cellvec2[1];
  const double zlen = pm->cellvec3[2];

  printf("#MSM SETUP\n");
  printf("#  cell lengths=  %g  %g  %g\n", xlen, ylen, zlen);
  printf("#  grid origin=  %g  %g  %g\n", pm->gx, pm->gy, pm->gz);
  if (ispall) {
    printf("#  periodic boundaries\n");
  }
  else if (!ispany) {
    printf("#  non-periodic boundaries\n");
  }
  else {
    printf("#  periodic boundaries in:%s%s%s\n",
        (ispx ? "  x" : ""), (ispy ? "  y" : ""), (ispz ? "  z" : ""));
  }
  printf("#  cutoff= %g\n", pm->a);
  printf("#  grid spacing= %g\n", pm->gridspacing);
  printf("#  hx, hy, hz= %g, %g, %g\n", pm->hx, pm->hy, pm->hz);
  printf("#  h-grid size= %d x %d x %d\n", pm->nx, pm->ny, pm->nz);
  printf("#  number of levels= %d\n", pm->nlevels);
  printf("#  approximation= %s\n", NL_msm_approx_name(pm->approx));
  printf("#  splitting= %s\n", NL_msm_split_name(pm->split));
  printf("#  grid hierarchy:\n");
  if (pm->msmflags & NL_MSM_COMPUTE_SPREC) {
    for (k = 0;  k < pm->nlevels;  k++) {
      NL_Msmgrid_float *qh = &(pm->qh_f[k]);
      int ia = qh->i0;
      int ib = ia + qh->ni - 1;
      int ja = qh->j0;
      int jb = ja + qh->nj - 1;
      int ka = qh->k0;
      int kb = ka + qh->nk - 1;
      printf("#  level= %d:  [%d..%d] x [%d..%d] x [%d..%d]\n",
          k, ia, ib, ja, jb, ka, kb);
    }
  }
  else {
    for (k = 0;  k < pm->nlevels;  k++) {
      NL_Msmgrid_double *qh = &(pm->qh[k]);
      int ia = qh->i0;
      int ib = ia + qh->ni - 1;
      int ja = qh->j0;
      int jb = ja + qh->nj - 1;
      int ka = qh->k0;
      int kb = ka + qh->nk - 1;
      printf("#  level= %d:  [%d..%d] x [%d..%d] x [%d..%d]\n",
          k, ia, ib, ja, jb, ka, kb);
    }
  }
  return NL_MSM_SUCCESS;
}

setup_bins_1away()

int setup_bins_1away ( NL_Msm *  pm )

static

Definition at line 318 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::bin, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::maxbins, NL_Msm_t::msmflags, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, NL_MSM_ERROR_MALLOC, NL_MSM_ERROR_PARAM, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::numbins, NL_Msm_t::recipvec1, NL_Msm_t::recipvec2, NL_Msm_t::recipvec3, X, Y, and Z.

Referenced by NL_msm_setup().

                                 {
  const double *u = pm->cellvec1;
  const double *v = pm->cellvec2;
  const double *w = pm->cellvec3;
  double *ru = pm->recipvec1;
  double *rv = pm->recipvec2;
  double *rw = pm->recipvec3;
  double p[3];
  double pulen, pvlen, pwlen, s;
  double cutoff = pm->a;
  int nbx, nby, nbz, numbins;
  int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
  int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
  int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);

  /* calculate the number of atom bins in each basis vector dimension,
   * such that each bin (a parallelepiped) inscribes the cutoff-cube;
   * for periodic boundaries, we have to choose equal sized bins with
   * length >= cutoff that tile the cell; for non-periodic boundaries,
   * we can have bins of length = cutoff */

  /* find the largest orthogonal box inscribed within parallelepiped cell
   * by taking orthogonal projections onto cross products of basis vectors */

  p[X] = v[Y]*w[Z] - v[Z]*w[Y];  /* p = v CROSS w */
  p[Y] = v[Z]*w[X] - v[X]*w[Z];
  p[Z] = v[X]*w[Y] - v[Y]*w[X];
  /* s = (u DOT p) / (p DOT p) */
  s = (u[X]*p[X] + u[Y]*p[Y] + u[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  p[X] *= s;  /* p is orthogonal projection of u onto v CROSS w */
  p[Y] *= s;
  p[Z] *= s;
  pulen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);

  p[X] = w[Y]*u[Z] - w[Z]*u[Y];  /* p = w CROSS u */
  p[Y] = w[Z]*u[X] - w[X]*u[Z];
  p[Z] = w[X]*u[Y] - w[Y]*u[X];
  /* s = (v DOT p) / (p DOT p) */
  s = (v[X]*p[X] + v[Y]*p[Y] + v[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  p[X] *= s;  /* p is orthogonal projection of v onto w CROSS u */
  p[Y] *= s;
  p[Z] *= s;
  pvlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);

  p[X] = u[Y]*v[Z] - u[Z]*v[Y];  /* p = u CROSS v */
  p[Y] = u[Z]*v[X] - u[X]*v[Z];
  p[Z] = u[X]*v[Y] - u[Y]*v[X];
  /* s = (w DOT p) / (p DOT p) */
  s = (w[X]*p[X] + w[Y]*p[Y] + w[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  p[X] *= s;  /* p is orthogonal projection of w onto u CROSS v */
  p[Y] *= s;
  p[Z] *= s;
  pwlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);

  /* find the number of cells along each dimension by dividing the
   * orthogonal projection vector lengths by the cutoff distance */
  nbx = (ispx ? (int) floor(pulen / cutoff) : (int) ceil(pulen / cutoff));
  nby = (ispy ? (int) floor(pvlen / cutoff) : (int) ceil(pvlen / cutoff));
  nbz = (ispz ? (int) floor(pwlen / cutoff) : (int) ceil(pwlen / cutoff));
  if (nbx == 0 || nby == 0 || nbz == 0) {
    return NL_MSM_ERROR_PARAM;  /* cutoff is too large for cell basis */
  }
  numbins = nbx * nby * nbz;
  /* we don't know the number of atoms until compute */

  /* allocate one atom index for each bin */
  if (pm->maxbins < numbins) {
    void *vptr = malloc(numbins * sizeof(int));
    if (vptr == NULL) return NL_MSM_ERROR_MALLOC;
    pm->bin = (int *) vptr;
    pm->maxbins = numbins;
  }
  pm->nbx = nbx;
  pm->nby = nby;
  pm->nbz = nbz;
  pm->numbins = numbins;

  /* must allocate next index array when we know number of atoms */

  /* rescale recip space vectors for non-periodic expansion */
  if (ispx == 0) {
    s = pulen / (nbx * cutoff);
    ru[X] *= s;
    ru[Y] *= s;
    ru[Z] *= s;
  }
  if (ispy == 0) {
    s = pvlen / (nby * cutoff);
    rv[X] *= s;
    rv[Y] *= s;
    rv[Z] *= s;
  }
  if (ispz == 0) {
    s = pwlen / (nbz * cutoff);
    rw[X] *= s;
    rw[Y] *= s;
    rw[Z] *= s;
  }
  return NL_MSM_SUCCESS;
}

setup_bins_k_away()

int setup_bins_k_away ( NL_Msm *  pm )

static

Definition at line 420 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::bin, NL_Msm_t::binfill, NL_Msm_t::bu, NL_Msm_t::bv, NL_Msm_t::bw, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::density, NL_Msm_t::maxbins, NL_Msm_t::msmflags, NL_Msm_t::nbinslots, NL_Msm_t::nbrhd, NL_Msm_t::nbrhdlen, NL_Msm_t::nbrhdmax, NL_Msm_t::nbx, NL_Msm_t::nby, NL_Msm_t::nbz, NL_MSM_ERROR_MALLOC, NL_MSM_ERROR_PARAM, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::nradius, NL_Msm_t::numbins, NL_Msm_t::recipvec1, NL_Msm_t::recipvec2, NL_Msm_t::recipvec3, X, Y, and Z.

Referenced by NL_msm_setup().

                                  {
  double *u = pm->cellvec1;
  double *v = pm->cellvec2;
  double *w = pm->cellvec3;
  double *ru = pm->recipvec1;
  double *rv = pm->recipvec2;
  double *rw = pm->recipvec3;
  double *bu = pm->bu;
  double *bv = pm->bv;
  double *bw = pm->bw;
  double p[3];
  double pulen, pvlen, pwlen, s;
  double cutoff = pm->a;
  int nbx, nby, nbz, numbins;
  int ispx = ((pm->msmflags & NL_MSM_PERIODIC_VEC1) != 0);
  int ispy = ((pm->msmflags & NL_MSM_PERIODIC_VEC2) != 0);
  int ispz = ((pm->msmflags & NL_MSM_PERIODIC_VEC3) != 0);
  int nbrhdmax, nradius, ndiameter, index, i, j, k;
  double volume;
  double binvol, abincnt, binlen, c;
  double min2;

  /* find the largest orthogonal box inscribed within parallelepiped cell
   * by taking orthogonal projections onto cross products of basis vectors */
  p[X] = v[Y]*w[Z] - v[Z]*w[Y];  /* p = v CROSS w */
  p[Y] = v[Z]*w[X] - v[X]*w[Z];
  p[Z] = v[X]*w[Y] - v[Y]*w[X];
  /* s = (u DOT p) / (p DOT p) */
  s = (u[X]*p[X] + u[Y]*p[Y] + u[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  p[X] *= s;  /* p is orthogonal projection of u onto v CROSS w */
  p[Y] *= s;
  p[Z] *= s;
  pulen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);

  p[X] = w[Y]*u[Z] - w[Z]*u[Y];  /* p = w CROSS u */
  p[Y] = w[Z]*u[X] - w[X]*u[Z];
  p[Z] = w[X]*u[Y] - w[Y]*u[X];
  /* s = (v DOT p) / (p DOT p) */
  s = (v[X]*p[X] + v[Y]*p[Y] + v[Z]*p[Z]) / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  p[X] *= s;  /* p is orthogonal projection of v onto w CROSS u */
  p[Y] *= s;
  p[Z] *= s;
  pvlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);

  p[X] = u[Y]*v[Z] - u[Z]*v[Y];  /* p = u CROSS v */
  p[Y] = u[Z]*v[X] - u[X]*v[Z];
  p[Z] = u[X]*v[Y] - u[Y]*v[X];
  volume = w[X]*p[X] + w[Y]*p[Y] + w[Z]*p[Z];
  /* s = (w DOT p) / (p DOT p) */
  s = volume / (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  p[X] *= s;  /* p is orthogonal projection of w onto u CROSS v */
  p[Y] *= s;
  p[Z] *= s;
  pwlen = sqrt(p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);

  if ((ispx && cutoff > pulen) ||
      (ispy && cutoff > pvlen) ||
      (ispz && cutoff > pwlen)) {
    printf("Cutoff %.3g too big for cell along basis%s%s%s\n",
        cutoff, (cutoff > pulen ? " x" : ""), (cutoff > pvlen ?  " y" : ""),
        (cutoff > pwlen ?  " z" : ""));
    return NL_MSM_ERROR_PARAM;
  }

  /* calculate the ideal bin volume based on a fixed bin size (nbinslots),
   * the particle density (density), and a desired fill ratio (binfill) */
  binvol = pm->binfill * pm->nbinslots / pm->density;
  abincnt = volume / binvol;
  binlen = pow(binvol, 1./3);  /* ideal side length - use for nonperiodic */

  /* factor "abincnt" into 3 parts, each part proportional to the
   * lengths of the orthogonal projection vectors calculated above */
  c = pow(abincnt / (pulen*pvlen*pwlen), 1./3);  /* proportionality const */
  nbx = (int) ceil(c * pulen);
  nby = (int) ceil(c * pvlen);
  nbz = (int) ceil(c * pwlen);
  numbins = nbx * nby * nbz;

  printf("nbx=%d  nby=%d  nbz=%d  numbins=%d\n", nbx, nby, nbz, numbins);

  /* allocate one atom index for each bin */
  if (pm->maxbins < numbins) {
    void *vptr = malloc(numbins * sizeof(int));
    if (vptr == NULL) return NL_MSM_ERROR_MALLOC;
    pm->bin = (int *) vptr;
    pm->maxbins = numbins;
  }
  pm->nbx = nbx;
  pm->nby = nby;
  pm->nbz = nbz;
  pm->numbins = numbins;

  /* rescale basis and recip space vectors for non-periodic expansion */
  if ( ! ispx) {
    s = pulen / (nbx * binlen);
    ru[X] *= s;
    ru[Y] *= s;
    ru[Z] *= s;
    s = 1./s;
    u[X] *= s;
    u[Y] *= s;
    u[Z] *= s;
  }
  if ( ! ispy) {
    s = pvlen / (nby * binlen);
    rv[X] *= s;
    rv[Y] *= s;
    rv[Z] *= s;
    s = 1./s;
    v[X] *= s;
    v[Y] *= s;
    v[Z] *= s;
  }
  if ( ! ispz) {
    s = pwlen / (nbz * binlen);
    rw[X] *= s;
    rw[Y] *= s;
    rw[Z] *= s;
    s = 1./s;
    w[X] *= s;
    w[Y] *= s;
    w[Z] *= s;
  }

  /* scale basis vectors by number of bins to get bin basis vectors */
  s = 1./nbx;
  bu[X] = s * u[X];
  bu[Y] = s * u[Y];
  bu[Z] = s * u[Z];
  s = 1./nby;
  bv[X] = s * v[X];
  bv[Y] = s * v[Y];
  bv[Z] = s * v[Z];
  s = 1./nbz;
  bw[X] = s * w[X];
  bw[Y] = s * w[Y];
  bw[Z] = s * w[Z];

  /* determine the neighborhood of bins */

  /* first find minimum width of cell */
  min2 = (bu[X]*bu[X] + bu[Y]*bu[Y] + bu[Z]*bu[Z]);
  s = (bv[X]*bv[X] + bv[Y]*bv[Y] + bv[Z]*bv[Z]);
  if (min2 > s) min2 = s;
  s = (bw[X]*bw[X] + bw[Y]*bw[Y] + bw[Z]*bw[Z]);
  if (min2 > s) min2 = s;

  /* also find the minimum of the four major diagonals */
  p[X] = bu[X] + bv[X] + bw[X];
  p[Y] = bu[Y] + bv[Y] + bw[Y];
  p[Z] = bu[Z] + bv[Z] + bw[Z];
  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  if (min2 > s) min2 = s;
  p[X] = -bu[X] + bv[X] + bw[X];
  p[Y] = -bu[Y] + bv[Y] + bw[Y];
  p[Z] = -bu[Z] + bv[Z] + bw[Z];
  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  if (min2 > s) min2 = s;
  p[X] = bu[X] - bv[X] + bw[X];
  p[Y] = bu[Y] - bv[Y] + bw[Y];
  p[Z] = bu[Z] - bv[Z] + bw[Z];
  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  if (min2 > s) min2 = s;
  p[X] = bu[X] + bv[X] - bw[X];
  p[Y] = bu[Y] + bv[Y] - bw[Y];
  p[Z] = bu[Z] + bv[Z] - bw[Z];
  s = (p[X]*p[X] + p[Y]*p[Y] + p[Z]*p[Z]);
  if (min2 > s) min2 = s;

  /* the neighborhood is contained in the cube of nradius bins,
   * store only upper half of bin neighborhood */
  nradius = (int) ceil(sqrt(cutoff*cutoff / min2));
  ndiameter = 2 * nradius + 1;
  nbrhdmax = (ndiameter * ndiameter * ndiameter) / 2 + 1;

#if 0
  printf("Neighborhood radius = %d\n", nradius);
  printf("min2 = %g\n", min2);
#endif

  /* allocate space for the entire cube */
  if (pm->nbrhd) free(pm->nbrhd);
  pm->nbrhd = (int *) calloc(3*nbrhdmax, sizeof(int));
  if (pm->nbrhd == NULL) return NL_MSM_ERROR_MALLOC;
  pm->nbrhdmax = nbrhdmax;
  pm->nradius = nradius;  /* save neighborhood radius for diagnostic purposes */

  /* trim the neighborhood */
  index = 0;
  for (k = -nradius;  k <= nradius;  k++) {
    for (j = -nradius;  j <= nradius;  j++) {
      for (i = -nradius;  i <= nradius;  i++) {
        /* XXX should we remove (0,0,0) from bin neighborhood? */
        int ip, jp, kp, iq, jq, kq;
        int is_within = 0;
        int binindex = (k * ndiameter + j) * ndiameter + i;
        if (binindex < 0) continue;
        for (kp = 0;  kp <= 1;  kp++) {
          for (jp = 0;  jp <= 1;  jp++) {
            for (ip = 0;  ip <= 1;  ip++) {
              p[X] = (i+ip)*bu[X] + (j+jp)*bv[X] + (k+kp)*bw[X];
              p[Y] = (i+ip)*bu[Y] + (j+jp)*bv[Y] + (k+kp)*bw[Y];
              p[Z] = (i+ip)*bu[Z] + (j+jp)*bv[Z] + (k+kp)*bw[Z];
              for (kq = 0;  kq <= 1;  kq++) {
                for (jq = 0;  jq <= 1;  jq++) {
                  for (iq = 0;  iq <= 1;  iq++) {
                    double q[3];
                    double dx, dy, dz, r2;
                    q[X] = iq*bu[X] + jq*bv[X] + kq*bw[X];
                    q[Y] = iq*bu[Y] + jq*bv[Y] + kq*bw[Y];
                    q[Z] = iq*bu[Z] + jq*bv[Z] + kq*bw[Z];
                    dx = p[X] - q[X];
                    dy = p[Y] - q[Y];
                    dz = p[Z] - q[Z];
                    r2 = dx*dx + dy*dy + dz*dz;
                    is_within |= (r2 < cutoff*cutoff);
                  }
                }
              } /* end q loop */

            }
          }
        } /* end p loop */
        if (is_within) {
          pm->nbrhd[index++] = i;
          pm->nbrhd[index++] = j;
          pm->nbrhd[index++] = k;
        }

      }
    }
  }
  pm->nbrhdlen = index / 3;

  return NL_MSM_SUCCESS;
}

setup_cell_vectors()

int setup_cell_vectors ( NL_Msm *  pm )

static

Definition at line 281 of file msm_setup.c.

References NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_MSM_SUCCESS, NL_Msm_t::recipvec1, NL_Msm_t::recipvec2, NL_Msm_t::recipvec3, X, Y, and Z.

Referenced by NL_msm_setup().

                                   {
  double *u = pm->cellvec1;
  double *v = pm->cellvec2;
  double *w = pm->cellvec3;
  double *bu = pm->recipvec1;
  double *bv = pm->recipvec2;
  double *bw = pm->recipvec3;
  double c[3], s;

  c[X] = v[Y]*w[Z] - v[Z]*w[Y];  /* v CROSS w */
  c[Y] = v[Z]*w[X] - v[X]*w[Z];
  c[Z] = v[X]*w[Y] - v[Y]*w[X];
  s = 1 / (u[X]*c[X] + u[Y]*c[Y] + u[Z]*c[Z]);  /* 1 / (c DOT u) */
  bu[X] = s*c[X];
  bu[Y] = s*c[Y];
  bu[Z] = s*c[Z];

  c[X] = w[Y]*u[Z] - w[Z]*u[Y];  /* w CROSS u */
  c[Y] = w[Z]*u[X] - w[X]*u[Z];
  c[Z] = w[X]*u[Y] - w[Y]*u[X];
  s = 1 / (v[X]*c[X] + v[Y]*c[Y] + v[Z]*c[Z]);  /* 1 / (c DOT v) */
  bv[X] = s*c[X];
  bv[Y] = s*c[Y];
  bv[Z] = s*c[Z];

  c[X] = u[Y]*v[Z] - u[Z]*v[Y];  /* u CROSS v */
  c[Y] = u[Z]*v[X] - u[X]*v[Z];
  c[Z] = u[X]*v[Y] - u[Y]*v[X];
  s = 1 / (w[X]*c[X] + w[Y]*c[Y] + w[Z]*c[Z]);  /* 1 / (c DOT w) */
  bw[X] = s*c[X];
  bw[Y] = s*c[Y];
  bw[Z] = s*c[Z];

  return NL_MSM_SUCCESS;
}

setup_grids()

int setup_grids ( NL_Msm *  pm )

static

Definition at line 707 of file msm_setup.c.

References NL_Msm_t::a, NL_Msm_t::approx, ASSERT, NL_Msm_t::cellcenter, NL_Msm_t::cellvec1, NL_Msm_t::cellvec2, NL_Msm_t::cellvec3, NL_Msm_t::eh, NL_Msm_t::eh_f, NL_Msm_t::gc, NL_Msm_t::gc_f, GRID_INDEX, GRID_INDEX_CHECK, GRID_INIT, GRID_RESIZE, NL_Msm_t::gx, NL_Msm_t::gy, NL_Msm_t::gz, NL_Msm_t::gzero, NL_Msm_t::gzero_f, NL_Msm_t::hx, NL_Msm_t::hy, NL_Msm_t::hz, INTERP_PARAMS, NL_Msm_t::lyzd, NL_Msm_t::lyzd_f, NL_Msm_t::lzd, NL_Msm_t::lzd_f, NL_Msm_t::max_lyzd, NL_Msm_t::max_lzd, NL_Msm_t::maxlevels, NL_Msm_t::msmflags, NL_MSM_COMPUTE_SPREC, NL_MSM_ERROR_MALLOC, NL_MSM_PERIODIC_ALL, NL_MSM_PERIODIC_VEC1, NL_MSM_PERIODIC_VEC2, NL_MSM_PERIODIC_VEC3, NL_MSM_SUCCESS, NL_Msm_t::nlevels, InterpParams_t::nu, NL_Msm_t::nx, NL_Msm_t::ny, NL_Msm_t::nz, InterpParams_t::omega, NL_Msm_t::qh, NL_Msm_t::qh_f, setup_hgrid_1d(), split(), NL_Msm_t::split, and SPOLY.

Referenced by NL_msm_setup().

                            {
  const int nu = INTERP_PARAMS[pm->approx].nu;
  const int omega = INTERP_PARAMS[pm->approx].omega;
  const int split = pm->split;
  const int ispx = (pm->msmflags & NL_MSM_PERIODIC_VEC1);
  const int ispy = (pm->msmflags & NL_MSM_PERIODIC_VEC2);
  const int ispz = (pm->msmflags & NL_MSM_PERIODIC_VEC3);
  const int ispany = (pm->msmflags & NL_MSM_PERIODIC_ALL);

  const int issprec = (pm->msmflags & NL_MSM_COMPUTE_SPREC);

  const double xlen = pm->cellvec1[0];  /* XXX */
  const double ylen = pm->cellvec2[1];
  const double zlen = pm->cellvec3[2];

  const double a = pm->a;
  double hx, hy, hz;
  double scaling;
  double d;  /* temporary for SPOLY derivative */

  NL_Msmgrid_double *p = NULL;
  NL_Msmgrid_float *p_f = NULL;
  int ia, ib, ja, jb, ka, kb, ni, nj, nk;
  int nx, ny, nz;  /* counts the grid points that span just the domain */

  int i, j, k, n;
  int index;
  int level, toplevel, nlevels, maxlevels;
  int lastnelems = 1;
  int isclamped = 0;
  int done, alldone;

  int rc = 0;

  rc = setup_hgrid_1d(pm, xlen, &hx, &nx, &ia, &ib, ispx);
  if (rc) return rc;

  rc = setup_hgrid_1d(pm, ylen, &hy, &ny, &ja, &jb, ispy);
  if (rc) return rc;

  rc = setup_hgrid_1d(pm, zlen, &hz, &nz, &ka, &kb, ispz);
  if (rc) return rc;

  pm->hx = hx;
  pm->hy = hy;
  pm->hz = hz;

  /* XXX set coordinate for h-grid (0,0,0) point */
  pm->gx = pm->cellcenter[0] - ((nx >> 1) * hx);
  pm->gy = pm->cellcenter[1] - ((ny >> 1) * hy);
  pm->gz = pm->cellcenter[2] - ((nz >> 1) * hz);

  pm->nx = nx;
  pm->ny = ny;
  pm->nz = nz;

  ni = ib - ia + 1;
  nj = jb - ja + 1;
  nk = kb - ka + 1;

  /* allocate temp buffer space for factored grid transfer */
  n = (nk > omega ? nk : omega);  /* row along z-dimension */
  if (pm->max_lzd < n) {
    if (issprec) {
      float *t;
      t = (float *) realloc(pm->lzd, n * sizeof(float));
      if (NULL == t) return NL_MSM_ERROR_MALLOC;
      pm->lzd_f = t;
    }
    else {
      double *t;
      t = (double *) realloc(pm->lzd, n * sizeof(double));
      if (NULL == t) return NL_MSM_ERROR_MALLOC;
      pm->lzd = t;
    }
    pm->max_lzd = n;
  }
  n *= (nj > omega ? nj : omega);  /* plane along yz-dimensions */
  if (pm->max_lyzd < n) {
    if (issprec) {
      float *t;
      t = (float *) realloc(pm->lyzd, n * sizeof(float));
      if (NULL == t) return NL_MSM_ERROR_MALLOC;
      pm->lyzd_f = t;
    }
    else {
      double *t;
      t = (double *) realloc(pm->lyzd, n * sizeof(double));
      if (NULL == t) return NL_MSM_ERROR_MALLOC;
      pm->lyzd = t;
    }
    pm->max_lyzd = n;
  }

  nlevels = pm->nlevels;
  if (nlevels <= 0) {
    /* automatically set number of levels */
    n = ni;
    if (n < nj) n = nj;
    if (n < nk) n = nk;
    for (maxlevels = 1;  n > 0;  n >>= 1)  maxlevels++;
    nlevels = maxlevels;
    if (ispany == 0) {  /* no periodicity */
      int omega3 = omega * omega * omega;
      int nhalf = (int) sqrtf(ni*nj*nk);  /* scale down for performance? */
      lastnelems = (nhalf > omega3 ? nhalf : omega3);
      isclamped = 1;
    }
  }
  else {
    /* user-defined number of levels */
    maxlevels = nlevels;
  }

  /* allocate any additional levels that may be needed */
  if (pm->maxlevels < maxlevels) {
    void *vqh, *veh, *vgc;
    if (issprec) {
      vqh = realloc(pm->qh, maxlevels * sizeof(NL_Msmgrid_float));
      if (NULL == vqh) return NL_MSM_ERROR_MALLOC;
      veh = realloc(pm->eh, maxlevels * sizeof(NL_Msmgrid_float));
      if (NULL == veh) return NL_MSM_ERROR_MALLOC;
      vgc = realloc(pm->gc, maxlevels * sizeof(NL_Msmgrid_float));
      if (NULL == vgc) return NL_MSM_ERROR_MALLOC;
      pm->qh_f = (NL_Msmgrid_float *) vqh;
      pm->eh_f = (NL_Msmgrid_float *) veh;
      pm->gc_f = (NL_Msmgrid_float *) vgc;
      /* initialize the newest grids appended to array */
      for (level = pm->maxlevels;  level < maxlevels;  level++) {
        GRID_INIT( &(pm->qh_f[level]) );
        GRID_INIT( &(pm->eh_f[level]) );
        GRID_INIT( &(pm->gc_f[level]) );
      }
    }
    else {
      vqh = realloc(pm->qh, maxlevels * sizeof(NL_Msmgrid_double));
      if (NULL == vqh) return NL_MSM_ERROR_MALLOC;
      veh = realloc(pm->eh, maxlevels * sizeof(NL_Msmgrid_double));
      if (NULL == veh) return NL_MSM_ERROR_MALLOC;
      vgc = realloc(pm->gc, maxlevels * sizeof(NL_Msmgrid_double));
      if (NULL == vgc) return NL_MSM_ERROR_MALLOC;
      pm->qh = (NL_Msmgrid_double *) vqh;
      pm->eh = (NL_Msmgrid_double *) veh;
      pm->gc = (NL_Msmgrid_double *) vgc;
      /* initialize the newest grids appended to array */
      for (level = pm->maxlevels;  level < maxlevels;  level++) {
        GRID_INIT( &(pm->qh[level]) );
        GRID_INIT( &(pm->eh[level]) );
        GRID_INIT( &(pm->gc[level]) );
      }
    }
    pm->maxlevels = maxlevels;
  }

  level = 0;
  done = 0;
  alldone = 0;
  do {
    if (issprec) {
      GRID_RESIZE( &(pm->qh_f[level]), float, ia, ni, ja, nj, ka, nk);
      GRID_RESIZE( &(pm->eh_f[level]), float, ia, ni, ja, nj, ka, nk);
    }
    else {
      GRID_RESIZE( &(pm->qh[level]), double, ia, ni, ja, nj, ka, nk);
      GRID_RESIZE( &(pm->eh[level]), double, ia, ni, ja, nj, ka, nk);
    }

    if (++level == nlevels)    done |= 0x07;  /* user limit on levels */

    alldone = (done == 0x07);  /* make sure all dimensions are done */

    if (isclamped) {
      int nelems = ni * nj * nk;
      if (nelems <= lastnelems)  done |= 0x07;
    }

    if (ispx) {
      ni >>= 1;
      ib = ni-1;
      if (ni & 1)              done |= 0x07;  /* == 3 or 1 */
      else if (ni == 2)        done |= 0x01;  /* can do one more */
    }
    else {
      ia = -((-ia+1)/2) - nu;
      ib = (ib+1)/2 + nu;
      ni = ib - ia + 1;
      if (ni <= omega)         done |= 0x01;  /* can do more restrictions */
    }

    if (ispy) {
      nj >>= 1;
      jb = nj-1;
      if (nj & 1)              done |= 0x07;  /* == 3 or 1 */
      else if (nj == 2)        done |= 0x02;  /* can do one more */
    }
    else {
      ja = -((-ja+1)/2) - nu;
      jb = (jb+1)/2 + nu;
      nj = jb - ja + 1;
      if (nj <= omega)         done |= 0x02;  /* can do more restrictions */
    }

    if (ispz) {
      nk >>= 1;
      kb = nk-1;
      if (nk & 1)              done |= 0x07;  /* == 3 or 1 */
      else if (nk == 2)        done |= 0x04;  /* can do one more */
    }
    else {
      ka = -((-ka+1)/2) - nu;
      kb = (kb+1)/2 + nu;
      nk = kb - ka + 1;
      if (nk <= omega)         done |= 0x04;  /* can do more restrictions */
    }

  } while ( ! alldone );
  pm->nlevels = level;

  toplevel = (ispany ? pm->nlevels : pm->nlevels - 1);

  /* ellipsoid axes for grid cutoff weights */
  ni = (int) ceil(2*a/hx) - 1;
  nj = (int) ceil(2*a/hy) - 1;
  nk = (int) ceil(2*a/hz) - 1;
  scaling = 1;  /* initially there is no scaling */
  for (level = 0;  level < toplevel;  level++) {
    if (issprec) {
      p_f = &(pm->gc_f[level]);
      GRID_RESIZE(p_f, float, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
    }
    else {
      p = &(pm->gc[level]);
      GRID_RESIZE(p, double, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
    }

    if (0 == level) {
      index = 0;
      for (k = -nk;  k <= nk;  k++) {
        for (j = -nj;  j <= nj;  j++) {
          for (i = -ni;  i <= ni;  i++) {
            double s, t, gs, gt, g;
            s = sqrt((i*hx)*(i*hx) + (j*hy)*(j*hy) + (k*hz)*(k*hz)) / a;
            t = 0.5 * s;
            if (t >= 1) {
              g = 0;
            }
            else if (s >= 1) {
              gs = 1/s;
              SPOLY(&gt, &d, t, split);
              g = (gs - 0.5 * gt) / a;
            }
            else {
              SPOLY(&gs, &d, s, split);
              SPOLY(&gt, &d, t, split);
              g = (gs - 0.5 * gt) / a;
            }
            if (issprec) {
              GRID_INDEX_CHECK(p_f, i, j, k);
              ASSERT(p_f->buffer + index == p_f->data + GRID_INDEX(p_f,i,j,k));
              p_f->buffer[index] = (float) g;
            }
            else {
              GRID_INDEX_CHECK(p, i, j, k);
              ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
              p->buffer[index] = g;
            }
            index++;
          }
        }
      } /* end loops over k-j-i */
    }
    else {
      /* set each level as scaling of h-level */
      if (issprec) {
        const NL_Msmgrid_float *first = &(pm->gc_f[0]);
        index = 0;
        for (k = -nk;  k <= nk;  k++) {
          for (j = -nj;  j <= nj;  j++) {
            for (i = -ni;  i <= ni;  i++) {
              GRID_INDEX_CHECK(p_f, i, j, k);
              ASSERT(p_f->buffer + index == p_f->data + GRID_INDEX(p_f,i,j,k));
              p_f->buffer[index] = (float) (scaling * first->buffer[index]);
              index++;
            }
          }
        } /* for loops */
      } /* if isprec */
      else {
        const NL_Msmgrid_double *first = &(pm->gc[0]);
        index = 0;
        for (k = -nk;  k <= nk;  k++) {
          for (j = -nj;  j <= nj;  j++) {
            for (i = -ni;  i <= ni;  i++) {
              GRID_INDEX_CHECK(p, i, j, k);
              ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
              p->buffer[index] = scaling * first->buffer[index];
              index++;
            }
          }
        } /* for loops */
      } /* else ! issprec */
    }
    scaling *= 0.5;  /* adjust scaling here to also accommodate toplevel */
  } /* end loop over levels */

  if (toplevel < pm->nlevels) {
    /* nonperiodic in all dimensions,
     * calculate top level weights, ellipsoid axes are length of grid */
    if (issprec) {
      const NL_Msmgrid_float *qhtop_f = &(pm->qh_f[toplevel]);
      ni = qhtop_f->ni - 1;
      nj = qhtop_f->nj - 1;
      nk = qhtop_f->nk - 1;
    }
    else {
      const NL_Msmgrid_double *qhtop = &(pm->qh[toplevel]);
      ni = qhtop->ni - 1;
      nj = qhtop->nj - 1;
      nk = qhtop->nk - 1;
    }
    if (issprec) {
      p_f = &(pm->gc_f[toplevel]);
      GRID_RESIZE(p_f, float, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
    }
    else {
      p = &(pm->gc[toplevel]);
      GRID_RESIZE(p, double, -ni, 2*ni+1, -nj, 2*nj+1, -nk, 2*nk+1);
    }
    index = 0;
    for (k = -nk;  k <= nk;  k++) {
      for (j = -nj;  j <= nj;  j++) {
        for (i = -ni;  i <= ni;  i++) {
          double s, gs;
          s = sqrt((i*hx)*(i*hx) + (j*hy)*(j*hy) + (k*hz)*(k*hz)) / a;
          if (s >= 1) {
            gs = 1/s;
          }
          else {
            SPOLY(&gs, &d, s, split);
          }
          if (issprec) {
            GRID_INDEX_CHECK(p_f, i, j, k);
            ASSERT(p_f->buffer + index == p_f->data + GRID_INDEX(p_f,i,j,k));
            p_f->buffer[index] = (float) (scaling * gs/a);
          }
          else {
            GRID_INDEX_CHECK(p, i, j, k);
            ASSERT( p->buffer + index == p->data + GRID_INDEX(p, i, j, k) );
            p->buffer[index] = scaling * gs/a;
          }
          index++;
        }
      }
    } /* end loops over k-j-i for coarsest level weights */
  }

  /* calculate self energy factor for splitting */
  if (1) {
    double s, gs;
    s = 0;
    SPOLY(&gs, &d, s, split);
    if (issprec) {
      pm->gzero_f = (float) (gs/a);
    }
    else {
      pm->gzero = gs/a;
    }
  }

  return NL_MSM_SUCCESS;
}

setup_hgrid_1d()

int setup_hgrid_1d ( NL_Msm *  pm, double  len, double *  hh, int *  nn, int *  aindex, int *  bindex, int  isperiodic  )

static

Definition at line 658 of file msm_setup.c.

References NL_Msm_t::approx, ASSERT, NL_Msm_t::gridspacing, INTERP_PARAMS, NL_MSM_ERROR_PARAM, NL_MSM_SUCCESS, and InterpParams_t::nu.

Referenced by setup_grids().

      {

  const int nu = INTERP_PARAMS[pm->approx].nu;  /* interp stencil radius */

  ASSERT(hmax > 0);
  if (isperiodic) {
    const double hmin = (4./5) * pm->gridspacing;  /* minimum bound on h */
    const double hmax = 1.5 * hmin;
    double h = len;
    int n = 1;    /* start with one grid point across domain */
    while (h >= hmax) {
      h *= 0.5;   /* halve h */
      n <<= 1;    /* double grid points */
    }
    if (h < hmin) {
      if (n < 4) {  /* either len is too small or hmin is too large */
        return NL_MSM_ERROR_PARAM;
      }
      h *= (4./3);  /* scale h by 4/3 */
      n >>= 2;      /* scale n by 3/4 */
      n *= 3;
    }
    /* now we have:  hmin <= h < hmax */
    /* now we have:  n is power of two times no more than one power of 3 */
    *hh = h;
    *nn = n;
    *aindex = 0;
    *bindex = n-1;
  }
  else {  /* non-periodic */
    double h = pm->gridspacing;
    int n = (int) floorf(len / h) + 1;
    *hh = h;
    *nn = n;
    *aindex = -nu;
    *bindex = n + nu;
  }
  return NL_MSM_SUCCESS;
}

Variable Documentation

INTERP_PARAMS

InterpParams INTERP_PARAMS[]

static

Initial value:

= {
  { 1, 4, 6 },    
  { 2, 6, 10 },   
  { 2, 6, 10 },   
  { 3, 8, 14 },   
  { 3, 8, 14 },   
  { 4, 10, 18 },  
  { 4, 10, 18 },  
  { 1, 4, 6 },    
}

Definition at line 206 of file msm_setup.c.

Referenced by setup_grids(), and setup_hgrid_1d().