NAMD
Public Types | Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | Friends | List of all members
ComputeMsmMgr Class Reference
Inheritance diagram for ComputeMsmMgr:

Public Types

enum  {
  VXX =0, VXY, VXZ, VYY,
  VYZ, VZZ, VMAX
}
 
enum  Approx {
  CUBIC =0, QUINTIC, QUINTIC2, SEPTIC,
  SEPTIC3, NONIC, NONIC4, C1HERMITE,
  NUM_APPROX
}
 
enum  Split {
  TAYLOR2 =0, TAYLOR3, TAYLOR4, TAYLOR5,
  TAYLOR6, TAYLOR7, TAYLOR8, TAYLOR2_DISP,
  TAYLOR3_DISP, TAYLOR4_DISP, TAYLOR5_DISP, TAYLOR6_DISP,
  TAYLOR7_DISP, TAYLOR8_DISP, NUM_SPLIT
}
 
enum  { MAX_POLY_DEGREE = 9, MAX_NSTENCIL_SIZE = (2*MAX_POLY_DEGREE + 1), MAX_NSTENCIL_SKIP_ZERO = (MAX_POLY_DEGREE + 2), NUM_APPROX_FORMS = (NONIC4 - CUBIC) + 1 }
 

Public Member Functions

 ComputeMsmMgr ()
 
 ~ComputeMsmMgr ()
 
void initialize (MsmInitMsg *)
 
void initialize_create ()
 
void recvMsmBlockProxy (MsmBlockProxyMsg *)
 
void recvMsmGridCutoffProxy (MsmGridCutoffProxyMsg *)
 
void recvMsmC1HermiteBlockProxy (MsmC1HermiteBlockProxyMsg *)
 
void recvMsmC1HermiteGridCutoffProxy (MsmC1HermiteGridCutoffProxyMsg *)
 
void update (CkQdMsg *)
 
void compute (msm::Array< int > &patchIDList)
 
void addPotential (GridMsg *)
 
void doneCompute ()
 
void setCompute (ComputeMsm *c)
 
msm::PatchPtrArraypatchPtrArray ()
 
msm::MapmapData ()
 
int numLevels () const
 
void setup_hgrid_1d (BigReal len, BigReal &hh, int &nn, int &ia, int &ib, int isperiodic)
 
void setup_periodic_blocksize (int &bsize, int n)
 
int blockFlatIndex (int level, int i, int j, int k)
 
float calcBlockWork (const msm::BlockDiagram &b)
 
float calcGcutWork (const msm::BlockSend &bs)
 
void initVirialContrib ()
 
void addVirialContrib ()
 
void subtractVirialContrib ()
 
void doneVirialContrib ()
 
void stencil_1d (Float phi[], Float t)
 
void d_stencil_1d (Float dphi[], Float phi[], Float t, Float h_1)
 
void stencil_1d_c1hermite (Float phi[], Float psi[], Float t, Float h)
 
void d_stencil_1d_c1hermite (Float dphi[], Float phi[], Float dpsi[], Float psi[], Float t, Float h, Float h_1)
 

Static Public Member Functions

static int sign (int n)
 
static void splitting (BigReal &g, BigReal &dg, BigReal r_a, int _split)
 
static void ndsplitting (BigReal pg[], BigReal s, int n, int _split)
 
static void gc_c1hermite_elem_accum (C1Matrix &matrix, BigReal _c, Vector rv, BigReal _a, int _split)
 

Public Attributes

CProxy_ComputeMsmMgr msmProxy
 
ComputeMsmmsmCompute
 
msm::Array< CProxy_MsmBlock > msmBlock
 
msm::Array
< CProxy_MsmC1HermiteBlock > 
msmC1HermiteBlock
 
CProxy_MsmGridCutoff msmGridCutoff
 
CProxy_MsmC1HermiteGridCutoff msmC1HermiteGridCutoff
 
int numGridCutoff
 
msm::Map map
 
msm::PatchPtrArray patchPtr
 
msm::Grid< Floatsubgrid
 
msm::Grid< C1Vectorsubgrid_c1hermite
 
msm::Array< int > blockAssign
 
msm::Array< int > gcutAssign
 
msm::Grid< Floatgvsum
 
int numVirialContrib
 
int cntVirialContrib
 
Float virial [VMAX]
 
Vector c
 
Vector u
 
Vector v
 
Vector w
 
Vector ru
 
Vector rv
 
Vector rw
 
int ispu
 
int ispv
 
int ispw
 
Lattice lattice
 
ScaledPosition smin
 
ScaledPosition smax
 
BigReal gridspacing
 
BigReal padding
 
BigReal gridScalingFactor
 
BigReal a
 
BigReal hxlen
 
BigReal hylen
 
BigReal hzlen
 
BigReal hxlen_1
 
BigReal hylen_1
 
BigReal hzlen_1
 
Vector hu
 
Vector hv
 
Vector hw
 
Float hufx
 
Float hufy
 
Float hufz
 
Float hvfx
 
Float hvfy
 
Float hvfz
 
Float hwfx
 
Float hwfy
 
Float hwfz
 
int nhx
 
int nhy
 
int nhz
 
int approx
 
int split
 
int nlevels
 
int dispersion
 
BigReal gzero
 
Vector sglower
 
BigReal shx
 
BigReal shy
 
BigReal shz
 
BigReal shx_1
 
BigReal shy_1
 
BigReal shz_1
 
Vector sx_shx
 
Vector sy_shy
 
Vector sz_shz
 
Float srx_x
 
Float srx_y
 
Float srx_z
 
Float sry_x
 
Float sry_y
 
Float sry_z
 
Float srz_x
 
Float srz_y
 
Float srz_z
 
int s_edge
 
int omega
 

Static Public Attributes

static const int PolyDegree [NUM_APPROX]
 
static const int Nstencil [NUM_APPROX]
 
static const int IndexOffset [NUM_APPROX][MAX_NSTENCIL_SKIP_ZERO]
 
static const Float PhiStencil [NUM_APPROX_FORMS][MAX_NSTENCIL_SKIP_ZERO]
 

Friends

struct msm::PatchData
 
class MsmBlock
 
class MsmBlockMap
 
class MsmGridCutoffMap
 

Detailed Description

Definition at line 364 of file ComputeMsm.C.

Member Enumeration Documentation

anonymous enum
Enumerator
VXX 
VXY 
VXZ 
VYY 
VYZ 
VZZ 
VMAX 

Definition at line 526 of file ComputeMsm.C.

anonymous enum
Enumerator
MAX_POLY_DEGREE 
MAX_NSTENCIL_SIZE 
MAX_NSTENCIL_SKIP_ZERO 
NUM_APPROX_FORMS 

Definition at line 633 of file ComputeMsm.C.

633  {
634  // Approximation formulas with up to degree 9 polynomials.
635  MAX_POLY_DEGREE = 9,
636 
637  // Max stencil length for polynomial approximation.
639 
640  // Max stencil length when skipping zeros
641  // (almost half entries are zero for interpolating polynomials).
643 
644  // Number of scalar approximation formulaes
645  NUM_APPROX_FORMS = (NONIC4 - CUBIC) + 1
646  };
Enumerator
CUBIC 
QUINTIC 
QUINTIC2 
SEPTIC 
SEPTIC3 
NONIC 
NONIC4 
C1HERMITE 
NUM_APPROX 

Definition at line 625 of file ComputeMsm.C.

Constructor & Destructor Documentation

ComputeMsmMgr::ComputeMsmMgr ( )

Definition at line 3850 of file ComputeMsm.C.

3850  :
3851  msmProxy(thisgroup), msmCompute(0)
3852 {
3853 #ifdef DEBUG_MSM_VERBOSE
3854  printf("ComputeMsmMgr: (constructor) PE %d\n", CkMyPe());
3855 #endif
3856  CkpvAccess(BOCclass_group).computeMsmMgr = thisgroup;
3857 
3858 #ifdef MSM_TIMING
3859  if (CkMyPe() == 0) {
3860  msmTimer = CProxy_MsmTimer::ckNew();
3861  }
3862  initTiming();
3863 #endif
3864 #ifdef MSM_PROFILING
3865  if (CkMyPe() == 0) {
3866  msmProfiler = CProxy_MsmProfiler::ckNew();
3867  }
3868  initProfiling();
3869 #endif
3870 }
ComputeMsm * msmCompute
Definition: ComputeMsm.C:462
CProxy_ComputeMsmMgr msmProxy
Definition: ComputeMsm.C:461
ComputeMsmMgr::~ComputeMsmMgr ( )

Definition at line 3872 of file ComputeMsm.C.

3873 {
3874 #ifdef DEBUG_MSM_VERBOSE
3875  printf("ComputeMsmMgr: (destructor) PE %d\n", CkMyPe());
3876 #endif
3877  // free memory?
3878 }

Member Function Documentation

void ComputeMsmMgr::addPotential ( GridMsg gm)

Definition at line 6014 of file ComputeMsm.C.

References approx, C1HERMITE, GridMsg::get(), NAMD_die(), patchPtr, subgrid, and subgrid_c1hermite.

Referenced by MsmBlock::sendDownPotential(), MsmC1HermiteBlock::sendDownPotential(), MsmBlock::sumReducedPotential(), and MsmC1HermiteBlock::sumReducedPotential().

6015 {
6016  int pid; // receive patch ID
6017  int pseq;
6018  if (approx == C1HERMITE) {
6019  gm->get(subgrid_c1hermite, pid, pseq);
6020  }
6021  else {
6022  gm->get(subgrid, pid, pseq);
6023  }
6024  delete gm;
6025  if (patchPtr[pid] == NULL) {
6026  char msg[100];
6027  snprintf(msg, sizeof(msg), "Expecting patch %d to exist on PE %d",
6028  pid, CkMyPe());
6029  NAMD_die(msg);
6030  }
6031  if (approx == C1HERMITE) {
6032  patchPtr[pid]->addPotentialC1Hermite(subgrid_c1hermite);
6033  }
6034  else {
6035  patchPtr[pid]->addPotential(subgrid);
6036  }
6037 }
msm::PatchPtrArray patchPtr
Definition: ComputeMsm.C:476
msm::Grid< Float > subgrid
Definition: ComputeMsm.C:480
void NAMD_die(const char *err_msg)
Definition: common.C:83
msm::Grid< C1Vector > subgrid_c1hermite
Definition: ComputeMsm.C:481
void get(msm::Grid< T > &g, int &id, int &seq)
Definition: ComputeMsm.C:141
void ComputeMsmMgr::addVirialContrib ( )
inline

Definition at line 533 of file ComputeMsm.C.

References numVirialContrib.

533  {
535  }
int numVirialContrib
Definition: ComputeMsm.C:524
int ComputeMsmMgr::blockFlatIndex ( int  level,
int  i,
int  j,
int  k 
)
inline

Definition at line 487 of file ComputeMsm.C.

References msm::Map::blockLevel, and map.

487  {
488  int n = 0;
489  for (int l = 0; l < level; l++) {
490  n += map.blockLevel[l].nn();
491  }
492  return (n + map.blockLevel[level].flatindex(i,j,k));
493  }
Array< Grid< BlockDiagram > > blockLevel
Definition: MsmMap.h:956
msm::Map map
Definition: ComputeMsm.C:471
float ComputeMsmMgr::calcBlockWork ( const msm::BlockDiagram b)
inline

Definition at line 494 of file ComputeMsm.C.

References approx, msm::Map::bsx, msm::Map::bsy, msm::Map::bsz, C1HERMITE, msm::IndexRange::extent(), msm::Map::gc, msm::Map::gc_c1hermite, msm::Ivec::i, msm::Ivec::j, msm::Ivec::k, map, msm::BlockDiagram::nrange, and msm::BlockDiagram::nrangeCutoff.

494  {
495  // XXX ratio of work for MsmBlock to MsmGridCutoff?
496  const float scalingFactor = 3;
497  const int volumeFullBlock = map.bsx[0] * map.bsy[0] * map.bsz[0];
498  msm::Ivec gn;
499  if (approx == C1HERMITE) {
500  gn = map.gc_c1hermite[0].extent();
501  }
502  else {
503  gn = map.gc[0].extent();
504  }
505  const int volumeFullCutoff = (map.bsx[0] + gn.i - 1) *
506  (map.bsy[0] + gn.j - 1) * (map.bsz[0] + gn.k - 1);
507  msm::Ivec n = b.nrange.extent();
508  int volumeBlock = n.i * n.j * n.k;
509  msm::Ivec nc = b.nrangeCutoff.extent();
510  int volumeCutoff = nc.i * nc.j * nc.k;
511  return( scalingFactor * (float(volumeBlock) / volumeFullBlock) *
512  (float(volumeCutoff) / volumeFullCutoff) );
513  }
Array< int > bsz
Definition: MsmMap.h:960
int k
Definition: MsmMap.h:411
Array< int > bsy
Definition: MsmMap.h:960
int i
Definition: MsmMap.h:411
Ivec extent() const
Definition: MsmMap.h:445
int j
Definition: MsmMap.h:411
IndexRange nrange
Definition: MsmMap.h:911
Array< int > bsx
Definition: MsmMap.h:960
IndexRange nrangeCutoff
Definition: MsmMap.h:912
Array< Grid< Float > > gc
Definition: MsmMap.h:944
msm::Map map
Definition: ComputeMsm.C:471
Array< Grid< C1Matrix > > gc_c1hermite
Definition: MsmMap.h:949
float ComputeMsmMgr::calcGcutWork ( const msm::BlockSend bs)
inline

Definition at line 514 of file ComputeMsm.C.

References msm::Map::bsx, msm::Map::bsy, msm::Map::bsz, msm::IndexRange::extent(), msm::Ivec::i, msm::Ivec::j, msm::Ivec::k, map, and msm::BlockSend::nrange_wrap.

514  {
515  const int volumeFullBlock = map.bsx[0] * map.bsy[0] * map.bsz[0];
516  msm::Ivec n = bs.nrange_wrap.extent();;
517  int volumeBlock = n.i * n.j * n.k;
518  return( float(volumeBlock) / volumeFullBlock );
519  }
Array< int > bsz
Definition: MsmMap.h:960
int k
Definition: MsmMap.h:411
Array< int > bsy
Definition: MsmMap.h:960
int i
Definition: MsmMap.h:411
Ivec extent() const
Definition: MsmMap.h:445
int j
Definition: MsmMap.h:411
Array< int > bsx
Definition: MsmMap.h:960
IndexRange nrange_wrap
Definition: MsmMap.h:870
msm::Map map
Definition: ComputeMsm.C:471
void ComputeMsmMgr::compute ( msm::Array< int > &  patchIDList)

Definition at line 5986 of file ComputeMsm.C.

References approx, C1HERMITE, msm::Array< T >::len(), NAMD_die(), and patchPtr.

Referenced by ComputeMsm::doWork().

5987 {
5988 #ifdef DEBUG_MSM_VERBOSE
5989  printf("ComputeMsmMgr: compute() PE=%d\n", CkMyPe());
5990 #endif
5991 
5992  int n;
5993  for (n = 0; n < patchIDList.len(); n++) {
5994  int patchID = patchIDList[n];
5995  if (patchPtr[patchID] == NULL) {
5996  char msg[100];
5997  snprintf(msg, sizeof(msg),
5998  "Expected MSM data for patch %d does not exist on PE %d",
5999  patchID, CkMyPe());
6000  NAMD_die(msg);
6001  }
6002  if (approx == C1HERMITE) {
6003  patchPtr[patchID]->anterpolationC1Hermite();
6004  }
6005  else {
6006  patchPtr[patchID]->anterpolation();
6007  }
6008  // all else should follow from here
6009  }
6010  return;
6011 }
int len() const
Definition: MsmMap.h:218
msm::PatchPtrArray patchPtr
Definition: ComputeMsm.C:476
void NAMD_die(const char *err_msg)
Definition: common.C:83
void ComputeMsmMgr::d_stencil_1d ( Float  dphi[],
Float  phi[],
Float  t,
Float  h_1 
)
inline

Definition at line 937 of file ComputeMsm.C.

References approx, CUBIC, NAMD_die(), NONIC, NONIC4, QUINTIC, QUINTIC2, SEPTIC, and SEPTIC3.

Referenced by msm::PatchData::interpolation().

937  {
938  switch (approx) {
939  case CUBIC:
940  phi[0] = 0.5f * (1 - t) * (2 - t) * (2 - t);
941  dphi[0] = (1.5f * t - 2) * (2 - t) * h_1;
942  t--;
943  phi[1] = (1 - t) * (1 + t - 1.5f * t * t);
944  dphi[1] = (-5 + 4.5f * t) * t * h_1;
945  t--;
946  phi[2] = (1 + t) * (1 - t - 1.5f * t * t);
947  dphi[2] = (-5 - 4.5f * t) * t * h_1;
948  t--;
949  phi[3] = 0.5f * (1 + t) * (2 + t) * (2 + t);
950  dphi[3] = (1.5f * t + 2) * (2 + t) * h_1;
951  break;
952  case QUINTIC:
953  phi[0] = (1.f/24) * (1-t) * (2-t) * (3-t) * (3-t) * (4-t);
954  dphi[0] = ((-1.f/24) * ((3-t) * (3-t) * (14 + t * (-14 + 3*t))
955  + 2 * (1-t) * (2-t) * (3-t) * (4-t))) * h_1;
956  t--;
957  phi[1] = (1-t) * (2-t) * (3-t) * ((1.f/6)
958  + t * (0.375f - (5.f/24)*t));
959  dphi[1] = (-((1.f/6) + t * (0.375f - (5.f/24)*t)) *
960  (11 + t * (-12 + 3*t)) + (1-t) * (2-t) * (3-t) *
961  (0.375f - (5.f/12)*t)) * h_1;
962  t--;
963  phi[2] = (1-t*t) * (2-t) * (0.5f + t * (0.25f - (5.f/12)*t));
964  dphi[2] = (-(0.5f + t * (0.25f - (5.f/12)*t)) * (1 + t * (4 - 3*t))
965  + (1-t*t) * (2-t) * (0.25f - (5.f/6)*t)) * h_1;
966  t--;
967  phi[3] = (1-t*t) * (2+t) * (0.5f - t * (0.25f + (5.f/12)*t));
968  dphi[3] = ((0.5f + t * (-0.25f - (5.f/12)*t)) * (1 + t * (-4 - 3*t))
969  - (1-t*t) * (2+t) * (0.25f + (5.f/6)*t)) * h_1;
970  t--;
971  phi[4] = (1+t) * (2+t) * (3+t) * ((1.f/6)
972  - t * (0.375f + (5.f/24)*t));
973  dphi[4] = (((1.f/6) + t * (-0.375f - (5.f/24)*t)) *
974  (11 + t * (12 + 3*t)) - (1+t) * (2+t) * (3+t) *
975  (0.375f + (5.f/12)*t)) * h_1;
976  t--;
977  phi[5] = (1.f/24) * (1+t) * (2+t) * (3+t) * (3+t) * (4+t);
978  dphi[5] = ((1.f/24) * ((3+t) * (3+t) * (14 + t * (14 + 3*t))
979  + 2 * (1+t) * (2+t) * (3+t) * (4+t))) * h_1;
980  break;
981  case QUINTIC2:
982  phi[0] = (1.f/24) * (3-t) * (3-t) * (3-t) * (t-2) * (5*t-8);
983  dphi[0] = ((1.f/24) * (3-t) * (3-t) * ((3-t)*(5*t-8)
984  - 3*(t-2)*(5*t-8) + 5*(t-2)*(3-t))) * h_1;
985  t--;
986  phi[1] = (-1.f/24) * (2-t) * (t-1) * (-48+t*(153+t*(-114+t*25)));
987  dphi[1] = ((-1.f/24) * ((2-t)*(-48+t*(153+t*(-114+t*25)))
988  - (t-1)* (-48+t*(153+t*(-114+t*25)))
989  + (2-t)*(t-1)*(153+t*(-228+t*75)))) * h_1;
990  t--;
991  phi[2] = (1.f/12) * (1-t) * (12+t*(12+t*(-3+t*(-38+t*25))));
992  dphi[2] = ((1.f/12) * (-(12+t*(12+t*(-3+t*(-38+t*25))))
993  + (1-t)*(12+t*(-6+t*(-114+t*100))))) * h_1;
994  t--;
995  phi[3] = (1.f/12) * (1+t) * (12+t*(-12+t*(-3+t*(38+t*25))));
996  dphi[3] = ((1.f/12) * ((12+t*(-12+t*(-3+t*(38+t*25))))
997  + (1+t)*(-12+t*(-6+t*(114+t*100))))) * h_1;
998  t--;
999  phi[4] = (-1.f/24) * (2+t) * (t+1) * (48+t*(153+t*(114+t*25)));
1000  dphi[4] = ((-1.f/24) * ((2+t)*(48+t*(153+t*(114+t*25)))
1001  + (t+1)* (48+t*(153+t*(114+t*25)))
1002  + (2+t)*(t+1)*(153+t*(228+t*75)))) * h_1;
1003  t--;
1004  phi[5] = (1.f/24) * (3+t) * (3+t) * (3+t) * (t+2) * (5*t+8);
1005  dphi[5] = ((1.f/24) * (3+t) * (3+t) * ((3+t)*(5*t+8)
1006  + 3*(t+2)*(5*t+8) + 5*(t+2)*(3+t))) * h_1;
1007  break;
1008  case SEPTIC:
1009  phi[0] = (-1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-4)*(t-5)*(t-6);
1010  dphi[0] = (-1.f/720)*(t-4)*(-1944+t*(3644+t*(-2512+t*(807
1011  +t*(-122+t*7))))) * h_1;
1012  t--;
1013  phi[1] = (1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-5)*(-6+t*(-20+7*t));
1014  dphi[1] = (1.f/720)*(756+t*(-9940+t*(17724+t*(-12740+t*(4445
1015  +t*(-750+t*49)))))) * h_1;
1016  t--;
1017  phi[2] = (-1.f/240)*(t*t-1)*(t-2)*(t-3)*(t-4)*(-10+t*(-12+7*t));
1018  dphi[2] = (-1.f/240)*(-28+t*(1260+t*(-756+t*(-1260+t*(1365
1019  +t*(-450+t*49)))))) * h_1;
1020  t--;
1021  phi[3] = (1.f/144)*(t*t-1)*(t*t-4)*(t-3)*(-12+t*(-4+7*t));
1022  dphi[3] = (1.f/144)*t*(-560+t*(84+t*(644+t*(-175
1023  +t*(-150+t*49))))) * h_1;
1024  t--;
1025  phi[4] = (-1.f/144)*(t*t-1)*(t*t-4)*(t+3)*(-12+t*(4+7*t));
1026  dphi[4] = (-1.f/144)*t*(560+t*(84+t*(-644+t*(-175
1027  +t*(150+t*49))))) * h_1;
1028  t--;
1029  phi[5] = (1.f/240)*(t*t-1)*(t+2)*(t+3)*(t+4)*(-10+t*(12+7*t));
1030  dphi[5] = (1.f/240)*(-28+t*(-1260+t*(-756+t*(1260+t*(1365
1031  +t*(450+t*49)))))) * h_1;
1032  t--;
1033  phi[6] = (-1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(-6+t*(20+7*t));
1034  dphi[6] = (-1.f/720)*(756+t*(9940+t*(17724+t*(12740+t*(4445
1035  +t*(750+t*49)))))) * h_1;
1036  t--;
1037  phi[7] = (1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+4)*(t+5)*(t+6);
1038  dphi[7] = (1.f/720)*(t+4)*(1944+t*(3644+t*(2512+t*(807
1039  +t*(122+t*7))))) * h_1;
1040  break;
1041  case SEPTIC3:
1042  phi[0] = (3632.f/5) + t*((-7456.f/5) + t*((58786.f/45) + t*(-633
1043  + t*((26383.f/144) + t*((-22807.f/720) + t*((727.f/240)
1044  + t*(-89.f/720)))))));
1045  dphi[0] = ((-7456.f/5) + t*((117572.f/45) + t*(-1899
1046  + t*((26383.f/36) + t*((-22807.f/144) + t*((727.f/40)
1047  + t*(-623.f/720))))))) * h_1;
1048  t--;
1049  phi[1] = -440 + t*((25949.f/20) + t*((-117131.f/72) + t*((2247.f/2)
1050  + t*((-66437.f/144) + t*((81109.f/720) + t*((-727.f/48)
1051  + t*(623.f/720)))))));
1052  dphi[1] = ((25949.f/20) + t*((-117131.f/36) + t*((6741.f/2)
1053  + t*((-66437.f/36) + t*((81109.f/144) + t*((-727.f/8)
1054  + t*(4361.f/720))))))) * h_1;
1055  t--;
1056  phi[2] = (138.f/5) + t*((-8617.f/60) + t*((12873.f/40) + t*((-791.f/2)
1057  + t*((4557.f/16) + t*((-9583.f/80) + t*((2181.f/80)
1058  + t*(-623.f/240)))))));
1059  dphi[2] = ((-8617.f/60) + t*((12873.f/20) + t*((-2373.f/2)
1060  + t*((4557.f/4) + t*((-9583.f/16) + t*((6543.f/40)
1061  + t*(-4361.f/240))))))) * h_1;
1062  t--;
1063  phi[3] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((2569.f/144)
1064  + t*((-727.f/48) + t*(623.f/144)))));
1065  dphi[3] = (t*((-49.f/18) + t*t*((-959.f/36) + t*((12845.f/144)
1066  + t*((-727.f/8) + t*(4361.f/144)))))) * h_1;
1067  t--;
1068  phi[4] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((-2569.f/144)
1069  + t*((-727.f/48) + t*(-623.f/144)))));
1070  dphi[4] = (t*((-49.f/18) + t*t*((-959.f/36) + t*((-12845.f/144)
1071  + t*((-727.f/8) + t*(-4361.f/144)))))) * h_1;
1072  t--;
1073  phi[5] = (138.f/5) + t*((8617.f/60) + t*((12873.f/40) + t*((791.f/2)
1074  + t*((4557.f/16) + t*((9583.f/80) + t*((2181.f/80)
1075  + t*(623.f/240)))))));
1076  dphi[5] = ((8617.f/60) + t*((12873.f/20) + t*((2373.f/2)
1077  + t*((4557.f/4) + t*((9583.f/16) + t*((6543.f/40)
1078  + t*(4361.f/240))))))) * h_1;
1079  t--;
1080  phi[6] = -440 + t*((-25949.f/20) + t*((-117131.f/72) + t*((-2247.f/2)
1081  + t*((-66437.f/144) + t*((-81109.f/720) + t*((-727.f/48)
1082  + t*(-623.f/720)))))));
1083  dphi[6] = ((-25949.f/20) + t*((-117131.f/36) + t*((-6741.f/2)
1084  + t*((-66437.f/36) + t*((-81109.f/144) + t*((-727.f/8)
1085  + t*(-4361.f/720))))))) * h_1;
1086  t--;
1087  phi[7] = (3632.f/5) + t*((7456.f/5) + t*((58786.f/45) + t*(633
1088  + t*((26383.f/144) + t*((22807.f/720) + t*((727.f/240)
1089  + t*(89.f/720)))))));
1090  dphi[7] = ((7456.f/5) + t*((117572.f/45) + t*(1899
1091  + t*((26383.f/36) + t*((22807.f/144) + t*((727.f/40)
1092  + t*(623.f/720))))))) * h_1;
1093  break;
1094  case NONIC:
1095  phi[0] = (-1.f/40320)*(t-8)*(t-7)*(t-6)*(t-5)*(t-5)*(t-4)*(t-3)*
1096  (t-2)*(t-1);
1097  dphi[0] = (-1.f/40320)*(t-5)*(-117648+t*(256552+t*(-221416
1098  +t*(99340+t*(-25261+t*(3667+t*(-283+t*9)))))))*h_1;
1099  t--;
1100  phi[1] = (1.f/40320)*(t-7)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*
1101  (-8+t*(-35+9*t));
1102  dphi[1] = (1.f/40320)*(71856+t*(-795368+t*(1569240+t*(-1357692
1103  +t*(634725+t*(-172116+t*(27090+t*(-2296+t*81))))))))*h_1;
1104  t--;
1105  phi[2] = (-1.f/10080)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*
1106  (-14+t*(-25+9*t));
1107  dphi[2] = (1.f/10080)*(3384+t*(-69080+t*(55026
1108  +t*(62580+t*(-99225+t*(51660+t*(-13104+t*(1640
1109  +t*(-81)))))))))*h_1;
1110  t--;
1111  phi[3] = (1.f/1440)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*
1112  (-6+t*(-5+3*t));
1113  dphi[3] = (1.f/1440)*(72+t*(-6344+t*(2070
1114  +t*(7644+t*(-4725+t*(-828+t*(1260+t*(-328+t*27))))))))*h_1;
1115  t--;
1116  phi[4] = (-1.f/2880)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*
1117  (-20+t*(-5+9*t));
1118  dphi[4] = (-1.f/2880)*t*(10792+t*(-972+t*(-12516
1119  +t*(2205+t*(3924+t*(-882+t*(-328+t*81)))))))*h_1;
1120  t--;
1121  phi[5] = (1.f/2880)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*
1122  (-20+t*(5+9*t));
1123  dphi[5] = (1.f/2880)*t*(-10792+t*(-972+t*(12516
1124  +t*(2205+t*(-3924+t*(-882+t*(328+t*81)))))))*h_1;
1125  t--;
1126  phi[6] = (-1.f/1440)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*
1127  (-6+t*(5+3*t));
1128  dphi[6] = (1.f/1440)*(-72+t*(-6344+t*(-2070
1129  +t*(7644+t*(4725+t*(-828+t*(-1260+t*(-328+t*(-27)))))))))*h_1;
1130  t--;
1131  phi[7] = (1.f/10080)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*
1132  (-14+t*(25+9*t));
1133  dphi[7] = (1.f/10080)*(-3384+t*(-69080+t*(-55026
1134  +t*(62580+t*(99225+t*(51660+t*(13104+t*(1640+t*81))))))))*h_1;
1135  t--;
1136  phi[8] = (-1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*(t+7)*
1137  (-8+t*(35+9*t));
1138  dphi[8] = (-1.f/40320)*(71856+t*(795368+t*(1569240
1139  +t*(1357692+t*(634725+t*(172116+t*(27090+t*(2296
1140  +t*81))))))))*h_1;
1141  t--;
1142  phi[9] = (1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+5)*(t+6)*
1143  (t+7)*(t+8);
1144  dphi[9] = (1.f/40320)*(t+5)*(117648+t*(256552+t*(221416
1145  +t*(99340+t*(25261+t*(3667+t*(283+t*9)))))))*h_1;
1146  break;
1147  case NONIC4:
1148  { // begin grouping to define local variables
1149  double Tphi[10], Tdphi[10];
1150  double T=t;
1151  Tphi[0] = 439375./7+T*(-64188125./504+T*(231125375./2016
1152  +T*(-17306975./288+T*(7761805./384+T*(-2895587./640
1153  +T*(129391./192+T*(-259715./4032+T*(28909./8064
1154  +T*(-3569./40320)))))))));
1155  Tdphi[0] = (-64188125./504+T*(231125375./1008
1156  +T*(-17306975./96+T*(7761805./96+T*(-2895587./128
1157  +T*(129391./32+T*(-259715./576+T*(28909./1008
1158  +T*(-3569./4480))))))))) * h_1;
1159  T--;
1160  Tphi[1] = -56375+T*(8314091./56+T*(-49901303./288+T*(3763529./32
1161  +T*(-19648027./384+T*(9469163./640+T*(-545977./192
1162  +T*(156927./448+T*(-28909./1152
1163  +T*(3569./4480)))))))));
1164  Tdphi[1] = (8314091./56+T*(-49901303./144+T*(11290587./32
1165  +T*(-19648027./96+T*(9469163./128+T*(-545977./32
1166  +T*(156927./64+T*(-28909./144
1167  +T*(32121./4480))))))))) * h_1;
1168  T--;
1169  Tphi[2] = 68776./7+T*(-1038011./28+T*(31157515./504+T*(-956669./16
1170  +T*(3548009./96+T*(-2422263./160+T*(197255./48
1171  +T*(-19959./28+T*(144545./2016
1172  +T*(-3569./1120)))))))));
1173  Tdphi[2] = (-1038011./28+T*(31157515./252+T*(-2870007./16
1174  +T*(3548009./24+T*(-2422263./32+T*(197255./8
1175  +T*(-19959./4+T*(144545./252
1176  +T*(-32121./1120))))))))) * h_1;
1177  T--;
1178  Tphi[3] = -154+T*(12757./12+T*(-230123./72+T*(264481./48
1179  +T*(-576499./96+T*(686147./160+T*(-96277./48
1180  +T*(14221./24+T*(-28909./288+T*(3569./480)))))))));
1181  Tdphi[3] = (12757./12+T*(-230123./36+T*(264481./16
1182  +T*(-576499./24+T*(686147./32+T*(-96277./8
1183  +T*(99547./24+T*(-28909./36
1184  +T*(10707./160))))))))) * h_1;
1185  T--;
1186  Tphi[4] = 1+T*T*(-205./144+T*T*(91./192+T*(-6181./320
1187  +T*(6337./96+T*(-2745./32+T*(28909./576
1188  +T*(-3569./320)))))));
1189  Tdphi[4] = T*(-205./72+T*T*(91./48+T*(-6181./64
1190  +T*(6337./16+T*(-19215./32+T*(28909./72
1191  +T*(-32121./320))))))) * h_1;
1192  T--;
1193  Tphi[5] = 1+T*T*(-205./144+T*T*(91./192+T*(6181./320
1194  +T*(6337./96+T*(2745./32+T*(28909./576
1195  +T*(3569./320)))))));
1196  Tdphi[5] = T*(-205./72+T*T*(91./48+T*(6181./64
1197  +T*(6337./16+T*(19215./32+T*(28909./72
1198  +T*(32121./320))))))) * h_1;
1199  T--;
1200  Tphi[6] = -154+T*(-12757./12+T*(-230123./72+T*(-264481./48
1201  +T*(-576499./96+T*(-686147./160+T*(-96277./48
1202  +T*(-14221./24+T*(-28909./288+T*(-3569./480)))))))));
1203  Tdphi[6] = (-12757./12+T*(-230123./36+T*(-264481./16
1204  +T*(-576499./24+T*(-686147./32+T*(-96277./8
1205  +T*(-99547./24+T*(-28909./36
1206  +T*(-10707./160))))))))) * h_1;
1207  T--;
1208  Tphi[7] = 68776./7+T*(1038011./28+T*(31157515./504+T*(956669./16
1209  +T*(3548009./96+T*(2422263./160+T*(197255./48
1210  +T*(19959./28+T*(144545./2016+T*(3569./1120)))))))));
1211  Tdphi[7] = (1038011./28+T*(31157515./252+T*(2870007./16
1212  +T*(3548009./24+T*(2422263./32+T*(197255./8
1213  +T*(19959./4+T*(144545./252
1214  +T*(32121./1120))))))))) * h_1;
1215  T--;
1216  Tphi[8] = -56375+T*(-8314091./56+T*(-49901303./288+T*(-3763529./32
1217  +T*(-19648027./384+T*(-9469163./640+T*(-545977./192
1218  +T*(-156927./448+T*(-28909./1152
1219  +T*(-3569./4480)))))))));
1220  Tdphi[8] = (-8314091./56+T*(-49901303./144+T*(-11290587./32
1221  +T*(-19648027./96+T*(-9469163./128+T*(-545977./32
1222  +T*(-156927./64+T*(-28909./144
1223  +T*(-32121./4480))))))))) * h_1;
1224  T--;
1225  Tphi[9] = 439375./7+T*(64188125./504+T*(231125375./2016
1226  +T*(17306975./288+T*(7761805./384+T*(2895587./640
1227  +T*(129391./192+T*(259715./4032+T*(28909./8064
1228  +T*(3569./40320)))))))));
1229  Tdphi[9] = (64188125./504+T*(231125375./1008
1230  +T*(17306975./96+T*(7761805./96+T*(2895587./128
1231  +T*(129391./32+T*(259715./576+T*(28909./1008
1232  +T*(3569./4480))))))))) * h_1;
1233  for (int i=0; i < 10; i++) {
1234  phi[i] = Float(Tphi[i]);
1235  dphi[i] = Float(Tdphi[i]);
1236  }
1237  } // end grouping to define local variables
1238  break;
1239  default:
1240  NAMD_die("Unknown MSM approximation.");
1241  } // switch
1242  } // d_stencil_1d()
void NAMD_die(const char *err_msg)
Definition: common.C:83
float Float
Definition: MsmMap.h:74
void ComputeMsmMgr::d_stencil_1d_c1hermite ( Float  dphi[],
Float  phi[],
Float  dpsi[],
Float  psi[],
Float  t,
Float  h,
Float  h_1 
)
inline

Definition at line 1252 of file ComputeMsm.C.

Referenced by msm::PatchData::interpolationC1Hermite().

1254  {
1255  phi[0] = (1 - t) * (1 - t) * (1 + 2*t);
1256  dphi[0] = -6 * t * (1 - t) * h_1;
1257  psi[0] = h * t * (1 - t) * (1 - t);
1258  dpsi[0] = (1 - t) * (1 - 3*t);
1259  t--;
1260  phi[1] = (1 + t) * (1 + t) * (1 - 2*t);
1261  dphi[1] = -6 * t * (1 + t) * h_1;
1262  psi[1] = h * t * (1 + t) * (1 + t);
1263  dpsi[1] = (1 + t) * (1 + 3*t);
1264  }
void ComputeMsmMgr::doneCompute ( )

Definition at line 6040 of file ComputeMsm.C.

References msmCompute, and ComputeMsm::saveResults().

Referenced by msm::PatchData::interpolation(), and msm::PatchData::interpolationC1Hermite().

6041 {
6043 }
ComputeMsm * msmCompute
Definition: ComputeMsm.C:462
void saveResults()
Definition: ComputeMsm.C:6157
void ComputeMsmMgr::doneVirialContrib ( )
inline

Definition at line 539 of file ComputeMsm.C.

References c, cntVirialContrib, gvsum, hufx, hufy, hufz, hvfx, hvfy, hvfz, hwfx, hwfy, hwfz, msm::IndexRange::ia(), msm::IndexRange::ib(), initVirialContrib(), msm::IndexRange::ja(), msm::IndexRange::jb(), msm::IndexRange::ka(), msm::IndexRange::kb(), numVirialContrib, virial, VMAX, VXX, VXY, VXZ, VYY, VYZ, and VZZ.

539  {
541  // reduce all gvsum contributions into virial tensor
542  for (int n = 0; n < VMAX; n++) { virial[n] = 0; }
543  int ia = gvsum.ia();
544  int ib = gvsum.ib();
545  int ja = gvsum.ja();
546  int jb = gvsum.jb();
547  int ka = gvsum.ka();
548  int kb = gvsum.kb();
549  for (int k = ka; k <= kb; k++) {
550  for (int j = ja; j <= jb; j++) {
551  for (int i = ia; i <= ib; i++) {
552  Float cu = Float(i);
553  Float cv = Float(j);
554  Float cw = Float(k);
555  Float c = gvsum(i,j,k);
556  Float vx = cu*hufx + cv*hvfx + cw*hwfx;
557  Float vy = cu*hufy + cv*hvfy + cw*hwfy;
558  Float vz = cu*hufz + cv*hvfz + cw*hwfz;
559  virial[VXX] -= c * vx * vx;
560  virial[VXY] -= c * vx * vy;
561  virial[VXZ] -= c * vx * vz;
562  virial[VYY] -= c * vy * vy;
563  virial[VYZ] -= c * vy * vz;
564  virial[VZZ] -= c * vz * vz;
565  }
566  }
567  }
569  }
570  }
int ka() const
Definition: MsmMap.h:438
int numVirialContrib
Definition: ComputeMsm.C:524
int ia() const
Definition: MsmMap.h:434
void initVirialContrib()
Definition: ComputeMsm.C:529
int cntVirialContrib
Definition: ComputeMsm.C:525
msm::Grid< Float > gvsum
Definition: ComputeMsm.C:523
int kb() const
Definition: MsmMap.h:439
Float virial[VMAX]
Definition: ComputeMsm.C:527
int jb() const
Definition: MsmMap.h:437
float Float
Definition: MsmMap.h:74
int ib() const
Definition: MsmMap.h:435
int ja() const
Definition: MsmMap.h:436
static void ComputeMsmMgr::gc_c1hermite_elem_accum ( C1Matrix matrix,
BigReal  _c,
Vector  rv,
BigReal  _a,
int  _split 
)
inlinestatic

Definition at line 1410 of file ComputeMsm.C.

References C1_VECTOR_SIZE, C1INDEX, D000, D001, D010, D011, D100, D101, D110, D111, C1Matrix::melem, ndsplitting(), rv, TAYLOR2, TAYLOR3, TAYLOR4, TAYLOR5, Vector::x, Vector::y, and Vector::z.

Referenced by initialize().

1411  {
1412  const BigReal a_1 = 1./_a;
1413  const BigReal a_2 = a_1 * a_1;
1414  const BigReal s = (rv * rv) * a_2;
1415  const BigReal dx = -2 * rv.x * a_2; // ds/dx
1416  const BigReal dy = -2 * rv.y * a_2; // ds/dy
1417  const BigReal dz = -2 * rv.z * a_2; // ds/dz
1418  const BigReal dd = 2 * a_2; // d^2s/dx^2 = d^2s/dy^2 = d^2s/dz^2
1419  BigReal tmp;
1420  enum { nderiv = C1_VECTOR_SIZE-1 };
1421  BigReal p[nderiv];
1422  Float *g = matrix.melem;
1423 
1424  // multiply entire matrix by this coefficient
1425  _c = _c * a_1;
1426 
1427  // compute derivatives (d/ds)^k of splitting g(s), s=r^2
1428  ndsplitting(p, s, nderiv, _split);
1429 
1430  // weight 0
1431  tmp = _c * p[0];
1432  g[C1INDEX(D000,D000)] += tmp;
1433 
1434  // weight 1
1435  tmp = _c * p[1] * dx;
1436  g[C1INDEX(D100,D000)] += tmp;
1437  g[C1INDEX(D000,D100)] -= tmp;
1438 
1439  tmp = _c * p[1] * dy;
1440  g[C1INDEX(D010,D000)] += tmp;
1441  g[C1INDEX(D000,D010)] -= tmp;
1442 
1443  tmp = _c * p[1] * dz;
1444  g[C1INDEX(D001,D000)] += tmp;
1445  g[C1INDEX(D000,D001)] -= tmp;
1446 
1447  // C1 splitting returns here
1448 
1449  // weight 2
1450  tmp = _c * p[2] * dx * dy;
1451  g[C1INDEX(D110,D000)] += tmp;
1452  g[C1INDEX(D000,D110)] += tmp;
1453  g[C1INDEX(D100,D010)] -= tmp;
1454  g[C1INDEX(D010,D100)] -= tmp;
1455 
1456  tmp = _c * p[2] * dx * dz;
1457  g[C1INDEX(D101,D000)] += tmp;
1458  g[C1INDEX(D000,D101)] += tmp;
1459  g[C1INDEX(D100,D001)] -= tmp;
1460  g[C1INDEX(D001,D100)] -= tmp;
1461 
1462  tmp = _c * p[2] * dy * dz;
1463  g[C1INDEX(D011,D000)] += tmp;
1464  g[C1INDEX(D000,D011)] += tmp;
1465  g[C1INDEX(D010,D001)] -= tmp;
1466  g[C1INDEX(D001,D010)] -= tmp;
1467 
1468  tmp = _c * (p[2] * dx*dx + p[1] * dd);
1469  g[C1INDEX(D100,D100)] -= tmp;
1470  tmp = _c * (p[2] * dy*dy + p[1] * dd);
1471  g[C1INDEX(D010,D010)] -= tmp;
1472  tmp = _c * (p[2] * dz*dz + p[1] * dd);
1473  g[C1INDEX(D001,D001)] -= tmp;
1474 
1475  // C2 splitting returns here
1476  if (_split == TAYLOR2) return;
1477 
1478  // weight 3
1479  tmp = _c * p[3] * dx * dy * dz;
1480  g[C1INDEX(D111,D000)] += tmp;
1481  g[C1INDEX(D110,D001)] -= tmp;
1482  g[C1INDEX(D101,D010)] -= tmp;
1483  g[C1INDEX(D011,D100)] -= tmp;
1484  g[C1INDEX(D100,D011)] += tmp;
1485  g[C1INDEX(D010,D101)] += tmp;
1486  g[C1INDEX(D001,D110)] += tmp;
1487  g[C1INDEX(D000,D111)] -= tmp;
1488 
1489  tmp = _c * (p[3] * dx*dx * dy + p[2] * dd * dy);
1490  g[C1INDEX(D110,D100)] -= tmp;
1491  g[C1INDEX(D100,D110)] += tmp;
1492 
1493  tmp = _c * (p[3] * dx*dx * dz + p[2] * dd * dz);
1494  g[C1INDEX(D101,D100)] -= tmp;
1495  g[C1INDEX(D100,D101)] += tmp;
1496 
1497  tmp = _c * (p[3] * dy*dy * dx + p[2] * dd * dx);
1498  g[C1INDEX(D110,D010)] -= tmp;
1499  g[C1INDEX(D010,D110)] += tmp;
1500 
1501  tmp = _c * (p[3] * dy*dy * dz + p[2] * dd * dz);
1502  g[C1INDEX(D011,D010)] -= tmp;
1503  g[C1INDEX(D010,D011)] += tmp;
1504 
1505  tmp = _c * (p[3] * dz*dz * dx + p[2] * dd * dx);
1506  g[C1INDEX(D101,D001)] -= tmp;
1507  g[C1INDEX(D001,D101)] += tmp;
1508 
1509  tmp = _c * (p[3] * dz*dz * dy + p[2] * dd * dy);
1510  g[C1INDEX(D011,D001)] -= tmp;
1511  g[C1INDEX(D001,D011)] += tmp;
1512 
1513  // C3 splitting returns here
1514  if (_split == TAYLOR3) return;
1515 
1516  // weight 4
1517  tmp = _c * (p[4] * dx*dx * dy * dz + p[3] * dd * dy * dz);
1518  g[C1INDEX(D111,D100)] -= tmp;
1519  g[C1INDEX(D100,D111)] -= tmp;
1520  g[C1INDEX(D110,D101)] += tmp;
1521  g[C1INDEX(D101,D110)] += tmp;
1522 
1523  tmp = _c * (p[4] * dy*dy * dx * dz + p[3] * dd * dx * dz);
1524  g[C1INDEX(D111,D010)] -= tmp;
1525  g[C1INDEX(D010,D111)] -= tmp;
1526  g[C1INDEX(D110,D011)] += tmp;
1527  g[C1INDEX(D011,D110)] += tmp;
1528 
1529  tmp = _c * (p[4] * dz*dz * dx * dy + p[3] * dd * dx * dy);
1530  g[C1INDEX(D111,D001)] -= tmp;
1531  g[C1INDEX(D001,D111)] -= tmp;
1532  g[C1INDEX(D101,D011)] += tmp;
1533  g[C1INDEX(D011,D101)] += tmp;
1534 
1535  tmp = _c * (p[4] * dx*dx * dy*dy + p[3] * dx*dx * dd
1536  + p[3] * dd * dy*dy + p[2] * dd * dd);
1537  g[C1INDEX(D110,D110)] += tmp;
1538  tmp = _c * (p[4] * dx*dx * dz*dz + p[3] * dx*dx * dd
1539  + p[3] * dd * dz*dz + p[2] * dd * dd);
1540  g[C1INDEX(D101,D101)] += tmp;
1541  tmp = _c * (p[4] * dy*dy * dz*dz + p[3] * dy*dy * dd
1542  + p[3] * dd * dz*dz + p[2] * dd * dd);
1543  g[C1INDEX(D011,D011)] += tmp;
1544 
1545  // C4 splitting returns here
1546  if (_split == TAYLOR4) return;
1547 
1548  // weight 5
1549  tmp = _c * (p[5] * dx*dx * dy*dy * dz + p[4] * dx*dx * dd * dz
1550  + p[4] * dd * dy*dy * dz + p[3] * dd * dd * dz);
1551  g[C1INDEX(D111,D110)] += tmp;
1552  g[C1INDEX(D110,D111)] -= tmp;
1553 
1554  tmp = _c * (p[5] * dx*dx * dz*dz * dy + p[4] * dx*dx * dd * dy
1555  + p[4] * dd * dz*dz * dy + p[3] * dd * dd * dy);
1556  g[C1INDEX(D111,D101)] += tmp;
1557  g[C1INDEX(D101,D111)] -= tmp;
1558 
1559  tmp = _c * (p[5] * dy*dy * dz*dz * dx + p[4] * dy*dy * dd * dx
1560  + p[4] * dd * dz*dz * dx + p[3] * dd * dd * dx);
1561  g[C1INDEX(D111,D011)] += tmp;
1562  g[C1INDEX(D011,D111)] -= tmp;
1563 
1564  // C5 splitting returns here
1565  if (_split == TAYLOR5) return;
1566 
1567  // weight 6
1568  tmp = _c * (p[6] * dx*dx * dy*dy * dz*dz + p[5] * dx*dx * dy*dy * dd
1569  + p[5] * dx*dx * dd * dz*dz + p[5] * dd * dy*dy * dz*dz
1570  + p[4] * dx*dx * dd * dd + p[4] * dd * dy*dy * dd
1571  + p[4] * dd * dd * dz*dz + p[3] * dd * dd * dd);
1572  g[C1INDEX(D111,D111)] -= tmp;
1573 
1574  // calculate full matrix for C6 or higher splitting
1575 
1576  } // gc_c1hermite_elem_accum()
Definition: MsmMap.h:187
Definition: MsmMap.h:187
BigReal z
Definition: Vector.h:66
Definition: MsmMap.h:187
#define C1INDEX(drj, dri)
Definition: MsmMap.h:191
Definition: MsmMap.h:187
BigReal x
Definition: Vector.h:66
Definition: MsmMap.h:187
float Float
Definition: MsmMap.h:74
static void ndsplitting(BigReal pg[], BigReal s, int n, int _split)
Definition: ComputeMsm.C:1266
Float melem[C1_MATRIX_SIZE]
Definition: MsmMap.h:111
Definition: MsmMap.h:187
BigReal y
Definition: Vector.h:66
Definition: MsmMap.h:187
Definition: MsmMap.h:187
double BigReal
Definition: common.h:112
void ComputeMsmMgr::initialize ( MsmInitMsg msg)

Definition at line 3965 of file ComputeMsm.C.

References Lattice::a(), a, Lattice::a_p(), Lattice::a_r(), approx, ASSERT, Lattice::b(), Lattice::b_p(), Lattice::b_r(), msm::Map::blockLevel, msm::Map::bsx, msm::Map::bsy, msm::Map::bsz, Lattice::c(), c, C1HERMITE, C1INDEX, Lattice::c_p(), Lattice::c_r(), cross(), CUBIC, SimParameters::cutoff, D000, D001, D010, D011, D100, D101, D110, D111, dispersion, endi(), msm::Map::foldfactor, msm::Map::gc, msm::Map::gc_c1hermite, gc_c1hermite_elem_accum(), msm::Map::gpro_c1hermite, msm::Map::gres_c1hermite, msm::Map::grespro, msm::Map::gridrange, gridScalingFactor, gridspacing, msm::Map::gvc, gvsum, gzero, hu, hufx, hufy, hufz, hv, hvfx, hvfy, hvfz, hw, hwfx, hwfy, hwfz, hxlen, hxlen_1, hylen, hylen_1, hzlen, hzlen_1, msm::IndexRange::ia(), msm::IndexRange::ib(), iINFO(), iout, msm::Map::ispx, msm::Map::ispy, msm::Map::ispz, msm::IndexRange::ja(), msm::IndexRange::jb(), msm::IndexRange::ka(), msm::IndexRange::kb(), SimParameters::lattice, lattice, Vector::length(), map, MSM_MAX_BLOCK_VOLUME, SimParameters::MSMApprox, SimParameters::MSMBlockSizeX, SimParameters::MSMBlockSizeY, SimParameters::MSMBlockSizeZ, SimParameters::MSMGridSpacing, SimParameters::MSMLevels, SimParameters::MSMPadding, SimParameters::MSMSplit, SimParameters::MSMxmax, SimParameters::MSMxmin, SimParameters::MSMymax, SimParameters::MSMymin, SimParameters::MSMzmax, SimParameters::MSMzmin, NAMD_die(), nhx, nhy, nhz, msm::IndexRange::ni(), msm::IndexRange::nj(), msm::IndexRange::nk(), nlevels, msm::IndexRange::nn(), NONIC, NONIC4, Nstencil, NUM_APPROX, NUM_SPLIT, PatchMap::numPatches(), PatchMap::Object(), Node::Object(), omega, padding, msm::Map::patchList, PhiStencil, PolyDegree, QUINTIC, QUINTIC2, msm::Grid< T >::reset(), msm::Array< T >::resize(), ru, rv, rw, s_edge, Lattice::scale(), SEPTIC, SEPTIC3, C1Matrix::set(), msm::Grid< T >::set(), msm::Grid< T >::setbounds(), setup_hgrid_1d(), setup_periodic_blocksize(), sglower, shx, shx_1, shy, shy_1, shz, shz_1, Node::simParameters, simParams, MsmInitMsg::smax, smax, MsmInitMsg::smin, smin, split, splitting(), srx_x, srx_y, srx_z, sry_x, sry_y, sry_z, srz_x, srz_y, srz_z, sx_shx, sy_shy, sz_shz, TAYLOR2, TAYLOR2_DISP, TAYLOR3, TAYLOR4, TAYLOR5, TAYLOR6, TAYLOR7, TAYLOR8, Vector::unit(), virial, VMAX, Vector::x, Vector::y, and Vector::z.

3966 {
3967 #ifdef DEBUG_MSM_VERBOSE
3968  printf("ComputeMsmMgr: initialize() PE %d\n", CkMyPe());
3969 #endif
3970 
3971  smin = msg->smin;
3972  smax = msg->smax;
3973  delete msg;
3974 
3975 #if 0
3976  printf("PE%d: initializing MSM\n", CkMyPe());
3977 #endif
3978 
3980 
3981  // get required sim params, check validity
3982  lattice = simParams->lattice;
3983 
3984  // set user-defined extent of system
3985  Vector rmin(simParams->MSMxmin, simParams->MSMymin, simParams->MSMzmin);
3986  Vector rmax(simParams->MSMxmax, simParams->MSMymax, simParams->MSMzmax);
3987  Vector sdmin = lattice.scale(rmin);
3988  Vector sdmax = lattice.scale(rmax);
3989  // swap coordinates between min and max to correct for possible rotation
3990  if (sdmin.x > sdmax.x) { double t=sdmin.x; sdmin.x=sdmax.x; sdmax.x=t; }
3991  if (sdmin.y > sdmax.y) { double t=sdmin.y; sdmin.y=sdmax.y; sdmax.y=t; }
3992  if (sdmin.z > sdmax.z) { double t=sdmin.z; sdmin.z=sdmax.z; sdmax.z=t; }
3993  // extend smin, smax by user-defined extent, where appropriate
3994  if ( ! lattice.a_p() && (sdmin.x != 0 || sdmax.x != 0)) {
3995  if (sdmin.x < smin.x) {
3996  smin.x = sdmin.x;
3997  if (CkMyPe() == 0) {
3998  iout << iINFO << "MSM extending minimum X to "
3999  << simParams->MSMxmin << " A\n" << endi;
4000  }
4001  }
4002  if (sdmax.x > smax.x) {
4003  smax.x = sdmax.x;
4004  if (CkMyPe() == 0) {
4005  iout << iINFO << "MSM extending maximum X to "
4006  << simParams->MSMxmax << " A\n" << endi;
4007  }
4008  }
4009  }
4010  if ( ! lattice.b_p() && (sdmin.y != 0 || sdmax.y != 0)) {
4011  if (sdmin.y < smin.y) {
4012  smin.y = sdmin.y;
4013  if (CkMyPe() == 0) {
4014  iout << iINFO << "MSM extending minimum Y to "
4015  << simParams->MSMymin << " A\n" << endi;
4016  }
4017  }
4018  if (sdmax.y > smax.y) {
4019  smax.y = sdmax.y;
4020  if (CkMyPe() == 0) {
4021  iout << iINFO << "MSM extending maximum Y to "
4022  << simParams->MSMymax << " A\n" << endi;
4023  }
4024  }
4025  }
4026  if ( ! lattice.c_p() && (sdmin.z != 0 || sdmax.z != 0)) {
4027  if (sdmin.z < smin.z) {
4028  smin.z = sdmin.z;
4029  if (CkMyPe() == 0) {
4030  iout << iINFO << "MSM extending minimum Z to "
4031  << simParams->MSMzmin << " A\n" << endi;
4032  }
4033  }
4034  if (sdmax.z > smax.z) {
4035  smax.z = sdmax.z;
4036  if (CkMyPe() == 0) {
4037  iout << iINFO << "MSM extending maximum Z to "
4038  << simParams->MSMzmax << " A\n" << endi;
4039  }
4040  }
4041  }
4042 
4043 #ifdef DEBUG_MSM_VERBOSE
4044  printf("smin = %g %g %g smax = %g %g %g\n",
4045  smin.x, smin.y, smin.z, smax.x, smax.y, smax.z);
4046 #endif
4047 
4048  approx = simParams->MSMApprox;
4049  if (approx < 0 || approx >= NUM_APPROX) {
4050  NAMD_die("MSM: unknown approximation requested (MSMApprox)");
4051  }
4052 
4053  split = simParams->MSMSplit;
4054  if (split < 0 || split >= NUM_SPLIT) {
4055  NAMD_die("MSM: unknown splitting requested (MSMSplit)");
4056  }
4057 
4058  if (CkMyPe() == 0) {
4059  const char *approx_str, *split_str;
4060  switch (approx) {
4061  case CUBIC: approx_str = "C1 cubic"; break;
4062  case QUINTIC: approx_str = "C1 quintic"; break;
4063  case QUINTIC2: approx_str = "C2 quintic"; break;
4064  case SEPTIC: approx_str = "C1 septic"; break;
4065  case SEPTIC3: approx_str = "C3 septic"; break;
4066  case NONIC: approx_str = "C1 nonic"; break;
4067  case NONIC4: approx_str = "C4 nonic"; break;
4068  case C1HERMITE: approx_str = "C1 Hermite"; break;
4069  default: approx_str = "unknown"; break;
4070  }
4071  switch (split) {
4072  case TAYLOR2: split_str = "C2 Taylor"; break;
4073  case TAYLOR3: split_str = "C3 Taylor"; break;
4074  case TAYLOR4: split_str = "C4 Taylor"; break;
4075  case TAYLOR5: split_str = "C5 Taylor"; break;
4076  case TAYLOR6: split_str = "C6 Taylor"; break;
4077  case TAYLOR7: split_str = "C7 Taylor"; break;
4078  case TAYLOR8: split_str = "C8 Taylor"; break;
4079  default: split_str = "unknown"; break;
4080  }
4081  iout << iINFO << "MSM using "
4082  << approx_str << " interpolation\n";
4083  iout << iINFO << "MSM using "
4084  << split_str << " splitting function\n";
4085  }
4086 
4087  a = simParams->cutoff;
4088 
4089  if (approx == C1HERMITE) {
4090  gridScalingFactor = 2;
4091  }
4092  else {
4093  gridScalingFactor = 1;
4094  }
4095 
4097  if (gridspacing <= 0) {
4098  NAMD_die("MSM: grid spacing must be greater than 0 (MSMGridSpacing)");
4099  }
4100  else if (gridspacing >= a) {
4101  NAMD_die("MSM: grid spacing must be less than cutoff (MSMGridSpacing)");
4102  }
4103 
4104  padding = gridScalingFactor * simParams->MSMPadding;
4105  if (padding < 0) {
4106  NAMD_die("MSM: padding must be non-negative (MSMPadding)");
4107  }
4108 
4109  // set maximum number of levels (default 0 adapts levels to system)
4110  nlevels = simParams->MSMLevels;
4111 
4112  // XXX dispersion unused for now
4113  dispersion = 0;
4114  if ( ! dispersion && split >= TAYLOR2_DISP) {
4115  NAMD_die("MSM: requested splitting for long-range dispersion "
4116  "(not implemented)");
4117  }
4118 
4119  // set block sizes for grid decomposition
4120  int bsx = simParams->MSMBlockSizeX / int(gridScalingFactor);
4121  int bsy = simParams->MSMBlockSizeY / int(gridScalingFactor);
4122  int bsz = simParams->MSMBlockSizeZ / int(gridScalingFactor);
4123  if (bsx <= 0 || bsy <= 0 || bsz <= 0) {
4124  NAMD_die("MSM: invalid block size requested (MSMBlockSize[XYZ])");
4125  }
4126 #ifdef MSM_FIXED_SIZE_GRID_MSG
4127  else if (bsx * bsy * bsz > MSM_MAX_BLOCK_VOLUME) {
4128  NAMD_die("MSM: requested block size (MSMBlockSize[XYZ]) too big");
4129  }
4130 #endif
4131  if (CkMyPe() == 0) {
4132  iout << iINFO << "MSM block size decomposition along X is "
4133  << bsx << " grid points\n";
4134  iout << iINFO << "MSM block size decomposition along Y is "
4135  << bsy << " grid points\n";
4136  iout << iINFO << "MSM block size decomposition along Z is "
4137  << bsz << " grid points\n";
4138  }
4139 
4140  s_edge = (PolyDegree[approx] - 1) / 2; // stencil edge size
4141  omega = 2 * PolyDegree[approx]; // smallest non-periodic grid length
4142 
4143  BigReal xlen, ylen, zlen;
4144  Vector sgmin, sgmax; // grid min and max, in scaled coordinates
4145  int ispx = lattice.a_p();
4146  int ispy = lattice.b_p();
4147  int ispz = lattice.c_p();
4148  int ispany = (ispx || ispy || ispz); // is there any periodicity?
4149 
4150  if (ispx) { // periodic along basis vector
4151  xlen = lattice.a().length();
4152  sgmax.x = 0.5;
4153  sgmin.x = -0.5;
4154  }
4155  else { // non-periodic
4156  sgmax.x = smax.x + padding; // pad the edges
4157  sgmin.x = smin.x - padding;
4158  ASSERT(gridspacing > 0);
4159  // restrict center to be on a grid point
4160  BigReal mupper = ceil(sgmax.x / (2*gridspacing));
4161  BigReal mlower = floor(sgmin.x / (2*gridspacing));
4162  sgmax.x = 2*gridspacing*mupper;
4163  sgmin.x = 2*gridspacing*mlower;
4164  xlen = sgmax.x - sgmin.x;
4165  }
4166 #ifdef DEBUG_MSM_VERBOSE
4167  printf("xlen = %g sgmin.x = %g sgmax.x = %g\n", xlen, sgmin.x, sgmax.x);
4168 #endif
4169 
4170  if (ispy) { // periodic along basis vector
4171  ylen = lattice.b().length();
4172  sgmax.y = 0.5;
4173  sgmin.y = -0.5;
4174  }
4175  else { // non-periodic
4176  sgmax.y = smax.y + padding; // pad the edges
4177  sgmin.y = smin.y - padding;
4178  ASSERT(gridspacing > 0);
4179  // restrict center to be on a grid point
4180  BigReal mupper = ceil(sgmax.y / (2*gridspacing));
4181  BigReal mlower = floor(sgmin.y / (2*gridspacing));
4182  sgmax.y = 2*gridspacing*mupper;
4183  sgmin.y = 2*gridspacing*mlower;
4184  ylen = sgmax.y - sgmin.y;
4185  }
4186 #ifdef DEBUG_MSM_VERBOSE
4187  printf("ylen = %g sgmin.y = %g sgmax.y = %g\n", ylen, sgmin.y, sgmax.y);
4188 #endif
4189 
4190  if (ispz) { // periodic along basis vector
4191  zlen = lattice.c().length();
4192  sgmax.z = 0.5;
4193  sgmin.z = -0.5;
4194  }
4195  else { // non-periodic
4196  sgmax.z = smax.z + padding; // pad the edges
4197  sgmin.z = smin.z - padding;
4198  ASSERT(gridspacing > 0);
4199  // restrict center to be on a grid point
4200  BigReal mupper = ceil(sgmax.z / (2*gridspacing));
4201  BigReal mlower = floor(sgmin.z / (2*gridspacing));
4202  sgmax.z = 2*gridspacing*mupper;
4203  sgmin.z = 2*gridspacing*mlower;
4204  zlen = sgmax.z - sgmin.z;
4205  }
4206 #ifdef DEBUG_MSM_VERBOSE
4207  printf("zlen = %g sgmin.z = %g sgmax.z = %g\n", zlen, sgmin.z, sgmax.z);
4208 #endif
4209  sglower = sgmin;
4210 
4211  int ia, ib, ja, jb, ka, kb;
4212  setup_hgrid_1d(xlen, hxlen, nhx, ia, ib, ispx);
4213  setup_hgrid_1d(ylen, hylen, nhy, ja, jb, ispy);
4214  setup_hgrid_1d(zlen, hzlen, nhz, ka, kb, ispz);
4215  hxlen_1 = 1 / hxlen;
4216  hylen_1 = 1 / hylen;
4217  hzlen_1 = 1 / hzlen;
4218  if (CkMyPe() == 0) {
4219  if (ispx || ispy || ispz) {
4220  iout << iINFO << "MSM grid spacing along X is "<< hxlen << " A\n";
4221  iout << iINFO << "MSM grid spacing along Y is "<< hylen << " A\n";
4222  iout << iINFO << "MSM grid spacing along Z is "<< hzlen << " A\n";
4223  }
4224  else {
4225  iout << iINFO << "MSM grid spacing is " << gridspacing << " A\n";
4226  }
4227  if ( ! ispx || ! ispy || ! ispz ) {
4228  iout << iINFO<<"MSM non-periodic padding is "<< padding << " A\n";
4229  }
4230  }
4231 
4232  int ni = ib - ia + 1;
4233  int nj = jb - ja + 1;
4234  int nk = kb - ka + 1;
4235  int n;
4236 
4237 #if 0
4238  // reserve temp space for factored grid transfer operation
4239  n = (nk > omega ? nk : omega); // row along z-dimension
4240  lzd.resize(n);
4241  n *= (nj > omega ? nj : omega); // plane along yz-dimensions
4242  lyzd.resize(n);
4243 #endif
4244 
4245  int lastnelems = 1;
4246  int smallestnbox = 1;
4247  int isclamped = 0;
4248  int maxlevels = nlevels; // user-defined number of levels
4249 
4250 #ifdef DEBUG_MSM_VERBOSE
4251  printf("maxlevels = %d\n", maxlevels);
4252 #endif
4253  if (nlevels <= 0) { // instead we set number of levels
4254  n = ni;
4255  if (n < nj) n = nj;
4256  if (n < nk) n = nk;
4257  for (maxlevels = 1; n > 0; n >>= 1) maxlevels++;
4258  if (ispany == 0) { // no periodicity
4259  // use rule of thumb 3/4 diameter of grid cutoff sphere
4260  int ngci = (int) ceil(3*a / hxlen) - 1;
4261  int ngcj = (int) ceil(3*a / hylen) - 1;
4262  int ngck = (int) ceil(3*a / hzlen) - 1;
4263  int omega3 = omega * omega * omega;
4264  int nhalf = (int) sqrt((double)ni * nj * nk);
4265  lastnelems = (nhalf > omega3 ? nhalf : omega3);
4266  smallestnbox = ngci * ngcj * ngck; // smaller grids don't reduce work
4267  isclamped = 1;
4268  }
4269  }
4270 #ifdef DEBUG_MSM_VERBOSE
4271  printf("maxlevels = %d\n", maxlevels);
4272 #endif
4273 
4274  // allocate space for storing grid dimensions for each level
4275  map.gridrange.resize(maxlevels);
4276 
4277  // set periodicity flags
4278  map.ispx = ispx;
4279  map.ispy = ispy;
4280  map.ispz = ispz;
4281 
4282  int level = 0;
4283  int done = 0;
4284  int alldone = 0;
4285  do {
4286  map.gridrange[level].setbounds(ia, ib, ja, jb, ka, kb);
4287 
4288  // Msm index?
4289 
4290  if (++level == nlevels) done |= 0x07; // user limit on levels
4291 
4292  if (isclamped) {
4293  int nelems = ni * nj * nk;
4294  if (nelems <= lastnelems) done |= 0x07;
4295  if (nelems <= smallestnbox) done |= 0x07;
4296  }
4297 
4298  alldone = (done == 0x07); // make sure all dimensions are done
4299 
4300  if (ispx) {
4301  ni >>= 1;
4302  ib = ni-1;
4303  if (ni & 1) done |= 0x07; // == 3 or 1
4304  else if (ni == 2) done |= 0x01; // can do one more
4305  }
4306  else {
4307  ia = -((-ia+1)/2) - s_edge;
4308  ib = (ib+1)/2 + s_edge;
4309  ni = ib - ia + 1;
4310  if (ni <= omega) done |= 0x01; // can do more restrictions
4311  }
4312 
4313  if (ispy) {
4314  nj >>= 1;
4315  jb = nj-1;
4316  if (nj & 1) done |= 0x07; // == 3 or 1
4317  else if (nj == 2) done |= 0x02; // can do one more
4318  }
4319  else {
4320  ja = -((-ja+1)/2) - s_edge;
4321  jb = (jb+1)/2 + s_edge;
4322  nj = jb - ja + 1;
4323  if (nj <= omega) done |= 0x02; // can do more restrictions
4324  }
4325 
4326  if (ispz) {
4327  nk >>= 1;
4328  kb = nk-1;
4329  if (nk & 1) done |= 0x07; // == 3 or 1
4330  else if (nk == 2) done |= 0x04; // can do one more
4331  }
4332  else {
4333  ka = -((-ka+1)/2) - s_edge;
4334  kb = (kb+1)/2 + s_edge;
4335  nk = kb - ka + 1;
4336  if (nk <= omega) done |= 0x04; // can do more restrictions
4337  }
4338  } while ( ! alldone );
4339  nlevels = level;
4340 
4341  // for periodic boundaries, don't visit top level (all 0)
4342  // toplevel visited only for all nonperiodic boundaries
4343  int toplevel = (ispany ? nlevels : nlevels - 1);
4344 
4345  // resize down to the actual number of levels (does not change alloc)
4346  map.gridrange.resize(nlevels);
4347 
4348  // print out some information about MSM
4349  if (CkMyPe() == 0) {
4350  iout << iINFO << "MSM using " << nlevels << " levels\n";
4351  for (n = 0; n < nlevels; n++) {
4352  char s[100];
4353  snprintf(s, sizeof(s), " level %d: "
4354  "[%d..%d] x [%d..%d] x [%d..%d]\n", n,
4355  map.gridrange[n].ia(), map.gridrange[n].ib(),
4356  map.gridrange[n].ja(), map.gridrange[n].jb(),
4357  map.gridrange[n].ka(), map.gridrange[n].kb());
4358  iout << iINFO << s;
4359  }
4360  iout << endi;
4361  }
4362 
4363  // find grid spacing basis vectors
4364  hu = hxlen * lattice.a().unit();
4365  hv = hylen * lattice.b().unit();
4366  hw = hzlen * lattice.c().unit();
4367  hufx = Float(hu.x);
4368  hufy = Float(hu.y);
4369  hufz = Float(hu.z);
4370  hvfx = Float(hv.x);
4371  hvfy = Float(hv.y);
4372  hvfz = Float(hv.z);
4373  hwfx = Float(hw.x);
4374  hwfy = Float(hw.y);
4375  hwfz = Float(hw.z);
4376 
4377  ru = lattice.a_r();
4378  rv = lattice.b_r();
4379  rw = lattice.c_r();
4380 
4381  // determine grid spacings in scaled space
4382  shx = ru * hu;
4383  shy = rv * hv;
4384  shz = rw * hw;
4385  shx_1 = 1 / shx;
4386  shy_1 = 1 / shy;
4387  shz_1 = 1 / shz;
4388 
4389  // row vectors to transform interpolated force back to real space
4390  // XXX Is not needed.
4391  sx_shx = shx_1 * Vector(ru.x, rv.x, rw.x);
4392  sy_shy = shy_1 * Vector(ru.y, rv.y, rw.y);
4393  sz_shz = shz_1 * Vector(ru.z, rv.z, rw.z);
4394  srx_x = Float(sx_shx.x);
4395  srx_y = Float(sx_shx.y);
4396  srx_z = Float(sx_shx.z);
4397  sry_x = Float(sy_shy.x);
4398  sry_y = Float(sy_shy.y);
4399  sry_z = Float(sy_shy.z);
4400  srz_x = Float(sz_shz.x);
4401  srz_y = Float(sz_shz.y);
4402  srz_z = Float(sz_shz.z);
4403 
4404  Vector pu = cross(hv, hw);
4405  BigReal s = (hu * pu) / (pu * pu);
4406  pu *= s; // pu is orthogonal projection of hu onto hv CROSS hw
4407 
4408  Vector pv = cross(hw, hu);
4409  s = (hv * pv) / (pv * pv);
4410  pv *= s; // pv is orthogonal projection of hv onto hw CROSS hu
4411 
4412  Vector pw = cross(hu, hv);
4413  s = (hw * pw) / (pw * pw);
4414  pw *= s; // pw is orthogonal projection of hw onto hu CROSS hv
4415 
4416  // radii for parallelepiped of weights enclosing grid cutoff sphere
4417  ni = (int) ceil(2*a / pu.length() ) - 1;
4418  nj = (int) ceil(2*a / pv.length() ) - 1;
4419  nk = (int) ceil(2*a / pw.length() ) - 1;
4420 
4421  Float scaling = 1;
4422  Float scaling_factor = 0.5f;
4423  BigReal a_1 = 1/a;
4424  BigReal a_p = a_1;
4425  if (dispersion) {
4426  a_p = a_p * a_p * a_p; // = 1/a^3
4427  a_p = a_p * a_p; // = 1/a^6
4428  scaling_factor = 1.f/64; // = 1/2^6
4429  }
4430  int i, j, k;
4431  if (approx < C1HERMITE) {
4432  // resize gc and gvc constants for number of levels
4433  map.gc.resize(nlevels);
4434  map.gvc.resize(nlevels);
4435 
4436  for (level = 0; level < toplevel; level++) {
4437  map.gc[level].setbounds(-ni, ni, -nj, nj, -nk, nk);
4438  map.gvc[level].setbounds(-ni, ni, -nj, nj, -nk, nk);
4439 
4440  for (k = -nk; k <= nk; k++) {
4441  for (j = -nj; j <= nj; j++) {
4442  for (i = -ni; i <= ni; i++) {
4443  if (level == 0) {
4444  BigReal s, t, gs=0, gt=0, g=0, dgs=0, dgt=0, dg=0;
4445  BigReal vlen = (i*hu + j*hv + k*hw).length();
4446  s = vlen * a_1;
4447  t = 0.5 * s;
4448  if (t >= 1) {
4449  g = 0;
4450  dg = 0;
4451  }
4452  else {
4453  splitting(gt, dgt, t, split);
4454  if (s >= 1) {
4455  BigReal s_1 = 1/s;
4456  if (dispersion) {
4457  gs = s_1 * s_1 * s_1; // = 1/s^3
4458  gs = gs * gs; // = 1/s^6
4459  dgs = -6 * gs * s_1;
4460  }
4461  else {
4462  gs = s_1;
4463  dgs = -gs * s_1;
4464  }
4465  }
4466  else {
4467  splitting(gs, dgs, s, split);
4468  }
4469  g = (gs - scaling_factor * gt) * a_p;
4470  BigReal c=0;
4471  if (i || j || k) {
4472  c = a_p * a_1 / vlen;
4473  }
4474  dg = 0.5 * (dgs - 0.5*scaling_factor * dgt) * c;
4475 
4476  // Msm index?
4477 
4478  }
4479  map.gc[0](i,j,k) = Float(g);
4480  map.gvc[0](i,j,k) = Float(dg);
4481  } // if level 0
4482  else {
4483  map.gc[level](i,j,k) = scaling * map.gc[0](i,j,k);
4484  map.gvc[level](i,j,k) = scaling * map.gvc[0](i,j,k);
4485  }
4486 
4487  } // for i
4488  } // for j
4489  } // for k
4490  scaling *= scaling_factor;
4491 
4492  } // for level
4493 
4494  // for summed virial factors
4495  gvsum.setbounds(-ni, ni, -nj, nj, -nk, nk);
4496  // make sure final virial sum is initialized to 0
4497  for (i = 0; i < VMAX; i++) { virial[i] = 0; }
4498 
4499  if (toplevel < nlevels) {
4500  // nonperiodic along all basis vector directions
4501  // calculate top level weights where all grid points
4502  // interact with each other
4503  ni = map.gridrange[toplevel].ni();
4504  nj = map.gridrange[toplevel].nj();
4505  nk = map.gridrange[toplevel].nk();
4506  map.gc[toplevel].setbounds(-ni, ni, -nj, nj, -nk, nk);
4507 
4508  // Msm index?
4509 
4510  for (k = -nk; k <= nk; k++) {
4511  for (j = -nj; j <= nj; j++) {
4512  for (i = -ni; i <= ni; i++) {
4513  BigReal s, gs, d;
4514  BigReal vlen = (i*hu + j*hv + k*hw).length();
4515  s = vlen * a_1;
4516  if (s >= 1) {
4517  BigReal s_1 = 1/s;
4518  if (dispersion) {
4519  gs = s_1 * s_1 * s_1; // = 1/s^3
4520  gs = gs * gs; // = 1/s^6
4521  }
4522  else {
4523  gs = s_1;
4524  }
4525  }
4526  else {
4527  splitting(gs, d, s, split);
4528  }
4529  map.gc[toplevel](i,j,k) = scaling * Float(gs * a_p);
4530  } // for i
4531  } // for j
4532  } // for k
4533  } // if toplevel
4534 
4535  // generate grespro stencil
4536  const int nstencil = Nstencil[approx];
4537  const Float *phi = PhiStencil[approx];
4538  map.grespro.set(0, nstencil, 0, nstencil, 0, nstencil);
4539  for (k = 0; k < nstencil; k++) {
4540  for (j = 0; j < nstencil; j++) {
4541  for (i = 0; i < nstencil; i++) {
4542  map.grespro(i,j,k) = phi[i] * phi[j] * phi[k];
4543  }
4544  }
4545  }
4546 
4547  } // end if approx < C1HERMITE
4548  else {
4549  // C1HERMITE
4550  // resize gc_c1hermite constants for number of levels
4551  map.gc_c1hermite.resize(nlevels);
4552  scaling = 1;
4553 
4554  for (level = 0; level < toplevel; level++) {
4555 
4556  Vector hmu = scaling * hu;
4557  Vector hmv = scaling * hv;
4558  Vector hmw = scaling * hw;
4559  BigReal am = scaling * a;
4560 
4561  map.gc_c1hermite[level].setbounds(-ni, ni, -nj, nj, -nk, nk);
4562 
4563  for (k = -nk; k <= nk; k++) {
4564  for (j = -nj; j <= nj; j++) {
4565  for (i = -ni; i <= ni; i++) {
4566  C1Matrix& m = map.gc_c1hermite[level](i,j,k);
4567  Vector rv = i*hmu + j*hmv + k*hmw;
4568  BigReal r2 = rv * rv;
4569  m.set(0);
4570  if (r2 < 4*am*am) {
4571  // accumulate D( g_{a}(0,r) ) term for this level
4572  gc_c1hermite_elem_accum(m, 1, rv, am, split);
4573  // accumulate D( -g_{2a}(0,r) ) term for this level
4574  gc_c1hermite_elem_accum(m, -1, rv, 2*am, split);
4575  } // if within cutoff
4576  }
4577  }
4578  } // end loop over gc_c1hermite elements for this level
4579  scaling *= 2; // double grid spacing and cutoff at each iteration
4580 
4581  } // end loop over levels
4582 
4583  if (toplevel < nlevels) {
4584  Vector hmu = scaling * hu;
4585  Vector hmv = scaling * hv;
4586  Vector hmw = scaling * hw;
4587  BigReal am = scaling * a;
4588 
4589  // nonperiodic along all basis vector directions
4590  // calculate top level weights where all grid points
4591  // interact with each other
4592  ni = map.gridrange[toplevel].ni();
4593  nj = map.gridrange[toplevel].nj();
4594  nk = map.gridrange[toplevel].nk();
4595  map.gc_c1hermite[toplevel].setbounds(-ni, ni, -nj, nj, -nk, nk);
4596 
4597  for (k = -nk; k <= nk; k++) {
4598  for (j = -nj; j <= nj; j++) {
4599  for (i = -ni; i <= ni; i++) {
4600  C1Matrix& m = map.gc_c1hermite[level](i,j,k);
4601  Vector rv = i*hmu + j*hmv + k*hmw;
4602  m.set(0);
4603  // accumulate D( g_{a}(0,r) ) term for this level
4604  gc_c1hermite_elem_accum(m, 1, rv, am, split);
4605  }
4606  }
4607  } // end loop over gc_c1hermite elements for top level
4608 
4609  } // end if top level
4610 
4611  // C1 Hermite restriction and prolongation stencils
4612  map.gres_c1hermite.resize(nlevels-1);
4613  map.gpro_c1hermite.resize(nlevels-1);
4614 
4615  enum {
4616  ND = 3, // stencil diameter
4617  NR = ND/2 // stencil radius
4618  };
4619 
4620  // the master basis functions PHI0 and PHI1 for the 3-point stencil
4621  // and their derivatives DPHI0 and DPHI1
4622  const double PHI0[ND] = { 0.5, 1, 0.5 };
4623  const double DPHI0[ND] = { 1.5, 0, -1.5 };
4624  const double PHI1[ND] = { -0.125, 0, 0.125 };
4625  const double DPHI1[ND] = { -0.25, 1, -0.25 };
4626 
4627  // for intermediate calculations
4628  double xphi0_base_array[ND];
4629  double dxphi0_base_array[ND];
4630  double yphi0_base_array[ND];
4631  double dyphi0_base_array[ND];
4632  double zphi0_base_array[ND];
4633  double dzphi0_base_array[ND];
4634  double xphi1_base_array[ND];
4635  double dxphi1_base_array[ND];
4636  double yphi1_base_array[ND];
4637  double dyphi1_base_array[ND];
4638  double zphi1_base_array[ND];
4639  double dzphi1_base_array[ND];
4640  // will point to center of stencil arrays
4641  double *xphi0, *dxphi0, *xphi1, *dxphi1;
4642  double *yphi0, *dyphi0, *yphi1, *dyphi1;
4643  double *zphi0, *dzphi0, *zphi1, *dzphi1;
4644 
4645  for (n = 0; n < ND; n++) {
4646  xphi0_base_array[n] = PHI0[n];
4647  dxphi0_base_array[n] = hxlen_1 * DPHI0[n]; // scale by grid spacing
4648  xphi1_base_array[n] = hxlen * PHI1[n]; // scale by grid spacing
4649  dxphi1_base_array[n] = DPHI1[n];
4650  yphi0_base_array[n] = PHI0[n];
4651  dyphi0_base_array[n] = hylen_1 * DPHI0[n]; // scale by grid spacing
4652  yphi1_base_array[n] = hylen * PHI1[n]; // scale by grid spacing
4653  dyphi1_base_array[n] = DPHI1[n];
4654  zphi0_base_array[n] = PHI0[n];
4655  dzphi0_base_array[n] = hzlen_1 * DPHI0[n]; // scale by grid spacing
4656  zphi1_base_array[n] = hzlen * PHI1[n]; // scale by grid spacing
4657  dzphi1_base_array[n] = DPHI1[n];
4658  }
4659  xphi0 = xphi0_base_array + NR; // point into center of arrays
4660  dxphi0 = dxphi0_base_array + NR;
4661  xphi1 = xphi1_base_array + NR;
4662  dxphi1 = dxphi1_base_array + NR;
4663  yphi0 = yphi0_base_array + NR;
4664  dyphi0 = dyphi0_base_array + NR;
4665  yphi1 = yphi1_base_array + NR;
4666  dyphi1 = dyphi1_base_array + NR;
4667  zphi0 = zphi0_base_array + NR;
4668  dzphi0 = dzphi0_base_array + NR;
4669  zphi1 = zphi1_base_array + NR;
4670  dzphi1 = dzphi1_base_array + NR;
4671 
4672  for (level = 0; level < nlevels-1; level++) {
4673  // allocate space for restriction and prolongation stencils
4674  map.gres_c1hermite[level].set(0, ND, 0, ND, 0, ND);
4675  map.gpro_c1hermite[level].set(0, ND, 0, ND, 0, ND);
4676 
4677  // scale up to next level grid spacing
4678  //
4679  // have to determine for each dimension whether or not
4680  // a periodic grid spacing has increased
4681  // (equivalent to if there are fewer grid points)
4682  for (n = -NR; n <= NR; n++) {
4683  if ( ! ispx ||
4684  map.gridrange[level+1].ni() < map.gridrange[level].ni() ) {
4685  dxphi0[n] *= 0.5;
4686  xphi1[n] *= 2;
4687  }
4688  if ( ! ispy ||
4689  map.gridrange[level+1].nj() < map.gridrange[level].nj() ) {
4690  dyphi0[n] *= 0.5;
4691  yphi1[n] *= 2;
4692  }
4693  if ( ! ispz ||
4694  map.gridrange[level+1].nk() < map.gridrange[level].nk() ) {
4695  dzphi0[n] *= 0.5;
4696  zphi1[n] *= 2;
4697  }
4698  }
4699 
4700  // loop over restriction stencil matrices
4701  // calculate from partial derivatives
4702  for (k = -NR; k <= NR; k++) {
4703  for (j = -NR; j <= NR; j++) {
4704  for (i = -NR; i <= NR; i++) {
4705  Float *t = map.gres_c1hermite[level](i+NR,j+NR,k+NR).melem;
4706 
4707  t[C1INDEX(D000,D000)] = xphi0[i] * yphi0[j] * zphi0[k];
4708  t[C1INDEX(D000,D100)] = dxphi0[i] * yphi0[j] * zphi0[k];
4709  t[C1INDEX(D000,D010)] = xphi0[i] * dyphi0[j] * zphi0[k];
4710  t[C1INDEX(D000,D001)] = xphi0[i] * yphi0[j] * dzphi0[k];
4711  t[C1INDEX(D000,D110)] = dxphi0[i] * dyphi0[j] * zphi0[k];
4712  t[C1INDEX(D000,D101)] = dxphi0[i] * yphi0[j] * dzphi0[k];
4713  t[C1INDEX(D000,D011)] = xphi0[i] * dyphi0[j] * dzphi0[k];
4714  t[C1INDEX(D000,D111)] = dxphi0[i] * dyphi0[j] * dzphi0[k];
4715 
4716  t[C1INDEX(D100,D000)] = xphi1[i] * yphi0[j] * zphi0[k];
4717  t[C1INDEX(D100,D100)] = dxphi1[i] * yphi0[j] * zphi0[k];
4718  t[C1INDEX(D100,D010)] = xphi1[i] * dyphi0[j] * zphi0[k];
4719  t[C1INDEX(D100,D001)] = xphi1[i] * yphi0[j] * dzphi0[k];
4720  t[C1INDEX(D100,D110)] = dxphi1[i] * dyphi0[j] * zphi0[k];
4721  t[C1INDEX(D100,D101)] = dxphi1[i] * yphi0[j] * dzphi0[k];
4722  t[C1INDEX(D100,D011)] = xphi1[i] * dyphi0[j] * dzphi0[k];
4723  t[C1INDEX(D100,D111)] = dxphi1[i] * dyphi0[j] * dzphi0[k];
4724 
4725  t[C1INDEX(D010,D000)] = xphi0[i] * yphi1[j] * zphi0[k];
4726  t[C1INDEX(D010,D100)] = dxphi0[i] * yphi1[j] * zphi0[k];
4727  t[C1INDEX(D010,D010)] = xphi0[i] * dyphi1[j] * zphi0[k];
4728  t[C1INDEX(D010,D001)] = xphi0[i] * yphi1[j] * dzphi0[k];
4729  t[C1INDEX(D010,D110)] = dxphi0[i] * dyphi1[j] * zphi0[k];
4730  t[C1INDEX(D010,D101)] = dxphi0[i] * yphi1[j] * dzphi0[k];
4731  t[C1INDEX(D010,D011)] = xphi0[i] * dyphi1[j] * dzphi0[k];
4732  t[C1INDEX(D010,D111)] = dxphi0[i] * dyphi1[j] * dzphi0[k];
4733 
4734  t[C1INDEX(D001,D000)] = xphi0[i] * yphi0[j] * zphi1[k];
4735  t[C1INDEX(D001,D100)] = dxphi0[i] * yphi0[j] * zphi1[k];
4736  t[C1INDEX(D001,D010)] = xphi0[i] * dyphi0[j] * zphi1[k];
4737  t[C1INDEX(D001,D001)] = xphi0[i] * yphi0[j] * dzphi1[k];
4738  t[C1INDEX(D001,D110)] = dxphi0[i] * dyphi0[j] * zphi1[k];
4739  t[C1INDEX(D001,D101)] = dxphi0[i] * yphi0[j] * dzphi1[k];
4740  t[C1INDEX(D001,D011)] = xphi0[i] * dyphi0[j] * dzphi1[k];
4741  t[C1INDEX(D001,D111)] = dxphi0[i] * dyphi0[j] * dzphi1[k];
4742 
4743  t[C1INDEX(D110,D000)] = xphi1[i] * yphi1[j] * zphi0[k];
4744  t[C1INDEX(D110,D100)] = dxphi1[i] * yphi1[j] * zphi0[k];
4745  t[C1INDEX(D110,D010)] = xphi1[i] * dyphi1[j] * zphi0[k];
4746  t[C1INDEX(D110,D001)] = xphi1[i] * yphi1[j] * dzphi0[k];
4747  t[C1INDEX(D110,D110)] = dxphi1[i] * dyphi1[j] * zphi0[k];
4748  t[C1INDEX(D110,D101)] = dxphi1[i] * yphi1[j] * dzphi0[k];
4749  t[C1INDEX(D110,D011)] = xphi1[i] * dyphi1[j] * dzphi0[k];
4750  t[C1INDEX(D110,D111)] = dxphi1[i] * dyphi1[j] * dzphi0[k];
4751 
4752  t[C1INDEX(D101,D000)] = xphi1[i] * yphi0[j] * zphi1[k];
4753  t[C1INDEX(D101,D100)] = dxphi1[i] * yphi0[j] * zphi1[k];
4754  t[C1INDEX(D101,D010)] = xphi1[i] * dyphi0[j] * zphi1[k];
4755  t[C1INDEX(D101,D001)] = xphi1[i] * yphi0[j] * dzphi1[k];
4756  t[C1INDEX(D101,D110)] = dxphi1[i] * dyphi0[j] * zphi1[k];
4757  t[C1INDEX(D101,D101)] = dxphi1[i] * yphi0[j] * dzphi1[k];
4758  t[C1INDEX(D101,D011)] = xphi1[i] * dyphi0[j] * dzphi1[k];
4759  t[C1INDEX(D101,D111)] = dxphi1[i] * dyphi0[j] * dzphi1[k];
4760 
4761  t[C1INDEX(D011,D000)] = xphi0[i] * yphi1[j] * zphi1[k];
4762  t[C1INDEX(D011,D100)] = dxphi0[i] * yphi1[j] * zphi1[k];
4763  t[C1INDEX(D011,D010)] = xphi0[i] * dyphi1[j] * zphi1[k];
4764  t[C1INDEX(D011,D001)] = xphi0[i] * yphi1[j] * dzphi1[k];
4765  t[C1INDEX(D011,D110)] = dxphi0[i] * dyphi1[j] * zphi1[k];
4766  t[C1INDEX(D011,D101)] = dxphi0[i] * yphi1[j] * dzphi1[k];
4767  t[C1INDEX(D011,D011)] = xphi0[i] * dyphi1[j] * dzphi1[k];
4768  t[C1INDEX(D011,D111)] = dxphi0[i] * dyphi1[j] * dzphi1[k];
4769 
4770  t[C1INDEX(D111,D000)] = xphi1[i] * yphi1[j] * zphi1[k];
4771  t[C1INDEX(D111,D100)] = dxphi1[i] * yphi1[j] * zphi1[k];
4772  t[C1INDEX(D111,D010)] = xphi1[i] * dyphi1[j] * zphi1[k];
4773  t[C1INDEX(D111,D001)] = xphi1[i] * yphi1[j] * dzphi1[k];
4774  t[C1INDEX(D111,D110)] = dxphi1[i] * dyphi1[j] * zphi1[k];
4775  t[C1INDEX(D111,D101)] = dxphi1[i] * yphi1[j] * dzphi1[k];
4776  t[C1INDEX(D111,D011)] = xphi1[i] * dyphi1[j] * dzphi1[k];
4777  t[C1INDEX(D111,D111)] = dxphi1[i] * dyphi1[j] * dzphi1[k];
4778  }
4779  }
4780  } // end loops over restriction stencil matrices
4781 
4782  // loop over prolongation stencil matrices
4783  // prolongation stencil matrices are the transpose of restriction
4784  for (k = -NR; k <= NR; k++) {
4785  for (j = -NR; j <= NR; j++) {
4786  for (i = -NR; i <= NR; i++) {
4787  Float *t = map.gres_c1hermite[level](i+NR,j+NR,k+NR).melem;
4788  Float *tt = map.gpro_c1hermite[level](i+NR,j+NR,k+NR).melem;
4789 
4790  tt[C1INDEX(D000,D000)] = t[C1INDEX(D000,D000)];
4791  tt[C1INDEX(D000,D100)] = t[C1INDEX(D100,D000)];
4792  tt[C1INDEX(D000,D010)] = t[C1INDEX(D010,D000)];
4793  tt[C1INDEX(D000,D001)] = t[C1INDEX(D001,D000)];
4794  tt[C1INDEX(D000,D110)] = t[C1INDEX(D110,D000)];
4795  tt[C1INDEX(D000,D101)] = t[C1INDEX(D101,D000)];
4796  tt[C1INDEX(D000,D011)] = t[C1INDEX(D011,D000)];
4797  tt[C1INDEX(D000,D111)] = t[C1INDEX(D111,D000)];
4798 
4799  tt[C1INDEX(D100,D000)] = t[C1INDEX(D000,D100)];
4800  tt[C1INDEX(D100,D100)] = t[C1INDEX(D100,D100)];
4801  tt[C1INDEX(D100,D010)] = t[C1INDEX(D010,D100)];
4802  tt[C1INDEX(D100,D001)] = t[C1INDEX(D001,D100)];
4803  tt[C1INDEX(D100,D110)] = t[C1INDEX(D110,D100)];
4804  tt[C1INDEX(D100,D101)] = t[C1INDEX(D101,D100)];
4805  tt[C1INDEX(D100,D011)] = t[C1INDEX(D011,D100)];
4806  tt[C1INDEX(D100,D111)] = t[C1INDEX(D111,D100)];
4807 
4808  tt[C1INDEX(D010,D000)] = t[C1INDEX(D000,D010)];
4809  tt[C1INDEX(D010,D100)] = t[C1INDEX(D100,D010)];
4810  tt[C1INDEX(D010,D010)] = t[C1INDEX(D010,D010)];
4811  tt[C1INDEX(D010,D001)] = t[C1INDEX(D001,D010)];
4812  tt[C1INDEX(D010,D110)] = t[C1INDEX(D110,D010)];
4813  tt[C1INDEX(D010,D101)] = t[C1INDEX(D101,D010)];
4814  tt[C1INDEX(D010,D011)] = t[C1INDEX(D011,D010)];
4815  tt[C1INDEX(D010,D111)] = t[C1INDEX(D111,D010)];
4816 
4817  tt[C1INDEX(D001,D000)] = t[C1INDEX(D000,D001)];
4818  tt[C1INDEX(D001,D100)] = t[C1INDEX(D100,D001)];
4819  tt[C1INDEX(D001,D010)] = t[C1INDEX(D010,D001)];
4820  tt[C1INDEX(D001,D001)] = t[C1INDEX(D001,D001)];
4821  tt[C1INDEX(D001,D110)] = t[C1INDEX(D110,D001)];
4822  tt[C1INDEX(D001,D101)] = t[C1INDEX(D101,D001)];
4823  tt[C1INDEX(D001,D011)] = t[C1INDEX(D011,D001)];
4824  tt[C1INDEX(D001,D111)] = t[C1INDEX(D111,D001)];
4825 
4826  tt[C1INDEX(D110,D000)] = t[C1INDEX(D000,D110)];
4827  tt[C1INDEX(D110,D100)] = t[C1INDEX(D100,D110)];
4828  tt[C1INDEX(D110,D010)] = t[C1INDEX(D010,D110)];
4829  tt[C1INDEX(D110,D001)] = t[C1INDEX(D001,D110)];
4830  tt[C1INDEX(D110,D110)] = t[C1INDEX(D110,D110)];
4831  tt[C1INDEX(D110,D101)] = t[C1INDEX(D101,D110)];
4832  tt[C1INDEX(D110,D011)] = t[C1INDEX(D011,D110)];
4833  tt[C1INDEX(D110,D111)] = t[C1INDEX(D111,D110)];
4834 
4835  tt[C1INDEX(D101,D000)] = t[C1INDEX(D000,D101)];
4836  tt[C1INDEX(D101,D100)] = t[C1INDEX(D100,D101)];
4837  tt[C1INDEX(D101,D010)] = t[C1INDEX(D010,D101)];
4838  tt[C1INDEX(D101,D001)] = t[C1INDEX(D001,D101)];
4839  tt[C1INDEX(D101,D110)] = t[C1INDEX(D110,D101)];
4840  tt[C1INDEX(D101,D101)] = t[C1INDEX(D101,D101)];
4841  tt[C1INDEX(D101,D011)] = t[C1INDEX(D011,D101)];
4842  tt[C1INDEX(D101,D111)] = t[C1INDEX(D111,D101)];
4843 
4844  tt[C1INDEX(D011,D000)] = t[C1INDEX(D000,D011)];
4845  tt[C1INDEX(D011,D100)] = t[C1INDEX(D100,D011)];
4846  tt[C1INDEX(D011,D010)] = t[C1INDEX(D010,D011)];
4847  tt[C1INDEX(D011,D001)] = t[C1INDEX(D001,D011)];
4848  tt[C1INDEX(D011,D110)] = t[C1INDEX(D110,D011)];
4849  tt[C1INDEX(D011,D101)] = t[C1INDEX(D101,D011)];
4850  tt[C1INDEX(D011,D011)] = t[C1INDEX(D011,D011)];
4851  tt[C1INDEX(D011,D111)] = t[C1INDEX(D111,D011)];
4852 
4853  tt[C1INDEX(D111,D000)] = t[C1INDEX(D000,D111)];
4854  tt[C1INDEX(D111,D100)] = t[C1INDEX(D100,D111)];
4855  tt[C1INDEX(D111,D010)] = t[C1INDEX(D010,D111)];
4856  tt[C1INDEX(D111,D001)] = t[C1INDEX(D001,D111)];
4857  tt[C1INDEX(D111,D110)] = t[C1INDEX(D110,D111)];
4858  tt[C1INDEX(D111,D101)] = t[C1INDEX(D101,D111)];
4859  tt[C1INDEX(D111,D011)] = t[C1INDEX(D011,D111)];
4860  tt[C1INDEX(D111,D111)] = t[C1INDEX(D111,D111)];
4861  }
4862  }
4863  } // end loops over prolongation stencil matrices
4864 
4865  } // end loop over levels
4866 
4867  } // end if C1HERMITE
4868 
4869  // calculate self energy factor for splitting
4870  BigReal gs=0, d=0;
4871  splitting(gs, d, 0, split);
4872  gzero = gs * a_p;
4873 
4874  if (CkMyPe() == 0) {
4875  iout << iINFO << "MSM finished calculating stencils\n" << endi;
4876  }
4877 
4878  // allocate map for patches
4879  PatchMap *pm = PatchMap::Object();
4880  int numpatches = pm->numPatches();
4881  map.patchList.resize(numpatches);
4882 #ifdef DEBUG_MSM_VERBOSE
4883  printf("numPatches = %d\n", numpatches);
4884 #endif
4885 
4886  // allocate map for blocks for each grid level
4887  map.blockLevel.resize(nlevels);
4888  map.bsx.resize(nlevels);
4889  map.bsy.resize(nlevels);
4890  map.bsz.resize(nlevels);
4891 #ifdef MSM_FOLD_FACTOR
4892  map.foldfactor.resize(nlevels);
4893 #endif
4894  for (level = 0; level < nlevels; level++) {
4895  msm::IndexRange& g = map.gridrange[level];
4897  int gia = g.ia();
4898  int gni = g.ni();
4899  int gja = g.ja();
4900  int gnj = g.nj();
4901  int gka = g.ka();
4902  int gnk = g.nk();
4903  map.bsx[level] = bsx;
4904  map.bsy[level] = bsy;
4905  map.bsz[level] = bsz;
4906  if (/* map.bsx[level] < gni ||
4907  map.bsy[level] < gnj ||
4908  map.bsz[level] < gnk */ 1) {
4909  // make sure that block sizes divide evenly into periodic dimensions
4910  if (ispx) setup_periodic_blocksize(map.bsx[level], gni);
4911  if (ispy) setup_periodic_blocksize(map.bsy[level], gnj);
4912  if (ispz) setup_periodic_blocksize(map.bsz[level], gnk);
4913 #ifdef MSM_DEBUG_VERBOSE
4914  if (CkMyPe() == 0) {
4915  printf("level = %d\n map.bs* = %d %d %d gn* = %d %d %d\n",
4916  level, map.bsx[level], map.bsy[level], map.bsz[level],gni,gnj,gnk);
4917  }
4918 #endif
4919  // subdivide grid into multiple blocks
4920  // == ceil(gni / bsx), etc.
4921  int bni = (gni / map.bsx[level]) + (gni % map.bsx[level] != 0);
4922  int bnj = (gnj / map.bsy[level]) + (gnj % map.bsy[level] != 0);
4923  int bnk = (gnk / map.bsz[level]) + (gnk % map.bsz[level] != 0);
4924 #ifdef MSM_FOLD_FACTOR
4925  if (/* level > 2 && */ (bni == 1 || bnj == 1 || bnk == 1)) {
4926  map.foldfactor[level].set(bsx / gni, bsy / gnj, bsz / gnk);
4927 #if 0
4928  if (CkMyPe() == 0) {
4929  printf("Setting MSM FoldFactor level %d: %d %d %d\n",
4930  level, bsx / gni, bsy / gnj, bsz / gnk);
4931  }
4932 #endif
4933  }
4934 #endif
4935  b.set(0, bni, 0, bnj, 0, bnk);
4936  for (k = 0; k < bnk; k++) {
4937  for (j = 0; j < bnj; j++) {
4938  for (i = 0; i < bni; i++) {
4939  b(i,j,k).reset();
4940  int ia = gia + i*map.bsx[level];
4941  int ib = ia + map.bsx[level] - 1;
4942  int ja = gja + j*map.bsy[level];
4943  int jb = ja + map.bsy[level] - 1;
4944  int ka = gka + k*map.bsz[level];
4945  int kb = ka + map.bsz[level] - 1;
4946  if (ib >= gia + gni) ib = gia + gni - 1;
4947  if (jb >= gja + gnj) jb = gja + gnj - 1;
4948  if (kb >= gka + gnk) kb = gka + gnk - 1;
4949  b(i,j,k).nrange.setbounds(ia, ib, ja, jb, ka, kb);
4950  }
4951  }
4952  }
4953  }
4954  /*
4955  else {
4956  // make entire grid into single block
4957  b.set(0, 1, 0, 1, 0, 1);
4958  b(0,0,0).reset();
4959  b(0,0,0).nrange.set(gia, gni, gja, gnj, gka, gnk);
4960  // set correct block dimensions
4961  map.bsx[level] = gni;
4962  map.bsy[level] = gnj;
4963  map.bsz[level] = gnk;
4964  }
4965  */
4966  }
4967  //CkExit();
4968 #ifdef DEBUG_MSM_VERBOSE
4969  printf("Done allocating map for grid levels\n");
4970  printf("Grid level decomposition:\n");
4971  for (level = 0; level < nlevels; level++) {
4973  int bia = b.ia();
4974  int bib = b.ib();
4975  int bja = b.ja();
4976  int bjb = b.jb();
4977  int bka = b.ka();
4978  int bkb = b.kb();
4979  for (k = bka; k <= bkb; k++) {
4980  for (j = bja; j <= bjb; j++) {
4981  for (i = bia; i <= bib; i++) {
4982  int ia = b(i,j,k).nrange.ia();
4983  int ib = b(i,j,k).nrange.ib();
4984  int ja = b(i,j,k).nrange.ja();
4985  int jb = b(i,j,k).nrange.jb();
4986  int ka = b(i,j,k).nrange.ka();
4987  int kb = b(i,j,k).nrange.kb();
4988  printf("level=%d id=%d %d %d [%d..%d] x [%d..%d] x [%d..%d]"
4989  " --> %d\n",
4990  level, i, j, k, ia, ib, ja, jb, ka, kb,
4991  b(i,j,k).nrange.nn());
4992  }
4993  }
4994  }
4995  }
4996 #endif
4997  if (CkMyPe() == 0) {
4998  iout << iINFO << "MSM finished creating map for grid levels\n" << endi;
4999  }
5000 
5001  initialize2();
5002 } // ComputeMsmMgr::initialize()
static Node * Object()
Definition: Node.h:86
Array< int > bsz
Definition: MsmMap.h:960
Array< Grid< C1Matrix > > gpro_c1hermite
Definition: MsmMap.h:952
Array< PatchDiagram > patchList
Definition: MsmMap.h:955
void reset(const T &t)
Definition: MsmMap.h:670
void resize(int n)
Definition: MsmMap.h:254
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
Definition: MsmMap.h:187
Vector a_r() const
Definition: Lattice.h:268
Definition: MsmMap.h:187
BigReal shz
Definition: ComputeMsm.C:613
#define MSM_MAX_BLOCK_VOLUME
Definition: MsmMap.h:37
Array< int > bsy
Definition: MsmMap.h:960
int ispx
Definition: MsmMap.h:958
int ka() const
Definition: MsmMap.h:438
static PatchMap * Object()
Definition: PatchMap.h:27
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
Vector c_r() const
Definition: Lattice.h:270
int ispy
Definition: MsmMap.h:958
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
int ia() const
Definition: MsmMap.h:434
Vector sx_shx
Definition: ComputeMsm.C:615
BigReal shy_1
Definition: ComputeMsm.C:614
ScaledPosition smin
Definition: ComputeMsm.C:593
Definition: MsmMap.h:187
Vector sy_shy
Definition: ComputeMsm.C:616
BigReal gridspacing
Definition: ComputeMsm.C:595
#define iout
Definition: InfoStream.h:87
int ispz
Definition: MsmMap.h:958
void setup_periodic_blocksize(int &bsize, int n)
Definition: ComputeMsm.C:3928
BigReal length(void) const
Definition: Vector.h:169
#define C1INDEX(drj, dri)
Definition: MsmMap.h:191
BigReal shz_1
Definition: ComputeMsm.C:614
BigReal shx
Definition: ComputeMsm.C:613
BigReal hylen_1
Definition: ComputeMsm.C:600
BigReal gridScalingFactor
Definition: ComputeMsm.C:597
Vector b_r() const
Definition: Lattice.h:269
void set(int pia, int pni, int pja, int pnj, int pka, int pnk)
Definition: MsmMap.h:608
Vector sglower
Definition: ComputeMsm.C:610
void setbounds(int pia, int pib, int pja, int pjb, int pka, int pkb)
Definition: MsmMap.h:612
Array< Grid< BlockDiagram > > blockLevel
Definition: MsmMap.h:956
BigReal hzlen_1
Definition: ComputeMsm.C:600
BigReal hxlen
Definition: ComputeMsm.C:599
Definition: MsmMap.h:187
static const Float PhiStencil[NUM_APPROX_FORMS][MAX_NSTENCIL_SKIP_ZERO]
Definition: ComputeMsm.C:660
ScaledPosition smax
Definition: ComputeMsm.C:594
Array< Grid< C1Matrix > > gres_c1hermite
Definition: MsmMap.h:951
static const int Nstencil[NUM_APPROX]
Definition: ComputeMsm.C:652
msm::Grid< Float > gvsum
Definition: ComputeMsm.C:523
#define ASSERT(E)
int kb() const
Definition: MsmMap.h:439
ScaledPosition smax
Definition: ComputeMsm.h:21
void setup_hgrid_1d(BigReal len, BigReal &hh, int &nn, int &ia, int &ib, int isperiodic)
Definition: ComputeMsm.C:3889
int nn() const
Definition: MsmMap.h:443
Float virial[VMAX]
Definition: ComputeMsm.C:527
BigReal hzlen
Definition: ComputeMsm.C:599
int jb() const
Definition: MsmMap.h:437
BigReal shy
Definition: ComputeMsm.C:613
int nj() const
Definition: MsmMap.h:441
BigReal x
Definition: Vector.h:66
Array< int > bsx
Definition: MsmMap.h:960
Array< IndexRange > gridrange
Definition: MsmMap.h:942
BigReal MSMGridSpacing
void NAMD_die(const char *err_msg)
Definition: common.C:83
Definition: MsmMap.h:187
float Float
Definition: MsmMap.h:74
Definition: MsmMap.h:187
#define simParams
Definition: Output.C:127
int numPatches(void) const
Definition: PatchMap.h:59
int nk() const
Definition: MsmMap.h:442
ScaledPosition smin
Definition: ComputeMsm.h:21
Array< FoldFactor > foldfactor
Definition: MsmMap.h:962
BigReal y
Definition: Vector.h:66
Vector b() const
Definition: Lattice.h:253
static const int PolyDegree[NUM_APPROX]
Definition: ComputeMsm.C:649
Definition: MsmMap.h:187
static void splitting(BigReal &g, BigReal &dg, BigReal r_a, int _split)
Definition: ComputeMsm.C:667
static void gc_c1hermite_elem_accum(C1Matrix &matrix, BigReal _c, Vector rv, BigReal _a, int _split)
Definition: ComputeMsm.C:1410
Array< Grid< Float > > gc
Definition: MsmMap.h:944
ScaledPosition scale(Position p) const
Definition: Lattice.h:83
int ib() const
Definition: MsmMap.h:435
Array< Grid< Float > > gvc
Definition: MsmMap.h:945
Vector sz_shz
Definition: ComputeMsm.C:617
Definition: MsmMap.h:187
msm::Map map
Definition: ComputeMsm.C:471
void set(Float r)
Definition: MsmMap.h:113
infostream & endi(infostream &s)
Definition: InfoStream.C:38
int ja() const
Definition: MsmMap.h:436
Grid< Float > grespro
Definition: MsmMap.h:946
int ni() const
Definition: MsmMap.h:440
int b_p() const
Definition: Lattice.h:274
BigReal MSMPadding
Array< Grid< C1Matrix > > gc_c1hermite
Definition: MsmMap.h:949
BigReal shx_1
Definition: ComputeMsm.C:614
int a_p() const
Definition: Lattice.h:273
BigReal padding
Definition: ComputeMsm.C:596
BigReal gzero
Definition: ComputeMsm.C:608
BigReal hylen
Definition: ComputeMsm.C:599
Vector a() const
Definition: Lattice.h:252
Vector unit(void) const
Definition: Vector.h:182
Vector c() const
Definition: Lattice.h:254
double BigReal
Definition: common.h:112
int c_p() const
Definition: Lattice.h:275
BigReal hxlen_1
Definition: ComputeMsm.C:600
Lattice lattice
Definition: ComputeMsm.C:592
void ComputeMsmMgr::initialize_create ( )

Definition at line 5794 of file ComputeMsm.C.

References approx, ASSERT, msm::Map::blockLevel, C1HERMITE, endi(), iINFO(), iout, map, msmBlock, msmC1HermiteBlock, msmC1HermiteGridCutoff, msmGridCutoff, msmProxy, msm::IndexRange::ni(), msm::IndexRange::nj(), msm::IndexRange::nk(), nlevels, MsmBlockProxyMsg::put(), MsmC1HermiteBlockProxyMsg::put(), MsmGridCutoffProxyMsg::put(), MsmC1HermiteGridCutoffProxyMsg::put(), and msm::Array< T >::resize().

5794  {
5795  int i, j, k, n, level;
5796 
5797  if (CkMyPe() == 0) {
5798 
5799  // on PE 0, create 3D chare array of MsmBlock for each level;
5800  // broadcast this array of proxies to the rest of the group
5801  if (approx == C1HERMITE) {
5803  }
5804  else {
5806  }
5807  for (level = 0; level < nlevels; level++) {
5808  int ni = map.blockLevel[level].ni();
5809  int nj = map.blockLevel[level].nj();
5810  int nk = map.blockLevel[level].nk();
5811 #ifdef MSM_NODE_MAPPING
5812  CkPrintf("Using MsmBlockMap for level %d\n", level);
5813  CProxy_MsmBlockMap blockMap = CProxy_MsmBlockMap::ckNew(level);
5814  CkArrayOptions opts(ni, nj, nk);
5815  opts.setMap(blockMap);
5816  if (approx == C1HERMITE) {
5817  msmC1HermiteBlock[level] =
5818  CProxy_MsmC1HermiteBlock::ckNew(level, opts);
5819  }
5820  else {
5821  msmBlock[level] = CProxy_MsmBlock::ckNew(level, opts);
5822  }
5823 #else
5824  if (approx == C1HERMITE) {
5825  msmC1HermiteBlock[level] =
5826  CProxy_MsmC1HermiteBlock::ckNew(level, ni, nj, nk);
5827  }
5828  else {
5829  msmBlock[level] = CProxy_MsmBlock::ckNew(level, ni, nj, nk);
5830  }
5831 #endif
5832 #ifdef DEBUG_MSM_VERBOSE
5833  printf("Create MsmBlock[%d] 3D chare array ( %d x %d x %d )\n",
5834  level, ni, nj, nk);
5835 #endif
5836  char msg[128];
5837  int nijk = ni * nj * nk;
5838  sprintf(msg, "MSM grid level %d decomposed into %d block%s"
5839  " ( %d x %d x %d )\n",
5840  level, nijk, (nijk==1 ? "" : "s"), ni, nj, nk);
5841  iout << iINFO << msg;
5842  }
5843  if (approx == C1HERMITE) {
5845  msg->put(msmC1HermiteBlock);
5846  msmProxy.recvMsmC1HermiteBlockProxy(msg); // broadcast
5847  }
5848  else {
5850  msg->put(msmBlock);
5851  msmProxy.recvMsmBlockProxy(msg); // broadcast
5852  }
5853 
5854 #ifdef MSM_GRID_CUTOFF_DECOMP
5855  // on PE 0, create 1D chare array of MsmGridCutoff
5856  // broadcast this array proxy to the rest of the group
5857 #ifdef MSM_NODE_MAPPING
5858  CkPrintf("Using MsmGridCutoffMap\n");
5859  CProxy_MsmGridCutoffMap gcutMap = CProxy_MsmGridCutoffMap::ckNew();
5860  CkArrayOptions optsgcut(numGridCutoff);
5861  optsgcut.setMap(gcutMap);
5862  if (approx == C1HERMITE) {
5863  msmC1HermiteGridCutoff = CProxy_MsmC1HermiteGridCutoff::ckNew(optsgcut);
5864  }
5865  else {
5866  msmGridCutoff = CProxy_MsmGridCutoff::ckNew(optsgcut);
5867  }
5868 #else
5869  if (approx == C1HERMITE) {
5871  CProxy_MsmC1HermiteGridCutoff::ckNew(numGridCutoff);
5872  }
5873  else {
5874  msmGridCutoff = CProxy_MsmGridCutoff::ckNew(numGridCutoff);
5875  }
5876 #endif
5877  if (approx == C1HERMITE) {
5880  gcmsg->put(&msmC1HermiteGridCutoff);
5881  msmProxy.recvMsmC1HermiteGridCutoffProxy(gcmsg);
5882  }
5883  else {
5885  gcmsg->put(&msmGridCutoff);
5886  msmProxy.recvMsmGridCutoffProxy(gcmsg);
5887  }
5888 
5889  // XXX PE 0 initializes each MsmGridCutoff
5890  // one-to-many
5891  // for M length chare array, better for each PE to initialize M/P?
5892  for (level = 0; level < nlevels; level++) { // for all levels
5894  int bni = b.ni();
5895  int bnj = b.nj();
5896  int bnk = b.nk();
5897  for (k = 0; k < bnk; k++) { // for all blocks
5898  for (j = 0; j < bnj; j++) {
5899  for (i = 0; i < bni; i++) {
5900  // source for charges
5901  msm::BlockIndex bi = msm::BlockIndex(level, msm::Ivec(i,j,k));
5902  int numSendAcross = b(i,j,k).sendAcross.len();
5903  ASSERT( numSendAcross == b(i,j,k).indexGridCutoff.len() );
5904  // for this source, loop over destinations for potentials
5905  for (n = 0; n < numSendAcross; n++) {
5906  msm::BlockSend &bs = b(i,j,k).sendAcross[n];
5907  int index = b(i,j,k).indexGridCutoff[n];
5908  MsmGridCutoffInitMsg *bsmsg = new MsmGridCutoffInitMsg(bi, bs);
5909  if (approx == C1HERMITE) {
5910  msmC1HermiteGridCutoff[index].setup(bsmsg);
5911  }
5912  else {
5913  msmGridCutoff[index].setup(bsmsg);
5914  }
5915  } // traverse sendAcross, indexGridCutoff arrays
5916 
5917  }
5918  }
5919  } // end for all blocks
5920 
5921  } // end for all levels
5922 
5923  iout << iINFO << "MSM grid cutoff calculation decomposed into "
5924  << numGridCutoff << " work objects\n";
5925 #endif
5926  iout << endi;
5927  }
5928 
5929 #ifdef DEBUG_MSM_VERBOSE
5930  printf("end of initialization\n");
5931 #endif
5932 } // ComputeMsmMgr::initialize_create()
void resize(int n)
Definition: MsmMap.h:254
std::ostream & iINFO(std::ostream &s)
Definition: InfoStream.C:107
void put(const msm::Array< CProxy_MsmC1HermiteBlock > &a)
Definition: ComputeMsm.C:182
#define iout
Definition: InfoStream.h:87
Array< Grid< BlockDiagram > > blockLevel
Definition: MsmMap.h:956
void put(const msm::Array< CProxy_MsmBlock > &a)
Definition: ComputeMsm.C:159
CProxy_ComputeMsmMgr msmProxy
Definition: ComputeMsm.C:461
#define ASSERT(E)
int nj() const
Definition: MsmMap.h:441
CProxy_MsmGridCutoff msmGridCutoff
Definition: ComputeMsm.C:467
CProxy_MsmC1HermiteGridCutoff msmC1HermiteGridCutoff
Definition: ComputeMsm.C:468
int nk() const
Definition: MsmMap.h:442
msm::Array< CProxy_MsmBlock > msmBlock
Definition: ComputeMsm.C:464
msm::Array< CProxy_MsmC1HermiteBlock > msmC1HermiteBlock
Definition: ComputeMsm.C:465
msm::Map map
Definition: ComputeMsm.C:471
void put(const CProxy_MsmC1HermiteGridCutoff *p)
Definition: ComputeMsm.C:223
infostream & endi(infostream &s)
Definition: InfoStream.C:38
int ni() const
Definition: MsmMap.h:440
void put(const CProxy_MsmGridCutoff *p)
Definition: ComputeMsm.C:205
void ComputeMsmMgr::initVirialContrib ( )
inline

Definition at line 529 of file ComputeMsm.C.

References cntVirialContrib, gvsum, and msm::Grid< T >::reset().

Referenced by doneVirialContrib().

529  {
530  gvsum.reset(0);
531  cntVirialContrib = 0;
532  }
void reset(const T &t)
Definition: MsmMap.h:670
int cntVirialContrib
Definition: ComputeMsm.C:525
msm::Grid< Float > gvsum
Definition: ComputeMsm.C:523
msm::Map& ComputeMsmMgr::mapData ( )
inline

Definition at line 447 of file ComputeMsm.C.

References map.

Referenced by ComputeMsm::doWork(), msm::PatchData::PatchData(), and ComputeMsm::saveResults().

447 { return map; }
msm::Map map
Definition: ComputeMsm.C:471
static void ComputeMsmMgr::ndsplitting ( BigReal  pg[],
BigReal  s,
int  n,
int  _split 
)
inlinestatic

Definition at line 1266 of file ComputeMsm.C.

References NAMD_die(), TAYLOR2, TAYLOR3, TAYLOR4, TAYLOR5, TAYLOR6, TAYLOR7, and TAYLOR8.

Referenced by gc_c1hermite_elem_accum().

1266  {
1267  int k = 0;
1268  if (k == n) return;
1269  if (s <= 1) {
1270  // compute derivatives of smoothed part
1271  switch (_split) {
1272  case TAYLOR2:
1273  pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8));
1274  if (k == n) break;
1275  pg[k++] = -1./2 + (s-1)*(3./4);
1276  if (k == n) break;
1277  pg[k++] = 3./4;
1278  break;
1279  case TAYLOR3:
1280  pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16)));
1281  if (k == n) break;
1282  pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16));
1283  if (k == n) break;
1284  pg[k++] = 3./4 + (s-1)*(-15./8);
1285  if (k == n) break;
1286  pg[k++] = -15./8;
1287  break;
1288  case TAYLOR4:
1289  pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
1290  + (s-1)*(35./128))));
1291  if (k == n) break;
1292  pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32)));
1293  if (k == n) break;
1294  pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32));
1295  if (k == n) break;
1296  pg[k++] = -15./8 + (s-1)*(105./16);
1297  if (k == n) break;
1298  pg[k++] = 105./16;
1299  break;
1300  case TAYLOR5:
1301  pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
1302  + (s-1)*(35./128 + (s-1)*(-63./256)))));
1303  if (k == n) break;
1304  pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
1305  + (s-1)*(-315./256))));
1306  if (k == n) break;
1307  pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64)));
1308  if (k == n) break;
1309  pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64));
1310  if (k == n) break;
1311  pg[k++] = 105./16 + (s-1)*(-945./32);
1312  if (k == n) break;
1313  pg[k++] = -945./32;
1314  break;
1315  case TAYLOR6:
1316  pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
1317  + (s-1)*(35./128 + (s-1)*(-63./256 + (s-1)*(231./1024))))));
1318  if (k == n) break;
1319  pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
1320  + (s-1)*(-315./256 + (s-1)*(693./512)))));
1321  if (k == n) break;
1322  pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64
1323  + (s-1)*(3465./512))));
1324  if (k == n) break;
1325  pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64
1326  + (s-1)*(3465./128)));
1327  if (k == n) break;
1328  pg[k++] = 105./16 + (s-1)*(-945./32 + (s-1)*(10395./128));
1329  if (k == n) break;
1330  pg[k++] = -945./32 + (s-1)*(10395./64);
1331  if (k == n) break;
1332  pg[k++] = 10395./64;
1333  break;
1334  case TAYLOR7:
1335  pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
1336  + (s-1)*(35./128 + (s-1)*(-63./256
1337  + (s-1)*(231./1024 + (s-1)*(-429./2048)))))));
1338  if (k == n) break;
1339  pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
1340  + (s-1)*(-315./256 + (s-1)*(693./512
1341  + (s-1)*(-3003./2048))))));
1342  if (k == n) break;
1343  pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64
1344  + (s-1)*(3465./512 + (s-1)*(-9009./1024)))));
1345  if (k == n) break;
1346  pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64 + (s-1)*(3465./128
1347  + (s-1)*(-45045./1024))));
1348  if (k == n) break;
1349  pg[k++] = 105./16 + (s-1)*(-945./32 + (s-1)*(10395./128
1350  + (s-1)*(-45045./256)));
1351  if (k == n) break;
1352  pg[k++] = -945./32 + (s-1)*(10395./64 + (s-1)*(-135135./256));
1353  if (k == n) break;
1354  pg[k++] = 10395./64 + (s-1)*(-135135./128);
1355  if (k == n) break;
1356  pg[k++] = -135135./128;
1357  break;
1358  case TAYLOR8:
1359  pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
1360  + (s-1)*(35./128 + (s-1)*(-63./256
1361  + (s-1)*(231./1024 + (s-1)*(-429./2048
1362  + (s-1)*(6435./32768))))))));
1363  if (k == n) break;
1364  pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
1365  + (s-1)*(-315./256 + (s-1)*(693./512
1366  + (s-1)*(-3003./2048 + (s-1)*(6435./4096)))))));
1367  if (k == n) break;
1368  pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64
1369  + (s-1)*(3465./512 + (s-1)*(-9009./1024
1370  + (s-1)*(45045./4096))))));
1371  if (k == n) break;
1372  pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64 + (s-1)*(3465./128
1373  + (s-1)*(-45045./1024 + (s-1)*(135135./2048)))));
1374  if (k == n) break;
1375  pg[k++] = 105./16 + (s-1)*(-945./32 + (s-1)*(10395./128
1376  + (s-1)*(-45045./256 + (s-1)*(675675./2048))));
1377  if (k == n) break;
1378  pg[k++] = -945./32 + (s-1)*(10395./64 + (s-1)*(-135135./256
1379  + (s-1)*(675675./512)));
1380  if (k == n) break;
1381  pg[k++] = 10395./64 + (s-1)*(-135135./128 + (s-1)*(2027025./512));
1382  if (k == n) break;
1383  pg[k++] = -135135./128 + (s-1)*(2027025./256);
1384  if (k == n) break;
1385  pg[k++] = 2027025./256;
1386  break;
1387  default:
1388  NAMD_die("Unknown MSM splitting.");
1389  }
1390  } // if (s <= 1)
1391  else { // (s > 1)
1392  // compute derivatives of s^(-1/2)
1393  const BigReal s_1 = 1./s;
1394  BigReal s_p = sqrt(s_1);
1395  BigReal p = -0.5;
1396  BigReal _c = 1;
1397  pg[k++] = _c * s_p;
1398  while (k < n) {
1399  s_p *= s_1;
1400  _c *= p;
1401  p -= 1;
1402  pg[k++] = _c * s_p;
1403  }
1404  } // else (s > 1)
1405  // higher derivatives are zero
1406  while (k < n) pg[k++] = 0;
1407  } // ndsplitting()
void NAMD_die(const char *err_msg)
Definition: common.C:83
double BigReal
Definition: common.h:112
int ComputeMsmMgr::numLevels ( ) const
inline

Definition at line 449 of file ComputeMsm.C.

References nlevels.

449 { return nlevels; }
msm::PatchPtrArray& ComputeMsmMgr::patchPtrArray ( )
inline

Definition at line 445 of file ComputeMsm.C.

References patchPtr.

Referenced by ComputeMsm::doWork(), and ComputeMsm::saveResults().

445 { return patchPtr; }
msm::PatchPtrArray patchPtr
Definition: ComputeMsm.C:476
void ComputeMsmMgr::recvMsmBlockProxy ( MsmBlockProxyMsg msg)

Definition at line 5935 of file ComputeMsm.C.

References MsmBlockProxyMsg::get(), and msmBlock.

5936 {
5937  msg->get(msmBlock);
5938  delete(msg);
5939 }
void get(msm::Array< CProxy_MsmBlock > &a)
Definition: ComputeMsm.C:168
msm::Array< CProxy_MsmBlock > msmBlock
Definition: ComputeMsm.C:464
void ComputeMsmMgr::recvMsmC1HermiteBlockProxy ( MsmC1HermiteBlockProxyMsg msg)

Definition at line 5947 of file ComputeMsm.C.

References MsmC1HermiteBlockProxyMsg::get(), and msmC1HermiteBlock.

5950 {
5951  msg->get(msmC1HermiteBlock);
5952  delete(msg);
5953 }
void get(msm::Array< CProxy_MsmC1HermiteBlock > &a)
Definition: ComputeMsm.C:192
msm::Array< CProxy_MsmC1HermiteBlock > msmC1HermiteBlock
Definition: ComputeMsm.C:465
void ComputeMsmMgr::recvMsmC1HermiteGridCutoffProxy ( MsmC1HermiteGridCutoffProxyMsg msg)

Definition at line 5955 of file ComputeMsm.C.

References MsmC1HermiteGridCutoffProxyMsg::get(), and msmC1HermiteGridCutoff.

5958 {
5959  msg->get(&msmC1HermiteGridCutoff);
5960  delete(msg);
5961 }
void get(CProxy_MsmC1HermiteGridCutoff *p)
Definition: ComputeMsm.C:229
CProxy_MsmC1HermiteGridCutoff msmC1HermiteGridCutoff
Definition: ComputeMsm.C:468
void ComputeMsmMgr::recvMsmGridCutoffProxy ( MsmGridCutoffProxyMsg msg)

Definition at line 5941 of file ComputeMsm.C.

References MsmGridCutoffProxyMsg::get(), and msmGridCutoff.

5942 {
5943  msg->get(&msmGridCutoff);
5944  delete(msg);
5945 }
CProxy_MsmGridCutoff msmGridCutoff
Definition: ComputeMsm.C:467
void get(CProxy_MsmGridCutoff *p)
Definition: ComputeMsm.C:210
void ComputeMsmMgr::setCompute ( ComputeMsm c)
inline

Definition at line 443 of file ComputeMsm.C.

References c, msmCompute, and ComputeMsm::setMgr().

443 { msmCompute = c; c->setMgr(this); } // local
ComputeMsm * msmCompute
Definition: ComputeMsm.C:462
void setMgr(ComputeMsmMgr *mgr)
Definition: ComputeMsm.h:32
void ComputeMsmMgr::setup_hgrid_1d ( BigReal  len,
BigReal hh,
int &  nn,
int &  ia,
int &  ib,
int  isperiodic 
)

Definition at line 3889 of file ComputeMsm.C.

References ASSERT, gridspacing, NAMD_die(), and s_edge.

Referenced by initialize().

3891 {
3892  ASSERT(gridspacing > 0);
3893  if (isperiodic) {
3894  const BigReal hmin = (4./5) * gridspacing;
3895  const BigReal hmax = 1.5 * hmin;
3896  hh = len;
3897  nn = 1; // start with one grid point across length
3898  while (hh >= hmax) {
3899  hh *= 0.5; // halve spacing and double grid points
3900  nn <<= 1;
3901  }
3902  if (hh < hmin) {
3903  if (nn < 4) {
3904  NAMD_die("Basis vector is too short or MSM grid spacing is too large");
3905  }
3906  hh *= (4./3); // scale hh by 4/3 and nn by 3/4
3907  nn >>= 2;
3908  nn *= 3;
3909  }
3910  // now we have: hmin <= h < hmax,
3911  // where nn is a power of 2 times no more than one power of 3
3912  ia = 0;
3913  ib = nn-1;
3914  }
3915  else {
3916  hh = gridspacing;
3917  // Instead of "nn = (int) ceil(len / hh);"
3918  // len is divisible by hh, up to roundoff error, so round to closest nn
3919  nn = (int) floor(len/hh + 0.5);
3920  ia = -s_edge;
3921  ib = nn + s_edge;
3922  }
3923 } // ComputeMsmMgr::setup_hgrid_1d()
BigReal gridspacing
Definition: ComputeMsm.C:595
#define ASSERT(E)
void NAMD_die(const char *err_msg)
Definition: common.C:83
double BigReal
Definition: common.h:112
void ComputeMsmMgr::setup_periodic_blocksize ( int &  bsize,
int  n 
)

Definition at line 3928 of file ComputeMsm.C.

References NAMD_die().

Referenced by initialize().

3929 {
3930  if (n % bsize != 0) {
3931  // n is either 2^k or 3*2^k
3932  int newbsize = 1;
3933  if (n % 3 == 0) newbsize = 3;
3934  while (newbsize < bsize && newbsize < n) newbsize *= 2;
3935  if (bsize < newbsize) newbsize /= 2;
3936  if (n % newbsize != 0) {
3937  NAMD_die("MSM grid size for periodic dimensions must be "
3938  "a power of 2 times at most one power of 3");
3939  }
3940  bsize = newbsize;
3941  }
3942  return;
3943 }
void NAMD_die(const char *err_msg)
Definition: common.C:83
static int ComputeMsmMgr::sign ( int  n)
inlinestatic

Definition at line 452 of file ComputeMsm.C.

452  {
453  return (n < 0 ? -1 : (n > 0 ? 1 : 0));
454  }
static void ComputeMsmMgr::splitting ( BigReal g,
BigReal dg,
BigReal  r_a,
int  _split 
)
inlinestatic

Definition at line 667 of file ComputeMsm.C.

References NAMD_die(), TAYLOR2, TAYLOR2_DISP, TAYLOR3, TAYLOR3_DISP, TAYLOR4, TAYLOR4_DISP, TAYLOR5, TAYLOR5_DISP, TAYLOR6, TAYLOR6_DISP, TAYLOR7, TAYLOR7_DISP, TAYLOR8, and TAYLOR8_DISP.

Referenced by initialize().

667  {
668  BigReal s = r_a * r_a; // s = (r/a)^2, assuming 0 <= s <= 1
669  switch (_split) {
670  case TAYLOR2:
671  g = 1 + (s-1)*(-1./2 + (s-1)*(3./8));
672  dg = (2*r_a)*(-1./2 + (s-1)*(3./4));
673  break;
674  case TAYLOR3:
675  g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16)));
676  dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16)));
677  break;
678  case TAYLOR4:
679  g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
680  + (s-1)*(35./128))));
681  dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
682  + (s-1)*(35./32))));
683  break;
684  case TAYLOR5:
685  g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
686  + (s-1)*(35./128 + (s-1)*(-63./256)))));
687  dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
688  + (s-1)*(35./32 + (s-1)*(-315./256)))));
689  break;
690  case TAYLOR6:
691  g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
692  + (s-1)*(35./128 + (s-1)*(-63./256
693  + (s-1)*(231./1024))))));
694  dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
695  + (s-1)*(35./32 + (s-1)*(-315./256
696  + (s-1)*(693./512))))));
697  break;
698  case TAYLOR7:
699  g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
700  + (s-1)*(35./128 + (s-1)*(-63./256
701  + (s-1)*(231./1024 + (s-1)*(-429./2048)))))));
702  dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
703  + (s-1)*(35./32 + (s-1)*(-315./256
704  + (s-1)*(693./512 + (s-1)*(-3003./2048)))))));
705  break;
706  case TAYLOR8:
707  g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
708  + (s-1)*(35./128 + (s-1)*(-63./256
709  + (s-1)*(231./1024 + (s-1)*(-429./2048
710  + (s-1)*(6435./32768))))))));
711  dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
712  + (s-1)*(35./32 + (s-1)*(-315./256
713  + (s-1)*(693./512 + (s-1)*(-3003./2048
714  + (s-1)*(6435./4096))))))));
715  break;
716  case TAYLOR2_DISP:
717  g = 1 + (s-1)*(-3 + (s-1)*(6));
718  dg = (2*r_a)*(-3 + (s-1)*(12));
719  break;
720  case TAYLOR3_DISP:
721  g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10)));
722  dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30)));
723  break;
724  case TAYLOR4_DISP:
725  g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10 + (s-1)*(15))));
726  dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30 + (s-1)*(60))));
727  break;
728  case TAYLOR5_DISP:
729  g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
730  + (s-1)*(15 + (s-1)*(-21)))));
731  dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
732  + (s-1)*(60 + (s-1)*(-105)))));
733  break;
734  case TAYLOR6_DISP:
735  g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
736  + (s-1)*(15 + (s-1)*(-21 + (s-1)*(28))))));
737  dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
738  + (s-1)*(60 + (s-1)*(-105 + (s-1)*(168))))));
739  break;
740  case TAYLOR7_DISP:
741  g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
742  + (s-1)*(15 + (s-1)*(-21 + (s-1)*(28
743  + (s-1)*(-36)))))));
744  dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
745  + (s-1)*(60 + (s-1)*(-105 + (s-1)*(168
746  + (s-1)*(-252)))))));
747  break;
748  case TAYLOR8_DISP:
749  g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
750  + (s-1)*(15 + (s-1)*(-21 + (s-1)*(28
751  + (s-1)*(-36 + (s-1)*(45))))))));
752  dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
753  + (s-1)*(60 + (s-1)*(-105 + (s-1)*(168
754  + (s-1)*(-252 + (s-1)*(360))))))));
755  break;
756  default:
757  NAMD_die("Unknown MSM splitting.");
758  } // switch
759  } // splitting()
void NAMD_die(const char *err_msg)
Definition: common.C:83
double BigReal
Definition: common.h:112
void ComputeMsmMgr::stencil_1d ( Float  phi[],
Float  t 
)
inline

Definition at line 761 of file ComputeMsm.C.

References approx, CUBIC, NAMD_die(), NONIC, NONIC4, QUINTIC, QUINTIC2, SEPTIC, and SEPTIC3.

Referenced by msm::PatchData::anterpolation().

761  {
762  switch (approx) {
763  case CUBIC:
764  phi[0] = 0.5f * (1 - t) * (2 - t) * (2 - t);
765  t--;
766  phi[1] = (1 - t) * (1 + t - 1.5f * t * t);
767  t--;
768  phi[2] = (1 + t) * (1 - t - 1.5f * t * t);
769  t--;
770  phi[3] = 0.5f * (1 + t) * (2 + t) * (2 + t);
771  break;
772  case QUINTIC:
773  phi[0] = (1.f/24) * (1-t) * (2-t) * (3-t) * (3-t) * (4-t);
774  t--;
775  phi[1] = (1-t) * (2-t) * (3-t) * ((1.f/6)
776  + t * (0.375f - (5.f/24)*t));
777  t--;
778  phi[2] = (1-t*t) * (2-t) * (0.5f + t * (0.25f - (5.f/12)*t));
779  t--;
780  phi[3] = (1-t*t) * (2+t) * (0.5f - t * (0.25f + (5.f/12)*t));
781  t--;
782  phi[4] = (1+t) * (2+t) * (3+t) * ((1.f/6)
783  - t * (0.375f + (5.f/24)*t));
784  t--;
785  phi[5] = (1.f/24) * (1+t) * (2+t) * (3+t) * (3+t) * (4+t);
786  break;
787  case QUINTIC2:
788  phi[0] = (1.f/24) * (3-t) * (3-t) * (3-t) * (t-2) * (5*t-8);
789  t--;
790  phi[1] = (-1.f/24) * (2-t) * (t-1) * (-48+t*(153+t*(-114+t*25)));
791  t--;
792  phi[2] = (1.f/12) * (1-t) * (12+t*(12+t*(-3+t*(-38+t*25))));
793  t--;
794  phi[3] = (1.f/12) * (1+t) * (12+t*(-12+t*(-3+t*(38+t*25))));
795  t--;
796  phi[4] = (-1.f/24) * (2+t) * (t+1) * (48+t*(153+t*(114+t*25)));
797  t--;
798  phi[5] = (1.f/24) * (3+t) * (3+t) * (3+t) * (t+2) * (5*t+8);
799  break;
800  case SEPTIC:
801  phi[0] = (-1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-4)*(t-5)*(t-6);
802  t--;
803  phi[1] = (1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-5)*(-6+t*(-20+7*t));
804  t--;
805  phi[2] = (-1.f/240)*(t*t-1)*(t-2)*(t-3)*(t-4)*(-10+t*(-12+7*t));
806  t--;
807  phi[3] = (1.f/144)*(t*t-1)*(t*t-4)*(t-3)*(-12+t*(-4+7*t));
808  t--;
809  phi[4] = (-1.f/144)*(t*t-1)*(t*t-4)*(t+3)*(-12+t*(4+7*t));
810  t--;
811  phi[5] = (1.f/240)*(t*t-1)*(t+2)*(t+3)*(t+4)*(-10+t*(12+7*t));
812  t--;
813  phi[6] = (-1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(-6+t*(20+7*t));
814  t--;
815  phi[7] = (1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+4)*(t+5)*(t+6);
816  break;
817  case SEPTIC3:
818  phi[0] = (3632.f/5) + t*((-7456.f/5) + t*((58786.f/45) + t*(-633
819  + t*((26383.f/144) + t*((-22807.f/720) + t*((727.f/240)
820  + t*(-89.f/720)))))));
821  t--;
822  phi[1] = -440 + t*((25949.f/20) + t*((-117131.f/72) + t*((2247.f/2)
823  + t*((-66437.f/144) + t*((81109.f/720) + t*((-727.f/48)
824  + t*(623.f/720)))))));
825  t--;
826  phi[2] = (138.f/5) + t*((-8617.f/60) + t*((12873.f/40) + t*((-791.f/2)
827  + t*((4557.f/16) + t*((-9583.f/80) + t*((2181.f/80)
828  + t*(-623.f/240)))))));
829  t--;
830  phi[3] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((2569.f/144)
831  + t*((-727.f/48) + t*(623.f/144)))));
832  t--;
833  phi[4] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((-2569.f/144)
834  + t*((-727.f/48) + t*(-623.f/144)))));
835  t--;
836  phi[5] = (138.f/5) + t*((8617.f/60) + t*((12873.f/40) + t*((791.f/2)
837  + t*((4557.f/16) + t*((9583.f/80) + t*((2181.f/80)
838  + t*(623.f/240)))))));
839  t--;
840  phi[6] = -440 + t*((-25949.f/20) + t*((-117131.f/72) + t*((-2247.f/2)
841  + t*((-66437.f/144) + t*((-81109.f/720) + t*((-727.f/48)
842  + t*(-623.f/720)))))));
843  t--;
844  phi[7] = (3632.f/5) + t*((7456.f/5) + t*((58786.f/45) + t*(633
845  + t*((26383.f/144) + t*((22807.f/720) + t*((727.f/240)
846  + t*(89.f/720)))))));
847  break;
848  case NONIC:
849  phi[0] = (-1.f/40320)*(t-8)*(t-7)*(t-6)*(t-5)*(t-5)*(t-4)*(t-3)*
850  (t-2)*(t-1);
851  t--;
852  phi[1] = (1.f/40320)*(t-7)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*
853  (-8+t*(-35+9*t));
854  t--;
855  phi[2] = (-1.f/10080)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*
856  (-14+t*(-25+9*t));
857  t--;
858  phi[3] = (1.f/1440)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*
859  (-6+t*(-5+3*t));
860  t--;
861  phi[4] = (-1.f/2880)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*
862  (-20+t*(-5+9*t));
863  t--;
864  phi[5] = (1.f/2880)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*
865  (-20+t*(5+9*t));
866  t--;
867  phi[6] = (-1.f/1440)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*
868  (-6+t*(5+3*t));
869  t--;
870  phi[7] = (1.f/10080)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*
871  (-14+t*(25+9*t));
872  t--;
873  phi[8] = (-1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*(t+7)*
874  (-8+t*(35+9*t));
875  t--;
876  phi[9] = (1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+5)*(t+6)*
877  (t+7)*(t+8);
878  break;
879  case NONIC4:
880  { // begin grouping to define local variables
881  double Tphi[10];
882  double T=t;
883  Tphi[0] = 439375./7+T*(-64188125./504+T*(231125375./2016
884  +T*(-17306975./288+T*(7761805./384+T*(-2895587./640
885  +T*(129391./192+T*(-259715./4032+T*(28909./8064
886  +T*(-3569./40320)))))))));
887  T--;
888  Tphi[1] = -56375+T*(8314091./56+T*(-49901303./288+T*(3763529./32
889  +T*(-19648027./384+T*(9469163./640+T*(-545977./192
890  +T*(156927./448+T*(-28909./1152
891  +T*(3569./4480)))))))));
892  T--;
893  Tphi[2] = 68776./7+T*(-1038011./28+T*(31157515./504+T*(-956669./16
894  +T*(3548009./96+T*(-2422263./160+T*(197255./48
895  +T*(-19959./28+T*(144545./2016
896  +T*(-3569./1120)))))))));
897  T--;
898  Tphi[3] = -154+T*(12757./12+T*(-230123./72+T*(264481./48
899  +T*(-576499./96+T*(686147./160+T*(-96277./48
900  +T*(14221./24+T*(-28909./288+T*(3569./480)))))))));
901  T--;
902  Tphi[4] = 1+T*T*(-205./144+T*T*(91./192+T*(-6181./320
903  +T*(6337./96+T*(-2745./32+T*(28909./576
904  +T*(-3569./320)))))));
905  T--;
906  Tphi[5] = 1+T*T*(-205./144+T*T*(91./192+T*(6181./320
907  +T*(6337./96+T*(2745./32+T*(28909./576
908  +T*(3569./320)))))));
909  T--;
910  Tphi[6] = -154+T*(-12757./12+T*(-230123./72+T*(-264481./48
911  +T*(-576499./96+T*(-686147./160+T*(-96277./48
912  +T*(-14221./24+T*(-28909./288+T*(-3569./480)))))))));
913  T--;
914  Tphi[7] = 68776./7+T*(1038011./28+T*(31157515./504+T*(956669./16
915  +T*(3548009./96+T*(2422263./160+T*(197255./48
916  +T*(19959./28+T*(144545./2016+T*(3569./1120)))))))));
917  T--;
918  Tphi[8] = -56375+T*(-8314091./56+T*(-49901303./288+T*(-3763529./32
919  +T*(-19648027./384+T*(-9469163./640+T*(-545977./192
920  +T*(-156927./448+T*(-28909./1152
921  +T*(-3569./4480)))))))));
922  T--;
923  Tphi[9] = 439375./7+T*(64188125./504+T*(231125375./2016
924  +T*(17306975./288+T*(7761805./384+T*(2895587./640
925  +T*(129391./192+T*(259715./4032+T*(28909./8064
926  +T*(3569./40320)))))))));
927  for (int i=0; i < 10; i++) {
928  phi[i] = Float(Tphi[i]);
929  }
930  } // end grouping to define local variables
931  break;
932  default:
933  NAMD_die("Unknown MSM approximation.");
934  } // switch
935  } // stencil_1d()
void NAMD_die(const char *err_msg)
Definition: common.C:83
float Float
Definition: MsmMap.h:74
void ComputeMsmMgr::stencil_1d_c1hermite ( Float  phi[],
Float  psi[],
Float  t,
Float  h 
)
inline

Definition at line 1244 of file ComputeMsm.C.

Referenced by msm::PatchData::anterpolationC1Hermite().

1244  {
1245  phi[0] = (1 - t) * (1 - t) * (1 + 2*t);
1246  psi[0] = h * t * (1 - t) * (1 - t);
1247  t--;
1248  phi[1] = (1 + t) * (1 + t) * (1 - 2*t);
1249  psi[1] = h * t * (1 + t) * (1 + t);
1250  }
void ComputeMsmMgr::subtractVirialContrib ( )
inline

Definition at line 536 of file ComputeMsm.C.

References numVirialContrib.

536  {
538  }
int numVirialContrib
Definition: ComputeMsm.C:524
void ComputeMsmMgr::update ( CkQdMsg *  msg)

Definition at line 5963 of file ComputeMsm.C.

References approx, C1HERMITE, msmBlock, msmC1HermiteBlock, and nlevels.

5964 {
5965 #ifdef DEBUG_MSM_VERBOSE
5966  printf("ComputeMsmMgr: update() PE %d\n", CkMyPe());
5967 #endif
5968  delete msg;
5969 
5970  // have to setup sections AFTER initialization is finished
5971  if (CkMyPe() == 0) {
5972  for (int level = 0; level < nlevels; level++) {
5973  if (approx == C1HERMITE) {
5974  msmC1HermiteBlock[level].setupSections();
5975  }
5976  else {
5977  msmBlock[level].setupSections();
5978  }
5979  }
5980  }
5981 
5982  // XXX how do update for constant pressure simulation?
5983 }
msm::Array< CProxy_MsmBlock > msmBlock
Definition: ComputeMsm.C:464
msm::Array< CProxy_MsmC1HermiteBlock > msmC1HermiteBlock
Definition: ComputeMsm.C:465

Friends And Related Function Documentation

friend struct msm::PatchData
friend

Definition at line 365 of file ComputeMsm.C.

friend class MsmBlock
friend

Definition at line 366 of file ComputeMsm.C.

friend class MsmBlockMap
friend

Definition at line 368 of file ComputeMsm.C.

friend class MsmGridCutoffMap
friend

Definition at line 369 of file ComputeMsm.C.

Member Data Documentation

BigReal ComputeMsmMgr::a

Definition at line 598 of file ComputeMsm.C.

Referenced by initialize().

int ComputeMsmMgr::approx
msm::Array<int> ComputeMsmMgr::blockAssign

Definition at line 484 of file ComputeMsm.C.

Vector ComputeMsmMgr::c

Definition at line 588 of file ComputeMsm.C.

Referenced by doneVirialContrib(), initialize(), and setCompute().

int ComputeMsmMgr::cntVirialContrib

Definition at line 525 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initVirialContrib().

int ComputeMsmMgr::dispersion

Definition at line 607 of file ComputeMsm.C.

Referenced by initialize().

msm::Array<int> ComputeMsmMgr::gcutAssign

Definition at line 485 of file ComputeMsm.C.

Referenced by MsmGridCutoffMap::MsmGridCutoffMap().

BigReal ComputeMsmMgr::gridScalingFactor

Definition at line 597 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::gridspacing

Definition at line 595 of file ComputeMsm.C.

Referenced by initialize(), and setup_hgrid_1d().

msm::Grid<Float> ComputeMsmMgr::gvsum

Definition at line 523 of file ComputeMsm.C.

Referenced by doneVirialContrib(), initialize(), and initVirialContrib().

BigReal ComputeMsmMgr::gzero
Vector ComputeMsmMgr::hu

Definition at line 601 of file ComputeMsm.C.

Referenced by initialize().

Float ComputeMsmMgr::hufx

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Float ComputeMsmMgr::hufy

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Float ComputeMsmMgr::hufz

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Vector ComputeMsmMgr::hv

Definition at line 601 of file ComputeMsm.C.

Referenced by initialize().

Float ComputeMsmMgr::hvfx

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Float ComputeMsmMgr::hvfy

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Float ComputeMsmMgr::hvfz

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Vector ComputeMsmMgr::hw

Definition at line 601 of file ComputeMsm.C.

Referenced by initialize().

Float ComputeMsmMgr::hwfx

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Float ComputeMsmMgr::hwfy

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

Float ComputeMsmMgr::hwfz

Definition at line 602 of file ComputeMsm.C.

Referenced by doneVirialContrib(), and initialize().

BigReal ComputeMsmMgr::hxlen
BigReal ComputeMsmMgr::hxlen_1
BigReal ComputeMsmMgr::hylen
BigReal ComputeMsmMgr::hylen_1
BigReal ComputeMsmMgr::hzlen
BigReal ComputeMsmMgr::hzlen_1
const int ComputeMsmMgr::IndexOffset
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},
{-1, 0, 1},
}

Definition at line 656 of file ComputeMsm.C.

Referenced by MsmBlockKernel< Vtype, Mtype >::prolongationKernel(), and MsmBlockKernel< Vtype, Mtype >::restrictionKernel().

int ComputeMsmMgr::ispu

Definition at line 590 of file ComputeMsm.C.

int ComputeMsmMgr::ispv

Definition at line 590 of file ComputeMsm.C.

int ComputeMsmMgr::ispw

Definition at line 590 of file ComputeMsm.C.

Lattice ComputeMsmMgr::lattice
msm::Map ComputeMsmMgr::map
msm::Array<CProxy_MsmBlock> ComputeMsmMgr::msmBlock
msm::Array<CProxy_MsmC1HermiteBlock> ComputeMsmMgr::msmC1HermiteBlock
CProxy_MsmC1HermiteGridCutoff ComputeMsmMgr::msmC1HermiteGridCutoff

Definition at line 468 of file ComputeMsm.C.

Referenced by initialize_create(), and recvMsmC1HermiteGridCutoffProxy().

ComputeMsm* ComputeMsmMgr::msmCompute

Definition at line 462 of file ComputeMsm.C.

Referenced by doneCompute(), and setCompute().

CProxy_MsmGridCutoff ComputeMsmMgr::msmGridCutoff

Definition at line 467 of file ComputeMsm.C.

Referenced by initialize_create(), and recvMsmGridCutoffProxy().

CProxy_ComputeMsmMgr ComputeMsmMgr::msmProxy

Definition at line 461 of file ComputeMsm.C.

Referenced by initialize_create().

int ComputeMsmMgr::nhx

Definition at line 603 of file ComputeMsm.C.

Referenced by initialize().

int ComputeMsmMgr::nhy

Definition at line 603 of file ComputeMsm.C.

Referenced by initialize().

int ComputeMsmMgr::nhz

Definition at line 603 of file ComputeMsm.C.

Referenced by initialize().

int ComputeMsmMgr::nlevels
const int ComputeMsmMgr::Nstencil
static
Initial value:
= {
5, 7, 7, 9, 9, 11, 11, 3,
}

Definition at line 652 of file ComputeMsm.C.

Referenced by initialize(), MsmBlockKernel< Vtype, Mtype >::prolongationKernel(), and MsmBlockKernel< Vtype, Mtype >::restrictionKernel().

int ComputeMsmMgr::numGridCutoff

Definition at line 469 of file ComputeMsm.C.

int ComputeMsmMgr::numVirialContrib

Definition at line 524 of file ComputeMsm.C.

Referenced by addVirialContrib(), doneVirialContrib(), and subtractVirialContrib().

int ComputeMsmMgr::omega

Definition at line 623 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::padding

Definition at line 596 of file ComputeMsm.C.

Referenced by initialize().

msm::PatchPtrArray ComputeMsmMgr::patchPtr

Definition at line 476 of file ComputeMsm.C.

Referenced by addPotential(), compute(), and patchPtrArray().

const Float ComputeMsmMgr::PhiStencil
static
Initial value:
= {
{-1.f/16, 9.f/16, 1, 9.f/16, -1.f/16},
{3.f/256, -25.f/256, 75.f/128, 1, 75.f/128, -25.f/256, 3.f/256},
{3.f/256, -25.f/256, 75.f/128, 1, 75.f/128, -25.f/256, 3.f/256},
{ -5.f/2048, 49.f/2048, -245.f/2048, 1225.f/2048, 1, 1225.f/2048,
-245.f/2048, 49.f/2048, -5.f/2048 },
{ -5.f/2048, 49.f/2048, -245.f/2048, 1225.f/2048, 1, 1225.f/2048,
-245.f/2048, 49.f/2048, -5.f/2048 },
{ 35.f/65536, -405.f/65536, 567.f/16384, -2205.f/16384,
19845.f/32768, 1, 19845.f/32768, -2205.f/16384, 567.f/16384,
-405.f/65536, 35.f/65536 },
{ 35.f/65536, -405.f/65536, 567.f/16384, -2205.f/16384,
19845.f/32768, 1, 19845.f/32768, -2205.f/16384, 567.f/16384,
-405.f/65536, 35.f/65536 },
}

Definition at line 660 of file ComputeMsm.C.

Referenced by initialize().

const int ComputeMsmMgr::PolyDegree
static
Initial value:
= {
3, 5, 5, 7, 7, 9, 9, 1,
}

Definition at line 649 of file ComputeMsm.C.

Referenced by msm::PatchData::anterpolation(), initialize(), and msm::PatchData::interpolation().

Vector ComputeMsmMgr::ru

Definition at line 589 of file ComputeMsm.C.

Referenced by initialize().

Vector ComputeMsmMgr::rv

Definition at line 589 of file ComputeMsm.C.

Referenced by gc_c1hermite_elem_accum(), and initialize().

Vector ComputeMsmMgr::rw

Definition at line 589 of file ComputeMsm.C.

Referenced by initialize().

int ComputeMsmMgr::s_edge
Vector ComputeMsmMgr::sglower
BigReal ComputeMsmMgr::shx

Definition at line 613 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::shx_1
BigReal ComputeMsmMgr::shy

Definition at line 613 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::shy_1
BigReal ComputeMsmMgr::shz

Definition at line 613 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::shz_1
ScaledPosition ComputeMsmMgr::smax

Definition at line 594 of file ComputeMsm.C.

Referenced by initialize().

ScaledPosition ComputeMsmMgr::smin

Definition at line 593 of file ComputeMsm.C.

Referenced by initialize().

int ComputeMsmMgr::split

Definition at line 605 of file ComputeMsm.C.

Referenced by initialize().

Float ComputeMsmMgr::srx_x
Float ComputeMsmMgr::srx_y
Float ComputeMsmMgr::srx_z
Float ComputeMsmMgr::sry_x
Float ComputeMsmMgr::sry_y
Float ComputeMsmMgr::sry_z
Float ComputeMsmMgr::srz_x
Float ComputeMsmMgr::srz_y
Float ComputeMsmMgr::srz_z
msm::Grid<Float> ComputeMsmMgr::subgrid
msm::Grid<C1Vector> ComputeMsmMgr::subgrid_c1hermite

Definition at line 481 of file ComputeMsm.C.

Referenced by addPotential().

Vector ComputeMsmMgr::sx_shx

Definition at line 615 of file ComputeMsm.C.

Referenced by initialize().

Vector ComputeMsmMgr::sy_shy

Definition at line 616 of file ComputeMsm.C.

Referenced by initialize().

Vector ComputeMsmMgr::sz_shz

Definition at line 617 of file ComputeMsm.C.

Referenced by initialize().

Vector ComputeMsmMgr::u

Definition at line 588 of file ComputeMsm.C.

Vector ComputeMsmMgr::v

Definition at line 588 of file ComputeMsm.C.

Float ComputeMsmMgr::virial[VMAX]

Definition at line 527 of file ComputeMsm.C.

Referenced by doneVirialContrib(), initialize(), and ComputeMsm::saveResults().

Vector ComputeMsmMgr::w

Definition at line 588 of file ComputeMsm.C.


The documentation for this class was generated from the following file: