ComputeMsmMgr Class Reference

List of all members.

Public Types

 VXX = 0
 VXY
 VXZ
 VYY
 VYZ
 VZZ
 VMAX
 CUBIC = 0
 QUINTIC
 QUINTIC2
 SEPTIC
 SEPTIC3
 NONIC
 NONIC4
 C1HERMITE
 NUM_APPROX
 TAYLOR2 = 0
 TAYLOR3
 TAYLOR4
 TAYLOR5
 TAYLOR6
 TAYLOR7
 TAYLOR8
 TAYLOR2_DISP
 TAYLOR3_DISP
 TAYLOR4_DISP
 TAYLOR5_DISP
 TAYLOR6_DISP
 TAYLOR7_DISP
 TAYLOR8_DISP
 NUM_SPLIT
 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
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.

00526 { VXX=0, VXY, VXZ, VYY, VYZ, VZZ, VMAX };

anonymous enum

Enumerator:
MAX_POLY_DEGREE 
MAX_NSTENCIL_SIZE 
MAX_NSTENCIL_SKIP_ZERO 
NUM_APPROX_FORMS 

Definition at line 633 of file ComputeMsm.C.

00633        {
00634     // Approximation formulas with up to degree 9 polynomials.
00635     MAX_POLY_DEGREE = 9,
00636 
00637     // Max stencil length for polynomial approximation.
00638     MAX_NSTENCIL_SIZE = (2*MAX_POLY_DEGREE + 1),
00639 
00640     // Max stencil length when skipping zeros
00641     // (almost half entries are zero for interpolating polynomials).
00642     MAX_NSTENCIL_SKIP_ZERO = (MAX_POLY_DEGREE + 2),
00643 
00644     // Number of scalar approximation formulaes
00645     NUM_APPROX_FORMS = (NONIC4 - CUBIC) + 1
00646   };

enum ComputeMsmMgr::Approx

Enumerator:
CUBIC 
QUINTIC 
QUINTIC2 
SEPTIC 
SEPTIC3 
NONIC 
NONIC4 
C1HERMITE 
NUM_APPROX 

Definition at line 625 of file ComputeMsm.C.

enum ComputeMsmMgr::Split

Enumerator:
TAYLOR2 
TAYLOR3 
TAYLOR4 
TAYLOR5 
TAYLOR6 
TAYLOR7 
TAYLOR8 
TAYLOR2_DISP 
TAYLOR3_DISP 
TAYLOR4_DISP 
TAYLOR5_DISP 
TAYLOR6_DISP 
TAYLOR7_DISP 
TAYLOR8_DISP 
NUM_SPLIT 

Definition at line 628 of file ComputeMsm.C.


Constructor & Destructor Documentation

ComputeMsmMgr::ComputeMsmMgr (  ) 

Definition at line 3850 of file ComputeMsm.C.

03850                              :
03851   msmProxy(thisgroup), msmCompute(0)
03852 {
03853 #ifdef DEBUG_MSM_VERBOSE
03854   printf("ComputeMsmMgr:  (constructor) PE %d\n", CkMyPe());
03855 #endif
03856   CkpvAccess(BOCclass_group).computeMsmMgr = thisgroup;
03857 
03858 #ifdef MSM_TIMING
03859   if (CkMyPe() == 0) {
03860     msmTimer = CProxy_MsmTimer::ckNew();
03861   }
03862   initTiming();
03863 #endif
03864 #ifdef MSM_PROFILING
03865   if (CkMyPe() == 0) {
03866     msmProfiler = CProxy_MsmProfiler::ckNew();
03867   }
03868   initProfiling();
03869 #endif
03870 }

ComputeMsmMgr::~ComputeMsmMgr (  ) 

Definition at line 3872 of file ComputeMsm.C.

03873 {
03874 #ifdef DEBUG_MSM_VERBOSE
03875   printf("ComputeMsmMgr:  (destructor) PE %d\n", CkMyPe());
03876 #endif
03877   // free memory?
03878 }


Member Function Documentation

void ComputeMsmMgr::addPotential ( GridMsg  ) 

Definition at line 6014 of file ComputeMsm.C.

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

Referenced by MsmC1HermiteGridCutoff::compute(), and MsmGridCutoff::compute().

06015 {
06016   int pid;  // receive patch ID
06017   int pseq;
06018   if (approx == C1HERMITE) {
06019     gm->get(subgrid_c1hermite, pid, pseq);
06020   }
06021   else {
06022     gm->get(subgrid, pid, pseq);
06023   }
06024   delete gm;
06025   if (patchPtr[pid] == NULL) {
06026     char msg[100];
06027     snprintf(msg, sizeof(msg), "Expecting patch %d to exist on PE %d",
06028         pid, CkMyPe());
06029     NAMD_die(msg);
06030   }
06031   if (approx == C1HERMITE) {
06032     patchPtr[pid]->addPotentialC1Hermite(subgrid_c1hermite);
06033   }
06034   else {
06035     patchPtr[pid]->addPotential(subgrid);
06036   }
06037 }

void ComputeMsmMgr::addVirialContrib (  )  [inline]

Definition at line 533 of file ComputeMsm.C.

References numVirialContrib.

00533                           {
00534     numVirialContrib++;
00535   }

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.

00487                                                      {
00488     int n = 0;
00489     for (int l = 0;  l < level;  l++) {
00490       n += map.blockLevel[l].nn();
00491     }
00492     return (n + map.blockLevel[level].flatindex(i,j,k));
00493   }

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::Map::gc, msm::Map::gc_c1hermite, msm::Ivec::i, msm::Ivec::j, msm::Ivec::k, map, msm::BlockDiagram::nrange, and msm::BlockDiagram::nrangeCutoff.

00494                                                 {
00495     // XXX ratio of work for MsmBlock to MsmGridCutoff?
00496     const float scalingFactor = 3;
00497     const int volumeFullBlock = map.bsx[0] * map.bsy[0] * map.bsz[0];
00498     msm::Ivec gn;
00499     if (approx == C1HERMITE) {
00500       gn = map.gc_c1hermite[0].extent();
00501     }
00502     else {
00503       gn = map.gc[0].extent();
00504     }
00505     const int volumeFullCutoff = (map.bsx[0] + gn.i - 1) *
00506       (map.bsy[0] + gn.j - 1) * (map.bsz[0] + gn.k - 1);
00507     msm::Ivec n = b.nrange.extent();
00508     int volumeBlock = n.i * n.j * n.k;
00509     msm::Ivec nc = b.nrangeCutoff.extent();
00510     int volumeCutoff = nc.i * nc.j * nc.k;
00511     return( scalingFactor * (float(volumeBlock) / volumeFullBlock) *
00512         (float(volumeCutoff) / volumeFullCutoff) );
00513   }

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::Ivec::i, msm::Ivec::j, msm::Ivec::k, map, and msm::BlockSend::nrange_wrap.

00514                                              {
00515     const int volumeFullBlock = map.bsx[0] * map.bsy[0] * map.bsz[0];
00516     msm::Ivec n = bs.nrange_wrap.extent();;
00517     int volumeBlock = n.i * n.j * n.k;
00518     return( float(volumeBlock) / volumeFullBlock );
00519   }

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.

05987 {
05988 #ifdef DEBUG_MSM_VERBOSE
05989   printf("ComputeMsmMgr:  compute() PE=%d\n", CkMyPe());
05990 #endif
05991 
05992   int n; 
05993   for (n = 0;  n < patchIDList.len();  n++) {
05994     int patchID = patchIDList[n];
05995     if (patchPtr[patchID] == NULL) {
05996       char msg[100];
05997       snprintf(msg, sizeof(msg),
05998           "Expected MSM data for patch %d does not exist on PE %d",
05999           patchID, CkMyPe());
06000       NAMD_die(msg);
06001     }
06002     if (approx == C1HERMITE) {
06003       patchPtr[patchID]->anterpolationC1Hermite();
06004     }
06005     else {
06006       patchPtr[patchID]->anterpolation();
06007     }
06008     // all else should follow from here
06009   }
06010   return;
06011 }

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, f, NAMD_die(), NONIC, NONIC4, QUINTIC, QUINTIC2, SEPTIC, and SEPTIC3.

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

00937                                                                    {
00938     switch (approx) {
00939       case CUBIC:
00940         phi[0] = 0.5f * (1 - t) * (2 - t) * (2 - t);
00941         dphi[0] = (1.5f * t - 2) * (2 - t) * h_1;
00942         t--;
00943         phi[1] = (1 - t) * (1 + t - 1.5f * t * t);
00944         dphi[1] = (-5 + 4.5f * t) * t * h_1;
00945         t--;
00946         phi[2] = (1 + t) * (1 - t - 1.5f * t * t);
00947         dphi[2] = (-5 - 4.5f * t) * t * h_1;
00948         t--;
00949         phi[3] = 0.5f * (1 + t) * (2 + t) * (2 + t);
00950         dphi[3] = (1.5f * t + 2) * (2 + t) * h_1;
00951         break;
00952       case QUINTIC:
00953         phi[0] = (1.f/24) * (1-t) * (2-t) * (3-t) * (3-t) * (4-t);
00954         dphi[0] = ((-1.f/24) * ((3-t) * (3-t) * (14 + t * (-14 + 3*t))
00955               + 2 * (1-t) * (2-t) * (3-t) * (4-t))) * h_1;
00956         t--;
00957         phi[1] = (1-t) * (2-t) * (3-t) * ((1.f/6)
00958             + t * (0.375f - (5.f/24)*t));
00959         dphi[1] = (-((1.f/6) + t * (0.375f - (5.f/24)*t)) *
00960             (11 + t * (-12 + 3*t)) + (1-t) * (2-t) * (3-t) *
00961             (0.375f - (5.f/12)*t)) * h_1;
00962         t--;
00963         phi[2] = (1-t*t) * (2-t) * (0.5f + t * (0.25f - (5.f/12)*t));
00964         dphi[2] = (-(0.5f + t * (0.25f - (5.f/12)*t)) * (1 + t * (4 - 3*t))
00965             + (1-t*t) * (2-t) * (0.25f - (5.f/6)*t)) * h_1;
00966         t--;
00967         phi[3] = (1-t*t) * (2+t) * (0.5f - t * (0.25f + (5.f/12)*t));
00968         dphi[3] = ((0.5f + t * (-0.25f - (5.f/12)*t)) * (1 + t * (-4 - 3*t))
00969             - (1-t*t) * (2+t) * (0.25f + (5.f/6)*t)) * h_1;
00970         t--;
00971         phi[4] = (1+t) * (2+t) * (3+t) * ((1.f/6)
00972             - t * (0.375f + (5.f/24)*t));
00973         dphi[4] = (((1.f/6) + t * (-0.375f - (5.f/24)*t)) *
00974             (11 + t * (12 + 3*t)) - (1+t) * (2+t) * (3+t) *
00975             (0.375f + (5.f/12)*t)) * h_1;
00976         t--;
00977         phi[5] = (1.f/24) * (1+t) * (2+t) * (3+t) * (3+t) * (4+t);
00978         dphi[5] = ((1.f/24) * ((3+t) * (3+t) * (14 + t * (14 + 3*t))
00979               + 2 * (1+t) * (2+t) * (3+t) * (4+t))) * h_1;
00980         break;
00981       case QUINTIC2:
00982         phi[0] = (1.f/24) * (3-t) * (3-t) * (3-t) * (t-2) * (5*t-8);
00983         dphi[0] = ((1.f/24) * (3-t) * (3-t) * ((3-t)*(5*t-8)
00984               - 3*(t-2)*(5*t-8) + 5*(t-2)*(3-t))) * h_1;
00985         t--;
00986         phi[1] = (-1.f/24) * (2-t) * (t-1) * (-48+t*(153+t*(-114+t*25)));
00987         dphi[1] = ((-1.f/24) * ((2-t)*(-48+t*(153+t*(-114+t*25)))
00988               - (t-1)* (-48+t*(153+t*(-114+t*25)))
00989               + (2-t)*(t-1)*(153+t*(-228+t*75)))) * h_1;
00990         t--;
00991         phi[2] = (1.f/12) * (1-t) * (12+t*(12+t*(-3+t*(-38+t*25))));
00992         dphi[2] = ((1.f/12) * (-(12+t*(12+t*(-3+t*(-38+t*25))))
00993               + (1-t)*(12+t*(-6+t*(-114+t*100))))) * h_1;
00994         t--;
00995         phi[3] = (1.f/12) * (1+t) * (12+t*(-12+t*(-3+t*(38+t*25))));
00996         dphi[3] = ((1.f/12) * ((12+t*(-12+t*(-3+t*(38+t*25))))
00997               + (1+t)*(-12+t*(-6+t*(114+t*100))))) * h_1;
00998         t--;
00999         phi[4] = (-1.f/24) * (2+t) * (t+1) * (48+t*(153+t*(114+t*25)));
01000         dphi[4] = ((-1.f/24) * ((2+t)*(48+t*(153+t*(114+t*25)))
01001               + (t+1)* (48+t*(153+t*(114+t*25)))
01002               + (2+t)*(t+1)*(153+t*(228+t*75)))) * h_1;
01003         t--;
01004         phi[5] = (1.f/24) * (3+t) * (3+t) * (3+t) * (t+2) * (5*t+8);
01005         dphi[5] = ((1.f/24) * (3+t) * (3+t) * ((3+t)*(5*t+8)
01006               + 3*(t+2)*(5*t+8) + 5*(t+2)*(3+t))) * h_1;
01007         break;
01008       case SEPTIC:
01009         phi[0] = (-1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-4)*(t-5)*(t-6);
01010         dphi[0] = (-1.f/720)*(t-4)*(-1944+t*(3644+t*(-2512+t*(807
01011                   +t*(-122+t*7))))) * h_1;
01012         t--;
01013         phi[1] = (1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-5)*(-6+t*(-20+7*t));
01014         dphi[1] = (1.f/720)*(756+t*(-9940+t*(17724+t*(-12740+t*(4445
01015                     +t*(-750+t*49)))))) * h_1;
01016         t--;
01017         phi[2] = (-1.f/240)*(t*t-1)*(t-2)*(t-3)*(t-4)*(-10+t*(-12+7*t));
01018         dphi[2] = (-1.f/240)*(-28+t*(1260+t*(-756+t*(-1260+t*(1365
01019                     +t*(-450+t*49)))))) * h_1;
01020         t--;
01021         phi[3] = (1.f/144)*(t*t-1)*(t*t-4)*(t-3)*(-12+t*(-4+7*t));
01022         dphi[3] = (1.f/144)*t*(-560+t*(84+t*(644+t*(-175
01023                   +t*(-150+t*49))))) * h_1;
01024         t--;
01025         phi[4] = (-1.f/144)*(t*t-1)*(t*t-4)*(t+3)*(-12+t*(4+7*t));
01026         dphi[4] = (-1.f/144)*t*(560+t*(84+t*(-644+t*(-175
01027                   +t*(150+t*49))))) * h_1;
01028         t--;
01029         phi[5] = (1.f/240)*(t*t-1)*(t+2)*(t+3)*(t+4)*(-10+t*(12+7*t));
01030         dphi[5] = (1.f/240)*(-28+t*(-1260+t*(-756+t*(1260+t*(1365
01031                     +t*(450+t*49)))))) * h_1;
01032         t--;
01033         phi[6] = (-1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(-6+t*(20+7*t));
01034         dphi[6] = (-1.f/720)*(756+t*(9940+t*(17724+t*(12740+t*(4445
01035                     +t*(750+t*49)))))) * h_1;
01036         t--;
01037         phi[7] = (1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+4)*(t+5)*(t+6);
01038         dphi[7] = (1.f/720)*(t+4)*(1944+t*(3644+t*(2512+t*(807
01039                   +t*(122+t*7))))) * h_1;
01040         break;
01041       case SEPTIC3:
01042         phi[0] = (3632.f/5) + t*((-7456.f/5) + t*((58786.f/45) + t*(-633
01043                 + t*((26383.f/144) + t*((-22807.f/720) + t*((727.f/240)
01044                       + t*(-89.f/720)))))));
01045         dphi[0] = ((-7456.f/5) + t*((117572.f/45) + t*(-1899
01046                 + t*((26383.f/36) + t*((-22807.f/144) + t*((727.f/40)
01047                       + t*(-623.f/720))))))) * h_1;
01048         t--;
01049         phi[1] = -440 + t*((25949.f/20) + t*((-117131.f/72) + t*((2247.f/2)
01050                 + t*((-66437.f/144) + t*((81109.f/720) + t*((-727.f/48)
01051                       + t*(623.f/720)))))));
01052         dphi[1] = ((25949.f/20) + t*((-117131.f/36) + t*((6741.f/2)
01053                 + t*((-66437.f/36) + t*((81109.f/144) + t*((-727.f/8)
01054                       + t*(4361.f/720))))))) * h_1;
01055         t--;
01056         phi[2] = (138.f/5) + t*((-8617.f/60) + t*((12873.f/40) + t*((-791.f/2)
01057                 + t*((4557.f/16) + t*((-9583.f/80) + t*((2181.f/80)
01058                       + t*(-623.f/240)))))));
01059         dphi[2] = ((-8617.f/60) + t*((12873.f/20) + t*((-2373.f/2)
01060                 + t*((4557.f/4) + t*((-9583.f/16) + t*((6543.f/40)
01061                       + t*(-4361.f/240))))))) * h_1;
01062         t--;
01063         phi[3] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((2569.f/144)
01064                 + t*((-727.f/48) + t*(623.f/144)))));
01065         dphi[3] = (t*((-49.f/18) + t*t*((-959.f/36) + t*((12845.f/144)
01066                   + t*((-727.f/8) + t*(4361.f/144)))))) * h_1;
01067         t--;
01068         phi[4] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((-2569.f/144)
01069                 + t*((-727.f/48) + t*(-623.f/144)))));
01070         dphi[4] = (t*((-49.f/18) + t*t*((-959.f/36) + t*((-12845.f/144)
01071                   + t*((-727.f/8) + t*(-4361.f/144)))))) * h_1;
01072         t--;
01073         phi[5] = (138.f/5) + t*((8617.f/60) + t*((12873.f/40) + t*((791.f/2)
01074                 + t*((4557.f/16) + t*((9583.f/80) + t*((2181.f/80)
01075                       + t*(623.f/240)))))));
01076         dphi[5] = ((8617.f/60) + t*((12873.f/20) + t*((2373.f/2)
01077                 + t*((4557.f/4) + t*((9583.f/16) + t*((6543.f/40)
01078                       + t*(4361.f/240))))))) * h_1;
01079         t--;
01080         phi[6] = -440 + t*((-25949.f/20) + t*((-117131.f/72) + t*((-2247.f/2)
01081                 + t*((-66437.f/144) + t*((-81109.f/720) + t*((-727.f/48)
01082                       + t*(-623.f/720)))))));
01083         dphi[6] = ((-25949.f/20) + t*((-117131.f/36) + t*((-6741.f/2)
01084                 + t*((-66437.f/36) + t*((-81109.f/144) + t*((-727.f/8)
01085                       + t*(-4361.f/720))))))) * h_1;
01086         t--;
01087         phi[7] = (3632.f/5) + t*((7456.f/5) + t*((58786.f/45) + t*(633
01088                 + t*((26383.f/144) + t*((22807.f/720) + t*((727.f/240)
01089                       + t*(89.f/720)))))));
01090         dphi[7] = ((7456.f/5) + t*((117572.f/45) + t*(1899
01091                 + t*((26383.f/36) + t*((22807.f/144) + t*((727.f/40)
01092                       + t*(623.f/720))))))) * h_1;
01093         break;
01094       case NONIC:
01095         phi[0] = (-1.f/40320)*(t-8)*(t-7)*(t-6)*(t-5)*(t-5)*(t-4)*(t-3)*
01096           (t-2)*(t-1);
01097         dphi[0] = (-1.f/40320)*(t-5)*(-117648+t*(256552+t*(-221416
01098                 +t*(99340+t*(-25261+t*(3667+t*(-283+t*9)))))))*h_1;
01099         t--;
01100         phi[1] = (1.f/40320)*(t-7)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*
01101           (-8+t*(-35+9*t));
01102         dphi[1] = (1.f/40320)*(71856+t*(-795368+t*(1569240+t*(-1357692
01103                   +t*(634725+t*(-172116+t*(27090+t*(-2296+t*81))))))))*h_1;
01104         t--;
01105         phi[2] = (-1.f/10080)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*
01106           (-14+t*(-25+9*t));
01107         dphi[2] = (1.f/10080)*(3384+t*(-69080+t*(55026
01108                 +t*(62580+t*(-99225+t*(51660+t*(-13104+t*(1640
01109                           +t*(-81)))))))))*h_1;
01110         t--;
01111         phi[3] = (1.f/1440)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*
01112           (-6+t*(-5+3*t));
01113         dphi[3] = (1.f/1440)*(72+t*(-6344+t*(2070
01114                 +t*(7644+t*(-4725+t*(-828+t*(1260+t*(-328+t*27))))))))*h_1;
01115         t--;
01116         phi[4] = (-1.f/2880)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*
01117           (-20+t*(-5+9*t));
01118         dphi[4] = (-1.f/2880)*t*(10792+t*(-972+t*(-12516
01119                 +t*(2205+t*(3924+t*(-882+t*(-328+t*81)))))))*h_1;
01120         t--;
01121         phi[5] = (1.f/2880)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*
01122           (-20+t*(5+9*t));
01123         dphi[5] = (1.f/2880)*t*(-10792+t*(-972+t*(12516
01124                 +t*(2205+t*(-3924+t*(-882+t*(328+t*81)))))))*h_1;
01125         t--;
01126         phi[6] = (-1.f/1440)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*
01127           (-6+t*(5+3*t));
01128         dphi[6] = (1.f/1440)*(-72+t*(-6344+t*(-2070
01129                 +t*(7644+t*(4725+t*(-828+t*(-1260+t*(-328+t*(-27)))))))))*h_1;
01130         t--;
01131         phi[7] = (1.f/10080)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*
01132           (-14+t*(25+9*t));
01133         dphi[7] = (1.f/10080)*(-3384+t*(-69080+t*(-55026
01134                 +t*(62580+t*(99225+t*(51660+t*(13104+t*(1640+t*81))))))))*h_1;
01135         t--;
01136         phi[8] = (-1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*(t+7)*
01137           (-8+t*(35+9*t));
01138         dphi[8] = (-1.f/40320)*(71856+t*(795368+t*(1569240
01139                 +t*(1357692+t*(634725+t*(172116+t*(27090+t*(2296
01140                           +t*81))))))))*h_1;
01141         t--;
01142         phi[9] = (1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+5)*(t+6)*
01143           (t+7)*(t+8);
01144         dphi[9] = (1.f/40320)*(t+5)*(117648+t*(256552+t*(221416
01145                 +t*(99340+t*(25261+t*(3667+t*(283+t*9)))))))*h_1;
01146         break;
01147       case NONIC4:
01148       { // begin grouping to define local variables
01149         double Tphi[10], Tdphi[10];
01150         double T=t;
01151         Tphi[0] = 439375./7+T*(-64188125./504+T*(231125375./2016
01152               +T*(-17306975./288+T*(7761805./384+T*(-2895587./640
01153                     +T*(129391./192+T*(-259715./4032+T*(28909./8064
01154                           +T*(-3569./40320)))))))));
01155         Tdphi[0] = (-64188125./504+T*(231125375./1008
01156               +T*(-17306975./96+T*(7761805./96+T*(-2895587./128
01157                     +T*(129391./32+T*(-259715./576+T*(28909./1008
01158                           +T*(-3569./4480))))))))) * h_1;
01159         T--;
01160         Tphi[1] = -56375+T*(8314091./56+T*(-49901303./288+T*(3763529./32
01161                 +T*(-19648027./384+T*(9469163./640+T*(-545977./192
01162                       +T*(156927./448+T*(-28909./1152
01163                           +T*(3569./4480)))))))));
01164         Tdphi[1] = (8314091./56+T*(-49901303./144+T*(11290587./32
01165                 +T*(-19648027./96+T*(9469163./128+T*(-545977./32
01166                       +T*(156927./64+T*(-28909./144
01167                           +T*(32121./4480))))))))) * h_1;
01168         T--;
01169         Tphi[2] = 68776./7+T*(-1038011./28+T*(31157515./504+T*(-956669./16
01170                 +T*(3548009./96+T*(-2422263./160+T*(197255./48
01171                       +T*(-19959./28+T*(144545./2016
01172                           +T*(-3569./1120)))))))));
01173         Tdphi[2] = (-1038011./28+T*(31157515./252+T*(-2870007./16
01174                 +T*(3548009./24+T*(-2422263./32+T*(197255./8
01175                       +T*(-19959./4+T*(144545./252
01176                           +T*(-32121./1120))))))))) * h_1;
01177         T--;
01178         Tphi[3] = -154+T*(12757./12+T*(-230123./72+T*(264481./48
01179                 +T*(-576499./96+T*(686147./160+T*(-96277./48
01180                       +T*(14221./24+T*(-28909./288+T*(3569./480)))))))));
01181         Tdphi[3] = (12757./12+T*(-230123./36+T*(264481./16
01182                 +T*(-576499./24+T*(686147./32+T*(-96277./8
01183                       +T*(99547./24+T*(-28909./36
01184                           +T*(10707./160))))))))) * h_1;
01185         T--;
01186         Tphi[4] = 1+T*T*(-205./144+T*T*(91./192+T*(-6181./320
01187                 +T*(6337./96+T*(-2745./32+T*(28909./576
01188                       +T*(-3569./320)))))));
01189         Tdphi[4] = T*(-205./72+T*T*(91./48+T*(-6181./64
01190                 +T*(6337./16+T*(-19215./32+T*(28909./72
01191                       +T*(-32121./320))))))) * h_1;
01192         T--;
01193         Tphi[5] = 1+T*T*(-205./144+T*T*(91./192+T*(6181./320
01194                 +T*(6337./96+T*(2745./32+T*(28909./576
01195                       +T*(3569./320)))))));
01196         Tdphi[5] = T*(-205./72+T*T*(91./48+T*(6181./64
01197                 +T*(6337./16+T*(19215./32+T*(28909./72
01198                       +T*(32121./320))))))) * h_1;
01199         T--;
01200         Tphi[6] = -154+T*(-12757./12+T*(-230123./72+T*(-264481./48
01201                 +T*(-576499./96+T*(-686147./160+T*(-96277./48
01202                       +T*(-14221./24+T*(-28909./288+T*(-3569./480)))))))));
01203         Tdphi[6] = (-12757./12+T*(-230123./36+T*(-264481./16
01204                 +T*(-576499./24+T*(-686147./32+T*(-96277./8
01205                       +T*(-99547./24+T*(-28909./36
01206                           +T*(-10707./160))))))))) * h_1;
01207         T--;
01208         Tphi[7] = 68776./7+T*(1038011./28+T*(31157515./504+T*(956669./16
01209                 +T*(3548009./96+T*(2422263./160+T*(197255./48
01210                       +T*(19959./28+T*(144545./2016+T*(3569./1120)))))))));
01211         Tdphi[7] = (1038011./28+T*(31157515./252+T*(2870007./16
01212                 +T*(3548009./24+T*(2422263./32+T*(197255./8
01213                       +T*(19959./4+T*(144545./252
01214                           +T*(32121./1120))))))))) * h_1;
01215         T--;
01216         Tphi[8] = -56375+T*(-8314091./56+T*(-49901303./288+T*(-3763529./32
01217                 +T*(-19648027./384+T*(-9469163./640+T*(-545977./192
01218                       +T*(-156927./448+T*(-28909./1152
01219                           +T*(-3569./4480)))))))));
01220         Tdphi[8] = (-8314091./56+T*(-49901303./144+T*(-11290587./32
01221                 +T*(-19648027./96+T*(-9469163./128+T*(-545977./32
01222                       +T*(-156927./64+T*(-28909./144
01223                           +T*(-32121./4480))))))))) * h_1;
01224         T--;
01225         Tphi[9] = 439375./7+T*(64188125./504+T*(231125375./2016
01226               +T*(17306975./288+T*(7761805./384+T*(2895587./640
01227                     +T*(129391./192+T*(259715./4032+T*(28909./8064
01228                           +T*(3569./40320)))))))));
01229         Tdphi[9] = (64188125./504+T*(231125375./1008
01230               +T*(17306975./96+T*(7761805./96+T*(2895587./128
01231                     +T*(129391./32+T*(259715./576+T*(28909./1008
01232                           +T*(3569./4480))))))))) * h_1;
01233         for (int i=0;  i < 10;  i++) {
01234           phi[i] = Float(Tphi[i]);
01235           dphi[i] = Float(Tdphi[i]);
01236         }
01237       } // end grouping to define local variables
01238         break;
01239       default:
01240         NAMD_die("Unknown MSM approximation.");
01241     } // switch
01242   } // d_stencil_1d()

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

01254                                    {
01255     phi[0] = (1 - t) * (1 - t) * (1 + 2*t);
01256     dphi[0] = -6 * t * (1 - t) * h_1;
01257     psi[0] = h * t * (1 - t) * (1 - t);
01258     dpsi[0] = (1 - t) * (1 - 3*t);
01259     t--;
01260     phi[1] = (1 + t) * (1 + t) * (1 - 2*t);
01261     dphi[1] = -6 * t * (1 + t) * h_1;
01262     psi[1] = h * t * (1 + t) * (1 + t);
01263     dpsi[1] = (1 + t) * (1 + 3*t);
01264   }

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

06041 {
06042   msmCompute->saveResults();
06043 }

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, initVirialContrib(), j, numVirialContrib, virial, VMAX, VXX, VXY, VXZ, VYY, VYZ, and VZZ.

00539                            {
00540     if (++cntVirialContrib >= numVirialContrib) {
00541       // reduce all gvsum contributions into virial tensor
00542       for (int n = 0;  n < VMAX;  n++) { virial[n] = 0; }
00543       int ia = gvsum.ia();
00544       int ib = gvsum.ib();
00545       int ja = gvsum.ja();
00546       int jb = gvsum.jb();
00547       int ka = gvsum.ka();
00548       int kb = gvsum.kb();
00549       for (int k = ka;  k <= kb;  k++) {
00550         for (int j = ja;  j <= jb;  j++) {
00551           for (int i = ia;  i <= ib;  i++) {
00552             Float cu = Float(i);
00553             Float cv = Float(j);
00554             Float cw = Float(k);
00555             Float c = gvsum(i,j,k);
00556             Float vx = cu*hufx + cv*hvfx + cw*hwfx;
00557             Float vy = cu*hufy + cv*hvfy + cw*hwfy;
00558             Float vz = cu*hufz + cv*hvfz + cw*hwfz;
00559             virial[VXX] -= c * vx * vx;
00560             virial[VXY] -= c * vx * vy;
00561             virial[VXZ] -= c * vx * vz;
00562             virial[VYY] -= c * vy * vy;
00563             virial[VYZ] -= c * vy * vz;
00564             virial[VZZ] -= c * vz * vz;
00565           }
00566         }
00567       }
00568       initVirialContrib();
00569     }
00570   }

static void ComputeMsmMgr::gc_c1hermite_elem_accum ( C1Matrix matrix,
BigReal  _c,
Vector  rv,
BigReal  _a,
int  _split 
) [inline, static]

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

01411                                          {
01412     const BigReal a_1 = 1./_a;
01413     const BigReal a_2 = a_1 * a_1;
01414     const BigReal s = (rv * rv) * a_2;
01415     const BigReal dx = -2 * rv.x * a_2;  // ds/dx
01416     const BigReal dy = -2 * rv.y * a_2;  // ds/dy
01417     const BigReal dz = -2 * rv.z * a_2;  // ds/dz
01418     const BigReal dd = 2 * a_2;  // d^2s/dx^2 = d^2s/dy^2 = d^2s/dz^2
01419     BigReal tmp;
01420     enum { nderiv = C1_VECTOR_SIZE-1 };
01421     BigReal p[nderiv];
01422     Float *g = matrix.melem;
01423 
01424     // multiply entire matrix by this coefficient
01425     _c = _c * a_1;
01426 
01427     // compute derivatives (d/ds)^k of splitting g(s), s=r^2
01428     ndsplitting(p, s, nderiv, _split);
01429 
01430     // weight 0
01431     tmp = _c * p[0];
01432     g[C1INDEX(D000,D000)] += tmp;
01433 
01434     // weight 1
01435     tmp = _c * p[1] * dx;
01436     g[C1INDEX(D100,D000)] += tmp;
01437     g[C1INDEX(D000,D100)] -= tmp;
01438 
01439     tmp = _c * p[1] * dy;
01440     g[C1INDEX(D010,D000)] += tmp;
01441     g[C1INDEX(D000,D010)] -= tmp;
01442 
01443     tmp = _c * p[1] * dz;
01444     g[C1INDEX(D001,D000)] += tmp;
01445     g[C1INDEX(D000,D001)] -= tmp;
01446 
01447     // C1 splitting returns here
01448 
01449     // weight 2
01450     tmp = _c * p[2] * dx * dy;
01451     g[C1INDEX(D110,D000)] += tmp;
01452     g[C1INDEX(D000,D110)] += tmp;
01453     g[C1INDEX(D100,D010)] -= tmp;
01454     g[C1INDEX(D010,D100)] -= tmp;
01455 
01456     tmp = _c * p[2] * dx * dz;
01457     g[C1INDEX(D101,D000)] += tmp;
01458     g[C1INDEX(D000,D101)] += tmp;
01459     g[C1INDEX(D100,D001)] -= tmp;
01460     g[C1INDEX(D001,D100)] -= tmp;
01461 
01462     tmp = _c * p[2] * dy * dz;
01463     g[C1INDEX(D011,D000)] += tmp;
01464     g[C1INDEX(D000,D011)] += tmp;
01465     g[C1INDEX(D010,D001)] -= tmp;
01466     g[C1INDEX(D001,D010)] -= tmp;
01467 
01468     tmp = _c * (p[2] * dx*dx + p[1] * dd);
01469     g[C1INDEX(D100,D100)] -= tmp;
01470     tmp = _c * (p[2] * dy*dy + p[1] * dd);
01471     g[C1INDEX(D010,D010)] -= tmp;
01472     tmp = _c * (p[2] * dz*dz + p[1] * dd);
01473     g[C1INDEX(D001,D001)] -= tmp;
01474 
01475     // C2 splitting returns here
01476     if (_split == TAYLOR2) return;
01477 
01478     // weight 3
01479     tmp = _c * p[3] * dx * dy * dz;
01480     g[C1INDEX(D111,D000)] += tmp;
01481     g[C1INDEX(D110,D001)] -= tmp;
01482     g[C1INDEX(D101,D010)] -= tmp;
01483     g[C1INDEX(D011,D100)] -= tmp;
01484     g[C1INDEX(D100,D011)] += tmp;
01485     g[C1INDEX(D010,D101)] += tmp;
01486     g[C1INDEX(D001,D110)] += tmp;
01487     g[C1INDEX(D000,D111)] -= tmp;
01488 
01489     tmp = _c * (p[3] * dx*dx * dy + p[2] * dd * dy);
01490     g[C1INDEX(D110,D100)] -= tmp;
01491     g[C1INDEX(D100,D110)] += tmp;
01492 
01493     tmp = _c * (p[3] * dx*dx * dz + p[2] * dd * dz);
01494     g[C1INDEX(D101,D100)] -= tmp;
01495     g[C1INDEX(D100,D101)] += tmp;
01496 
01497     tmp = _c * (p[3] * dy*dy * dx + p[2] * dd * dx);
01498     g[C1INDEX(D110,D010)] -= tmp;
01499     g[C1INDEX(D010,D110)] += tmp;
01500 
01501     tmp = _c * (p[3] * dy*dy * dz + p[2] * dd * dz);
01502     g[C1INDEX(D011,D010)] -= tmp;
01503     g[C1INDEX(D010,D011)] += tmp;
01504 
01505     tmp = _c * (p[3] * dz*dz * dx + p[2] * dd * dx);
01506     g[C1INDEX(D101,D001)] -= tmp;
01507     g[C1INDEX(D001,D101)] += tmp;
01508 
01509     tmp = _c * (p[3] * dz*dz * dy + p[2] * dd * dy);
01510     g[C1INDEX(D011,D001)] -= tmp;
01511     g[C1INDEX(D001,D011)] += tmp;
01512 
01513     // C3 splitting returns here
01514     if (_split == TAYLOR3) return;
01515 
01516     // weight 4
01517     tmp = _c * (p[4] * dx*dx * dy * dz + p[3] * dd * dy * dz);
01518     g[C1INDEX(D111,D100)] -= tmp;
01519     g[C1INDEX(D100,D111)] -= tmp;
01520     g[C1INDEX(D110,D101)] += tmp;
01521     g[C1INDEX(D101,D110)] += tmp;
01522 
01523     tmp = _c * (p[4] * dy*dy * dx * dz + p[3] * dd * dx * dz);
01524     g[C1INDEX(D111,D010)] -= tmp;
01525     g[C1INDEX(D010,D111)] -= tmp;
01526     g[C1INDEX(D110,D011)] += tmp;
01527     g[C1INDEX(D011,D110)] += tmp;
01528 
01529     tmp = _c * (p[4] * dz*dz * dx * dy + p[3] * dd * dx * dy);
01530     g[C1INDEX(D111,D001)] -= tmp;
01531     g[C1INDEX(D001,D111)] -= tmp;
01532     g[C1INDEX(D101,D011)] += tmp;
01533     g[C1INDEX(D011,D101)] += tmp;
01534 
01535     tmp = _c * (p[4] * dx*dx * dy*dy + p[3] * dx*dx * dd
01536         + p[3] * dd * dy*dy + p[2] * dd * dd);
01537     g[C1INDEX(D110,D110)] += tmp;
01538     tmp = _c * (p[4] * dx*dx * dz*dz + p[3] * dx*dx * dd
01539         + p[3] * dd * dz*dz + p[2] * dd * dd);
01540     g[C1INDEX(D101,D101)] += tmp;
01541     tmp = _c * (p[4] * dy*dy * dz*dz + p[3] * dy*dy * dd
01542         + p[3] * dd * dz*dz + p[2] * dd * dd);
01543     g[C1INDEX(D011,D011)] += tmp;
01544 
01545     // C4 splitting returns here
01546     if (_split == TAYLOR4) return;
01547 
01548     // weight 5
01549     tmp = _c * (p[5] * dx*dx * dy*dy * dz + p[4] * dx*dx * dd * dz
01550         + p[4] * dd * dy*dy * dz + p[3] * dd * dd * dz);
01551     g[C1INDEX(D111,D110)] += tmp;
01552     g[C1INDEX(D110,D111)] -= tmp;
01553 
01554     tmp = _c * (p[5] * dx*dx * dz*dz * dy + p[4] * dx*dx * dd * dy
01555         + p[4] * dd * dz*dz * dy + p[3] * dd * dd * dy);
01556     g[C1INDEX(D111,D101)] += tmp;
01557     g[C1INDEX(D101,D111)] -= tmp;
01558 
01559     tmp = _c * (p[5] * dy*dy * dz*dz * dx + p[4] * dy*dy * dd * dx
01560         + p[4] * dd * dz*dz * dx + p[3] * dd * dd * dx);
01561     g[C1INDEX(D111,D011)] += tmp;
01562     g[C1INDEX(D011,D111)] -= tmp;
01563 
01564     // C5 splitting returns here
01565     if (_split == TAYLOR5) return;
01566 
01567     // weight 6
01568     tmp = _c * (p[6] * dx*dx * dy*dy * dz*dz + p[5] * dx*dx * dy*dy * dd
01569         + p[5] * dx*dx * dd * dz*dz + p[5] * dd * dy*dy * dz*dz
01570         + p[4] * dx*dx * dd * dd + p[4] * dd * dy*dy * dd
01571         + p[4] * dd * dd * dz*dz + p[3] * dd * dd * dd);
01572     g[C1INDEX(D111,D111)] -= tmp;
01573 
01574     // calculate full matrix for C6 or higher splitting
01575 
01576   } // gc_c1hermite_elem_accum()

void ComputeMsmMgr::initialize ( MsmInitMsg  ) 

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, c, Lattice::c(), C1HERMITE, C1INDEX, Lattice::c_p(), Lattice::c_r(), cross(), CUBIC, 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, iINFO(), iout, msm::Map::ispx, msm::Map::ispy, msm::Map::ispz, j, lattice, Vector::length(), map, MSM_MAX_BLOCK_VOLUME, NAMD_die(), nhx, nhy, nhz, nlevels, NONIC, NONIC4, Nstencil, NUM_APPROX, NUM_SPLIT, PatchMap::numPatches(), PatchMap::Object(), Node::Object(), omega, padding, msm::Map::patchList, PhiStencil, PolyDegree, QUINTIC, QUINTIC2, ru, rv, rw, s_edge, Lattice::scale(), SEPTIC, SEPTIC3, Vector::set(), 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.

03966 {
03967 #ifdef DEBUG_MSM_VERBOSE
03968   printf("ComputeMsmMgr:  initialize() PE %d\n", CkMyPe());
03969 #endif
03970 
03971   smin = msg->smin;
03972   smax = msg->smax;
03973   delete msg;
03974 
03975 #if 0
03976   printf("PE%d: initializing MSM\n", CkMyPe());
03977 #endif
03978 
03979   SimParameters *simParams = Node::Object()->simParameters;
03980 
03981   // get required sim params, check validity
03982   lattice = simParams->lattice;
03983 
03984   // set user-defined extent of system
03985   Vector rmin(simParams->MSMxmin, simParams->MSMymin, simParams->MSMzmin);
03986   Vector rmax(simParams->MSMxmax, simParams->MSMymax, simParams->MSMzmax);
03987   Vector sdmin = lattice.scale(rmin);
03988   Vector sdmax = lattice.scale(rmax);
03989   // swap coordinates between min and max to correct for possible rotation
03990   if (sdmin.x > sdmax.x) { double t=sdmin.x; sdmin.x=sdmax.x; sdmax.x=t; }
03991   if (sdmin.y > sdmax.y) { double t=sdmin.y; sdmin.y=sdmax.y; sdmax.y=t; }
03992   if (sdmin.z > sdmax.z) { double t=sdmin.z; sdmin.z=sdmax.z; sdmax.z=t; }
03993   // extend smin, smax by user-defined extent, where appropriate
03994   if ( ! lattice.a_p() && (sdmin.x != 0 || sdmax.x != 0)) {
03995     if (sdmin.x < smin.x) {
03996       smin.x = sdmin.x;
03997       if (CkMyPe() == 0) {
03998         iout << iINFO << "MSM extending minimum X to "
03999           << simParams->MSMxmin << " A\n" << endi;
04000       }
04001     }
04002     if (sdmax.x > smax.x) {
04003       smax.x = sdmax.x;
04004       if (CkMyPe() == 0) {
04005         iout << iINFO << "MSM extending maximum X to "
04006           << simParams->MSMxmax << " A\n" << endi;
04007       }
04008     }
04009   }
04010   if ( ! lattice.b_p() && (sdmin.y != 0 || sdmax.y != 0)) {
04011     if (sdmin.y < smin.y) {
04012       smin.y = sdmin.y;
04013       if (CkMyPe() == 0) {
04014         iout << iINFO << "MSM extending minimum Y to "
04015           << simParams->MSMymin << " A\n" << endi;
04016       }
04017     }
04018     if (sdmax.y > smax.y) {
04019       smax.y = sdmax.y;
04020       if (CkMyPe() == 0) {
04021         iout << iINFO << "MSM extending maximum Y to "
04022           << simParams->MSMymax << " A\n" << endi;
04023       }
04024     }
04025   }
04026   if ( ! lattice.c_p() && (sdmin.z != 0 || sdmax.z != 0)) {
04027     if (sdmin.z < smin.z) {
04028       smin.z = sdmin.z;
04029       if (CkMyPe() == 0) {
04030         iout << iINFO << "MSM extending minimum Z to "
04031           << simParams->MSMzmin << " A\n" << endi;
04032       }
04033     }
04034     if (sdmax.z > smax.z) {
04035       smax.z = sdmax.z;
04036       if (CkMyPe() == 0) {
04037         iout << iINFO << "MSM extending maximum Z to "
04038           << simParams->MSMzmax << " A\n" << endi;
04039       }
04040     }
04041   }
04042 
04043 #ifdef DEBUG_MSM_VERBOSE
04044   printf("smin = %g %g %g  smax = %g %g %g\n",
04045       smin.x, smin.y, smin.z, smax.x, smax.y, smax.z);
04046 #endif
04047 
04048   approx = simParams->MSMApprox;
04049   if (approx < 0 || approx >= NUM_APPROX) {
04050     NAMD_die("MSM: unknown approximation requested (MSMApprox)");
04051   }
04052 
04053   split = simParams->MSMSplit;
04054   if (split < 0 || split >= NUM_SPLIT) {
04055     NAMD_die("MSM: unknown splitting requested (MSMSplit)");
04056   }
04057 
04058   if (CkMyPe() == 0) {
04059     const char *approx_str, *split_str;
04060     switch (approx) {
04061       case CUBIC:      approx_str = "C1 cubic";    break;
04062       case QUINTIC:    approx_str = "C1 quintic";  break;
04063       case QUINTIC2:   approx_str = "C2 quintic";  break;
04064       case SEPTIC:     approx_str = "C1 septic";   break;
04065       case SEPTIC3:    approx_str = "C3 septic";   break;
04066       case NONIC:      approx_str = "C1 nonic";    break;
04067       case NONIC4:     approx_str = "C4 nonic";    break;
04068       case C1HERMITE:  approx_str = "C1 Hermite";  break;
04069       default:         approx_str = "unknown";     break;
04070     }
04071     switch (split) {
04072       case TAYLOR2:  split_str = "C2 Taylor";   break;
04073       case TAYLOR3:  split_str = "C3 Taylor";   break;
04074       case TAYLOR4:  split_str = "C4 Taylor";   break;
04075       case TAYLOR5:  split_str = "C5 Taylor";   break;
04076       case TAYLOR6:  split_str = "C6 Taylor";   break;
04077       case TAYLOR7:  split_str = "C7 Taylor";   break;
04078       case TAYLOR8:  split_str = "C8 Taylor";   break;
04079       default:       split_str = "unknown";     break;
04080     }
04081     iout << iINFO << "MSM using "
04082                   << approx_str << " interpolation\n";
04083     iout << iINFO << "MSM using "
04084                   << split_str << " splitting function\n";
04085   }
04086 
04087   a = simParams->cutoff;
04088 
04089   if (approx == C1HERMITE) {
04090     gridScalingFactor = 2;
04091   }
04092   else {
04093     gridScalingFactor = 1;
04094   }
04095 
04096   gridspacing = gridScalingFactor * simParams->MSMGridSpacing;
04097   if (gridspacing <= 0) {
04098     NAMD_die("MSM: grid spacing must be greater than 0 (MSMGridSpacing)");
04099   }
04100   else if (gridspacing >= a) {
04101     NAMD_die("MSM: grid spacing must be less than cutoff (MSMGridSpacing)");
04102   }
04103 
04104   padding = gridScalingFactor * simParams->MSMPadding;
04105   if (padding < 0) {
04106     NAMD_die("MSM: padding must be non-negative (MSMPadding)");
04107   }
04108 
04109   // set maximum number of levels (default 0 adapts levels to system)
04110   nlevels = simParams->MSMLevels;
04111 
04112   // XXX dispersion unused for now
04113   dispersion = 0;
04114   if ( ! dispersion && split >= TAYLOR2_DISP) {
04115     NAMD_die("MSM: requested splitting for long-range dispersion "
04116         "(not implemented)");
04117   }
04118 
04119   // set block sizes for grid decomposition
04120   int bsx = simParams->MSMBlockSizeX / int(gridScalingFactor);
04121   int bsy = simParams->MSMBlockSizeY / int(gridScalingFactor);
04122   int bsz = simParams->MSMBlockSizeZ / int(gridScalingFactor);
04123   if (bsx <= 0 || bsy <= 0 || bsz <= 0) {
04124     NAMD_die("MSM: invalid block size requested (MSMBlockSize[XYZ])");
04125   }
04126 #ifdef MSM_FIXED_SIZE_GRID_MSG
04127   else if (bsx * bsy * bsz > MSM_MAX_BLOCK_VOLUME) {
04128     NAMD_die("MSM: requested block size (MSMBlockSize[XYZ]) too big");
04129   }
04130 #endif
04131   if (CkMyPe() == 0) {
04132     iout << iINFO << "MSM block size decomposition along X is "
04133                   << bsx << " grid points\n";
04134     iout << iINFO << "MSM block size decomposition along Y is "
04135                   << bsy << " grid points\n";
04136     iout << iINFO << "MSM block size decomposition along Z is "
04137                   << bsz << " grid points\n";
04138   }
04139 
04140   s_edge = (PolyDegree[approx] - 1) / 2;  // stencil edge size
04141   omega = 2 * PolyDegree[approx];         // smallest non-periodic grid length
04142 
04143   BigReal xlen, ylen, zlen;
04144   Vector sgmin, sgmax;  // grid min and max, in scaled coordinates
04145   int ispx = lattice.a_p();
04146   int ispy = lattice.b_p();
04147   int ispz = lattice.c_p();
04148   int ispany = (ispx || ispy || ispz);  // is there any periodicity?
04149 
04150   if (ispx) {  // periodic along basis vector
04151     xlen = lattice.a().length();
04152     sgmax.x = 0.5;
04153     sgmin.x = -0.5;
04154   }
04155   else {  // non-periodic
04156     sgmax.x = smax.x + padding;  // pad the edges
04157     sgmin.x = smin.x - padding;
04158     ASSERT(gridspacing > 0);
04159     // restrict center to be on a grid point
04160     BigReal mupper = ceil(sgmax.x / (2*gridspacing));
04161     BigReal mlower = floor(sgmin.x / (2*gridspacing));
04162     sgmax.x = 2*gridspacing*mupper;
04163     sgmin.x = 2*gridspacing*mlower;
04164     xlen = sgmax.x - sgmin.x;
04165   }
04166 #ifdef DEBUG_MSM_VERBOSE
04167   printf("xlen = %g   sgmin.x = %g   sgmax.x = %g\n", xlen, sgmin.x, sgmax.x);
04168 #endif
04169 
04170   if (ispy) {  // periodic along basis vector
04171     ylen = lattice.b().length();
04172     sgmax.y = 0.5;
04173     sgmin.y = -0.5;
04174   }
04175   else {  // non-periodic
04176     sgmax.y = smax.y + padding;  // pad the edges
04177     sgmin.y = smin.y - padding;
04178     ASSERT(gridspacing > 0);
04179     // restrict center to be on a grid point
04180     BigReal mupper = ceil(sgmax.y / (2*gridspacing));
04181     BigReal mlower = floor(sgmin.y / (2*gridspacing));
04182     sgmax.y = 2*gridspacing*mupper;
04183     sgmin.y = 2*gridspacing*mlower;
04184     ylen = sgmax.y - sgmin.y;
04185   }
04186 #ifdef DEBUG_MSM_VERBOSE
04187   printf("ylen = %g   sgmin.y = %g   sgmax.y = %g\n", ylen, sgmin.y, sgmax.y);
04188 #endif
04189 
04190   if (ispz) {  // periodic along basis vector
04191     zlen = lattice.c().length();
04192     sgmax.z = 0.5;
04193     sgmin.z = -0.5;
04194   }
04195   else {  // non-periodic
04196     sgmax.z = smax.z + padding;  // pad the edges
04197     sgmin.z = smin.z - padding;
04198     ASSERT(gridspacing > 0);
04199     // restrict center to be on a grid point
04200     BigReal mupper = ceil(sgmax.z / (2*gridspacing));
04201     BigReal mlower = floor(sgmin.z / (2*gridspacing));
04202     sgmax.z = 2*gridspacing*mupper;
04203     sgmin.z = 2*gridspacing*mlower;
04204     zlen = sgmax.z - sgmin.z;
04205   }
04206 #ifdef DEBUG_MSM_VERBOSE
04207   printf("zlen = %g   sgmin.z = %g   sgmax.z = %g\n", zlen, sgmin.z, sgmax.z);
04208 #endif
04209   sglower = sgmin;
04210 
04211   int ia, ib, ja, jb, ka, kb;
04212   setup_hgrid_1d(xlen, hxlen, nhx, ia, ib, ispx);
04213   setup_hgrid_1d(ylen, hylen, nhy, ja, jb, ispy);
04214   setup_hgrid_1d(zlen, hzlen, nhz, ka, kb, ispz);
04215   hxlen_1 = 1 / hxlen;
04216   hylen_1 = 1 / hylen;
04217   hzlen_1 = 1 / hzlen;
04218   if (CkMyPe() == 0) {
04219     if (ispx || ispy || ispz) {
04220       iout << iINFO << "MSM grid spacing along X is "<< hxlen << " A\n";
04221       iout << iINFO << "MSM grid spacing along Y is "<< hylen << " A\n";
04222       iout << iINFO << "MSM grid spacing along Z is "<< hzlen << " A\n";
04223     }
04224     else {
04225       iout << iINFO << "MSM grid spacing is " << gridspacing << " A\n";
04226     }
04227     if ( ! ispx || ! ispy || ! ispz ) {
04228       iout << iINFO<<"MSM non-periodic padding is "<< padding << " A\n";
04229     }
04230   }
04231 
04232   int ni = ib - ia + 1;
04233   int nj = jb - ja + 1;
04234   int nk = kb - ka + 1;
04235   int n;
04236 
04237 #if 0
04238   // reserve temp space for factored grid transfer operation
04239   n = (nk > omega ? nk : omega);  // row along z-dimension
04240   lzd.resize(n);
04241   n *= (nj > omega ? nj : omega);  // plane along yz-dimensions
04242   lyzd.resize(n);
04243 #endif
04244 
04245   int lastnelems = 1;
04246   int smallestnbox = 1;
04247   int isclamped = 0;
04248   int maxlevels = nlevels;  // user-defined number of levels
04249 
04250 #ifdef DEBUG_MSM_VERBOSE
04251   printf("maxlevels = %d\n", maxlevels);
04252 #endif
04253   if (nlevels <= 0) {  // instead we set number of levels
04254     n = ni;
04255     if (n < nj) n = nj;
04256     if (n < nk) n = nk;
04257     for (maxlevels = 1;  n > 0;  n >>= 1)  maxlevels++;
04258     if (ispany == 0) {  // no periodicity
04259       // use rule of thumb 3/4 diameter of grid cutoff sphere
04260       int ngci = (int) ceil(3*a / hxlen) - 1;
04261       int ngcj = (int) ceil(3*a / hylen) - 1;
04262       int ngck = (int) ceil(3*a / hzlen) - 1;
04263       int omega3 = omega * omega * omega;
04264       int nhalf = (int) sqrt((double)ni * nj * nk);
04265       lastnelems = (nhalf > omega3 ? nhalf : omega3);
04266       smallestnbox = ngci * ngcj * ngck;  // smaller grids don't reduce work
04267       isclamped = 1;
04268     }
04269   }
04270 #ifdef DEBUG_MSM_VERBOSE
04271   printf("maxlevels = %d\n", maxlevels);
04272 #endif
04273 
04274   // allocate space for storing grid dimensions for each level
04275   map.gridrange.resize(maxlevels);
04276 
04277   // set periodicity flags
04278   map.ispx = ispx;
04279   map.ispy = ispy;
04280   map.ispz = ispz;
04281 
04282   int level = 0;
04283   int done = 0;
04284   int alldone = 0;
04285   do {
04286     map.gridrange[level].setbounds(ia, ib, ja, jb, ka, kb);
04287 
04288     // Msm index?
04289 
04290     if (++level == nlevels)  done |= 0x07;  // user limit on levels
04291 
04292     if (isclamped) {
04293       int nelems = ni * nj * nk;
04294       if (nelems <= lastnelems)    done |= 0x07;
04295       if (nelems <= smallestnbox)  done |= 0x07;
04296     }
04297 
04298     alldone = (done == 0x07);  // make sure all dimensions are done
04299 
04300     if (ispx) {
04301       ni >>= 1;
04302       ib = ni-1;
04303       if (ni & 1)        done |= 0x07;  // == 3 or 1
04304       else if (ni == 2)  done |= 0x01;  // can do one more
04305     }
04306     else {
04307       ia = -((-ia+1)/2) - s_edge;
04308       ib = (ib+1)/2 + s_edge;
04309       ni = ib - ia + 1;
04310       if (ni <= omega)   done |= 0x01;  // can do more restrictions
04311     }
04312 
04313     if (ispy) {
04314       nj >>= 1;
04315       jb = nj-1;
04316       if (nj & 1)        done |= 0x07;  // == 3 or 1
04317       else if (nj == 2)  done |= 0x02;  // can do one more
04318     }
04319     else {
04320       ja = -((-ja+1)/2) - s_edge;
04321       jb = (jb+1)/2 + s_edge;
04322       nj = jb - ja + 1;
04323       if (nj <= omega)   done |= 0x02;  // can do more restrictions
04324     }
04325 
04326     if (ispz) {
04327       nk >>= 1;
04328       kb = nk-1;
04329       if (nk & 1)        done |= 0x07;  // == 3 or 1
04330       else if (nk == 2)  done |= 0x04;  // can do one more
04331     }
04332     else {
04333       ka = -((-ka+1)/2) - s_edge;
04334       kb = (kb+1)/2 + s_edge;
04335       nk = kb - ka + 1;
04336       if (nk <= omega)   done |= 0x04;  // can do more restrictions
04337     }
04338   } while ( ! alldone );
04339   nlevels = level;
04340 
04341   // for periodic boundaries, don't visit top level (all 0)
04342   // toplevel visited only for all nonperiodic boundaries
04343   int toplevel = (ispany ? nlevels : nlevels - 1);
04344 
04345   // resize down to the actual number of levels (does not change alloc)
04346   map.gridrange.resize(nlevels);
04347 
04348   // print out some information about MSM
04349   if (CkMyPe() == 0) {
04350     iout << iINFO << "MSM using " << nlevels << " levels\n";
04351     for (n = 0;  n < nlevels;  n++) {
04352       char s[100];
04353       snprintf(s, sizeof(s), "    level %d:  "
04354           "[%d..%d] x [%d..%d] x [%d..%d]\n", n,
04355           map.gridrange[n].ia(), map.gridrange[n].ib(),
04356           map.gridrange[n].ja(), map.gridrange[n].jb(),
04357           map.gridrange[n].ka(), map.gridrange[n].kb());
04358       iout << iINFO << s;
04359     }
04360     iout << endi;
04361   }
04362 
04363   // find grid spacing basis vectors
04364   hu = hxlen * lattice.a().unit();
04365   hv = hylen * lattice.b().unit();
04366   hw = hzlen * lattice.c().unit();
04367   hufx = Float(hu.x);
04368   hufy = Float(hu.y);
04369   hufz = Float(hu.z);
04370   hvfx = Float(hv.x);
04371   hvfy = Float(hv.y);
04372   hvfz = Float(hv.z);
04373   hwfx = Float(hw.x);
04374   hwfy = Float(hw.y);
04375   hwfz = Float(hw.z);
04376 
04377   ru = lattice.a_r();
04378   rv = lattice.b_r();
04379   rw = lattice.c_r();
04380 
04381   // determine grid spacings in scaled space
04382   shx = ru * hu;
04383   shy = rv * hv;
04384   shz = rw * hw;
04385   shx_1 = 1 / shx;
04386   shy_1 = 1 / shy;
04387   shz_1 = 1 / shz;
04388 
04389   // row vectors to transform interpolated force back to real space
04390   // XXX Is not needed.
04391   sx_shx = shx_1 * Vector(ru.x, rv.x, rw.x);
04392   sy_shy = shy_1 * Vector(ru.y, rv.y, rw.y);
04393   sz_shz = shz_1 * Vector(ru.z, rv.z, rw.z);
04394   srx_x = Float(sx_shx.x);
04395   srx_y = Float(sx_shx.y);
04396   srx_z = Float(sx_shx.z);
04397   sry_x = Float(sy_shy.x);
04398   sry_y = Float(sy_shy.y);
04399   sry_z = Float(sy_shy.z);
04400   srz_x = Float(sz_shz.x);
04401   srz_y = Float(sz_shz.y);
04402   srz_z = Float(sz_shz.z);
04403 
04404   Vector pu = cross(hv, hw);
04405   BigReal s = (hu * pu) / (pu * pu);
04406   pu *= s;  // pu is orthogonal projection of hu onto hv CROSS hw
04407 
04408   Vector pv = cross(hw, hu);
04409   s = (hv * pv) / (pv * pv);
04410   pv *= s;  // pv is orthogonal projection of hv onto hw CROSS hu
04411 
04412   Vector pw = cross(hu, hv);
04413   s = (hw * pw) / (pw * pw);
04414   pw *= s;  // pw is orthogonal projection of hw onto hu CROSS hv
04415 
04416   // radii for parallelepiped of weights enclosing grid cutoff sphere
04417   ni = (int) ceil(2*a / pu.length() ) - 1;
04418   nj = (int) ceil(2*a / pv.length() ) - 1;
04419   nk = (int) ceil(2*a / pw.length() ) - 1;
04420 
04421   Float scaling = 1;
04422   Float scaling_factor = 0.5f;
04423   BigReal a_1 = 1/a;
04424   BigReal a_p = a_1;
04425   if (dispersion) {
04426     a_p = a_p * a_p * a_p;   // = 1/a^3
04427     a_p = a_p * a_p;         // = 1/a^6
04428     scaling_factor = 1.f/64;  // = 1/2^6
04429   }
04430   int i, j, k;
04431   if (approx < C1HERMITE) {
04432     // resize gc and gvc constants for number of levels
04433     map.gc.resize(nlevels);
04434     map.gvc.resize(nlevels);
04435 
04436     for (level = 0;  level < toplevel;  level++) {
04437       map.gc[level].setbounds(-ni, ni, -nj, nj, -nk, nk);
04438       map.gvc[level].setbounds(-ni, ni, -nj, nj, -nk, nk);
04439 
04440       for (k = -nk;  k <= nk;  k++) {
04441         for (j = -nj;  j <= nj;  j++) {
04442           for (i = -ni;  i <= ni;  i++) {
04443             if (level == 0) {
04444               BigReal s, t, gs=0, gt=0, g=0, dgs=0, dgt=0, dg=0;
04445               BigReal vlen = (i*hu + j*hv + k*hw).length();
04446               s = vlen * a_1;
04447               t = 0.5 * s;
04448               if (t >= 1) {
04449                 g = 0;
04450                 dg = 0;
04451               }
04452               else {
04453                 splitting(gt, dgt, t, split);
04454                 if (s >= 1) {
04455                   BigReal s_1 = 1/s;
04456                   if (dispersion) {
04457                     gs = s_1 * s_1 * s_1;  // = 1/s^3
04458                     gs = gs * gs;  // = 1/s^6
04459                     dgs = -6 * gs * s_1;
04460                   }
04461                   else {
04462                     gs = s_1;
04463                     dgs = -gs * s_1;
04464                   }
04465                 }
04466                 else {
04467                   splitting(gs, dgs, s, split);
04468                 }
04469                 g = (gs - scaling_factor * gt) * a_p;
04470                 BigReal c=0;
04471                 if (i || j || k) {
04472                   c = a_p * a_1 / vlen;
04473                 }
04474                 dg = 0.5 * (dgs - 0.5*scaling_factor * dgt) * c;
04475 
04476                 // Msm index?
04477 
04478               }
04479               map.gc[0](i,j,k) = Float(g);
04480               map.gvc[0](i,j,k) = Float(dg);
04481             } // if level 0
04482             else {
04483               map.gc[level](i,j,k) = scaling * map.gc[0](i,j,k);
04484               map.gvc[level](i,j,k) = scaling * map.gvc[0](i,j,k);
04485             }
04486 
04487           } // for i
04488         } // for j
04489       } // for k
04490       scaling *= scaling_factor;
04491 
04492     } // for level
04493 
04494     // for summed virial factors
04495     gvsum.setbounds(-ni, ni, -nj, nj, -nk, nk);
04496     // make sure final virial sum is initialized to 0
04497     for (i = 0;  i < VMAX;  i++) { virial[i] = 0; }
04498 
04499     if (toplevel < nlevels) {
04500       // nonperiodic along all basis vector directions
04501       // calculate top level weights where all grid points
04502       // interact with each other
04503       ni = map.gridrange[toplevel].ni();
04504       nj = map.gridrange[toplevel].nj();
04505       nk = map.gridrange[toplevel].nk();
04506       map.gc[toplevel].setbounds(-ni, ni, -nj, nj, -nk, nk);
04507 
04508       // Msm index?
04509 
04510       for (k = -nk;  k <= nk;  k++) {
04511         for (j = -nj;  j <= nj;  j++) {
04512           for (i = -ni;  i <= ni;  i++) {
04513             BigReal s, gs, d;
04514             BigReal vlen = (i*hu + j*hv + k*hw).length();
04515             s = vlen * a_1;
04516             if (s >= 1) {
04517               BigReal s_1 = 1/s;
04518               if (dispersion) {
04519                 gs = s_1 * s_1 * s_1;  // = 1/s^3
04520                 gs = gs * gs;  // = 1/s^6
04521               }
04522               else {
04523                 gs = s_1;
04524               }
04525             }
04526             else {
04527               splitting(gs, d, s, split);
04528             }
04529             map.gc[toplevel](i,j,k) = scaling * Float(gs * a_p);
04530           } // for i
04531         } // for j
04532       } // for k
04533     } // if toplevel
04534 
04535     // generate grespro stencil
04536     const int nstencil = Nstencil[approx];
04537     const Float *phi = PhiStencil[approx];
04538     map.grespro.set(0, nstencil, 0, nstencil, 0, nstencil);
04539     for (k = 0;  k < nstencil;  k++) {
04540       for (j = 0;  j < nstencil;  j++) {
04541         for (i = 0;  i < nstencil;  i++) {
04542           map.grespro(i,j,k) = phi[i] * phi[j] * phi[k];
04543         }
04544       }
04545     }
04546 
04547   } // end if approx < C1HERMITE
04548   else {
04549     // C1HERMITE
04550     // resize gc_c1hermite constants for number of levels
04551     map.gc_c1hermite.resize(nlevels);
04552     scaling = 1;
04553 
04554     for (level = 0;  level < toplevel;  level++) {
04555 
04556       Vector hmu = scaling * hu;
04557       Vector hmv = scaling * hv;
04558       Vector hmw = scaling * hw;
04559       BigReal am = scaling * a;
04560 
04561       map.gc_c1hermite[level].setbounds(-ni, ni, -nj, nj, -nk, nk);
04562 
04563       for (k = -nk;  k <= nk;  k++) {
04564         for (j = -nj;  j <= nj;  j++) {
04565           for (i = -ni;  i <= ni;  i++) {
04566             C1Matrix& m = map.gc_c1hermite[level](i,j,k);
04567             Vector rv = i*hmu + j*hmv + k*hmw;
04568             BigReal r2 = rv * rv;
04569             m.set(0);
04570             if (r2 < 4*am*am) {
04571               // accumulate D( g_{a}(0,r) ) term for this level
04572               gc_c1hermite_elem_accum(m, 1, rv, am, split);
04573               // accumulate D( -g_{2a}(0,r) ) term for this level
04574               gc_c1hermite_elem_accum(m, -1, rv, 2*am, split);
04575             } // if within cutoff
04576           }
04577         }
04578       } // end loop over gc_c1hermite elements for this level
04579       scaling *= 2;  // double grid spacing and cutoff at each iteration
04580 
04581     } // end loop over levels
04582 
04583     if (toplevel < nlevels) {
04584       Vector hmu = scaling * hu;
04585       Vector hmv = scaling * hv;
04586       Vector hmw = scaling * hw;
04587       BigReal am = scaling * a;
04588 
04589       // nonperiodic along all basis vector directions
04590       // calculate top level weights where all grid points
04591       // interact with each other
04592       ni = map.gridrange[toplevel].ni();
04593       nj = map.gridrange[toplevel].nj();
04594       nk = map.gridrange[toplevel].nk();
04595       map.gc_c1hermite[toplevel].setbounds(-ni, ni, -nj, nj, -nk, nk);
04596 
04597       for (k = -nk;  k <= nk;  k++) {
04598         for (j = -nj;  j <= nj;  j++) {
04599           for (i = -ni;  i <= ni;  i++) {
04600             C1Matrix& m = map.gc_c1hermite[level](i,j,k);
04601             Vector rv = i*hmu + j*hmv + k*hmw;
04602             m.set(0);
04603             // accumulate D( g_{a}(0,r) ) term for this level
04604             gc_c1hermite_elem_accum(m, 1, rv, am, split);
04605           }
04606         }
04607       } // end loop over gc_c1hermite elements for top level
04608 
04609     } // end if top level
04610 
04611     // C1 Hermite restriction and prolongation stencils
04612     map.gres_c1hermite.resize(nlevels-1);
04613     map.gpro_c1hermite.resize(nlevels-1);
04614 
04615     enum {
04616       ND = 3,    // stencil diameter
04617       NR = ND/2  // stencil radius
04618     };
04619 
04620     // the master basis functions PHI0 and PHI1 for the 3-point stencil
04621     // and their derivatives DPHI0 and DPHI1
04622     const double  PHI0[ND] = { 0.5, 1, 0.5 };
04623     const double DPHI0[ND] = { 1.5, 0, -1.5 };
04624     const double  PHI1[ND] = { -0.125, 0, 0.125 };
04625     const double DPHI1[ND] = { -0.25, 1, -0.25 };
04626 
04627     // for intermediate calculations
04628     double  xphi0_base_array[ND];
04629     double dxphi0_base_array[ND];
04630     double  yphi0_base_array[ND];
04631     double dyphi0_base_array[ND];
04632     double  zphi0_base_array[ND];
04633     double dzphi0_base_array[ND];
04634     double  xphi1_base_array[ND];
04635     double dxphi1_base_array[ND];
04636     double  yphi1_base_array[ND];
04637     double dyphi1_base_array[ND];
04638     double  zphi1_base_array[ND];
04639     double dzphi1_base_array[ND];
04640     // will point to center of stencil arrays
04641     double *xphi0, *dxphi0, *xphi1, *dxphi1;
04642     double *yphi0, *dyphi0, *yphi1, *dyphi1;
04643     double *zphi0, *dzphi0, *zphi1, *dzphi1;
04644 
04645     for (n = 0;  n < ND;  n++) {
04646       xphi0_base_array[n]  = PHI0[n];
04647       dxphi0_base_array[n] = hxlen_1 * DPHI0[n];  // scale by grid spacing
04648       xphi1_base_array[n]  = hxlen * PHI1[n];     // scale by grid spacing
04649       dxphi1_base_array[n] = DPHI1[n];
04650       yphi0_base_array[n]  = PHI0[n];
04651       dyphi0_base_array[n] = hylen_1 * DPHI0[n];  // scale by grid spacing
04652       yphi1_base_array[n]  = hylen * PHI1[n];     // scale by grid spacing
04653       dyphi1_base_array[n] = DPHI1[n];
04654       zphi0_base_array[n]  = PHI0[n];
04655       dzphi0_base_array[n] = hzlen_1 * DPHI0[n];  // scale by grid spacing
04656       zphi1_base_array[n]  = hzlen * PHI1[n];     // scale by grid spacing
04657       dzphi1_base_array[n] = DPHI1[n];
04658     }
04659     xphi0  =  xphi0_base_array + NR;  // point into center of arrays
04660     dxphi0 = dxphi0_base_array + NR;
04661     xphi1  =  xphi1_base_array + NR;
04662     dxphi1 = dxphi1_base_array + NR;
04663     yphi0  =  yphi0_base_array + NR;
04664     dyphi0 = dyphi0_base_array + NR;
04665     yphi1  =  yphi1_base_array + NR;
04666     dyphi1 = dyphi1_base_array + NR;
04667     zphi0  =  zphi0_base_array + NR;
04668     dzphi0 = dzphi0_base_array + NR;
04669     zphi1  =  zphi1_base_array + NR;
04670     dzphi1 = dzphi1_base_array + NR;
04671 
04672     for (level = 0;  level < nlevels-1;  level++) {
04673       // allocate space for restriction and prolongation stencils
04674       map.gres_c1hermite[level].set(0, ND, 0, ND, 0, ND);
04675       map.gpro_c1hermite[level].set(0, ND, 0, ND, 0, ND);
04676 
04677       // scale up to next level grid spacing
04678       //
04679       // have to determine for each dimension whether or not 
04680       // a periodic grid spacing has increased 
04681       // (equivalent to if there are fewer grid points)
04682       for (n = -NR;  n <= NR;  n++) {
04683         if ( ! ispx ||
04684               map.gridrange[level+1].ni() < map.gridrange[level].ni() ) {
04685           dxphi0[n] *= 0.5;
04686           xphi1[n] *= 2;
04687         }
04688         if ( ! ispy ||
04689               map.gridrange[level+1].nj() < map.gridrange[level].nj() ) {
04690           dyphi0[n] *= 0.5;
04691           yphi1[n] *= 2;
04692         }
04693         if ( ! ispz ||
04694               map.gridrange[level+1].nk() < map.gridrange[level].nk() ) {
04695           dzphi0[n] *= 0.5;
04696           zphi1[n] *= 2;
04697         }
04698       }
04699 
04700       // loop over restriction stencil matrices
04701       // calculate from partial derivatives
04702       for (k = -NR;  k <= NR;  k++) {
04703         for (j = -NR;  j <= NR;  j++) {
04704           for (i = -NR;  i <= NR;  i++) {
04705             Float *t = map.gres_c1hermite[level](i+NR,j+NR,k+NR).melem;
04706 
04707             t[C1INDEX(D000,D000)] =  xphi0[i] *  yphi0[j]  *  zphi0[k];
04708             t[C1INDEX(D000,D100)] = dxphi0[i] *  yphi0[j]  *  zphi0[k];
04709             t[C1INDEX(D000,D010)] =  xphi0[i] * dyphi0[j]  *  zphi0[k];
04710             t[C1INDEX(D000,D001)] =  xphi0[i] *  yphi0[j]  * dzphi0[k];
04711             t[C1INDEX(D000,D110)] = dxphi0[i] * dyphi0[j]  *  zphi0[k];
04712             t[C1INDEX(D000,D101)] = dxphi0[i] *  yphi0[j]  * dzphi0[k];
04713             t[C1INDEX(D000,D011)] =  xphi0[i] * dyphi0[j]  * dzphi0[k];
04714             t[C1INDEX(D000,D111)] = dxphi0[i] * dyphi0[j]  * dzphi0[k];
04715 
04716             t[C1INDEX(D100,D000)] =  xphi1[i] *  yphi0[j]  *  zphi0[k];
04717             t[C1INDEX(D100,D100)] = dxphi1[i] *  yphi0[j]  *  zphi0[k];
04718             t[C1INDEX(D100,D010)] =  xphi1[i] * dyphi0[j]  *  zphi0[k];
04719             t[C1INDEX(D100,D001)] =  xphi1[i] *  yphi0[j]  * dzphi0[k];
04720             t[C1INDEX(D100,D110)] = dxphi1[i] * dyphi0[j]  *  zphi0[k];
04721             t[C1INDEX(D100,D101)] = dxphi1[i] *  yphi0[j]  * dzphi0[k];
04722             t[C1INDEX(D100,D011)] =  xphi1[i] * dyphi0[j]  * dzphi0[k];
04723             t[C1INDEX(D100,D111)] = dxphi1[i] * dyphi0[j]  * dzphi0[k];
04724 
04725             t[C1INDEX(D010,D000)] =  xphi0[i] *  yphi1[j]  *  zphi0[k];
04726             t[C1INDEX(D010,D100)] = dxphi0[i] *  yphi1[j]  *  zphi0[k];
04727             t[C1INDEX(D010,D010)] =  xphi0[i] * dyphi1[j]  *  zphi0[k];
04728             t[C1INDEX(D010,D001)] =  xphi0[i] *  yphi1[j]  * dzphi0[k];
04729             t[C1INDEX(D010,D110)] = dxphi0[i] * dyphi1[j]  *  zphi0[k];
04730             t[C1INDEX(D010,D101)] = dxphi0[i] *  yphi1[j]  * dzphi0[k];
04731             t[C1INDEX(D010,D011)] =  xphi0[i] * dyphi1[j]  * dzphi0[k];
04732             t[C1INDEX(D010,D111)] = dxphi0[i] * dyphi1[j]  * dzphi0[k];
04733 
04734             t[C1INDEX(D001,D000)] =  xphi0[i] *  yphi0[j]  *  zphi1[k];
04735             t[C1INDEX(D001,D100)] = dxphi0[i] *  yphi0[j]  *  zphi1[k];
04736             t[C1INDEX(D001,D010)] =  xphi0[i] * dyphi0[j]  *  zphi1[k];
04737             t[C1INDEX(D001,D001)] =  xphi0[i] *  yphi0[j]  * dzphi1[k];
04738             t[C1INDEX(D001,D110)] = dxphi0[i] * dyphi0[j]  *  zphi1[k];
04739             t[C1INDEX(D001,D101)] = dxphi0[i] *  yphi0[j]  * dzphi1[k];
04740             t[C1INDEX(D001,D011)] =  xphi0[i] * dyphi0[j]  * dzphi1[k];
04741             t[C1INDEX(D001,D111)] = dxphi0[i] * dyphi0[j]  * dzphi1[k];
04742 
04743             t[C1INDEX(D110,D000)] =  xphi1[i] *  yphi1[j]  *  zphi0[k];
04744             t[C1INDEX(D110,D100)] = dxphi1[i] *  yphi1[j]  *  zphi0[k];
04745             t[C1INDEX(D110,D010)] =  xphi1[i] * dyphi1[j]  *  zphi0[k];
04746             t[C1INDEX(D110,D001)] =  xphi1[i] *  yphi1[j]  * dzphi0[k];
04747             t[C1INDEX(D110,D110)] = dxphi1[i] * dyphi1[j]  *  zphi0[k];
04748             t[C1INDEX(D110,D101)] = dxphi1[i] *  yphi1[j]  * dzphi0[k];
04749             t[C1INDEX(D110,D011)] =  xphi1[i] * dyphi1[j]  * dzphi0[k];
04750             t[C1INDEX(D110,D111)] = dxphi1[i] * dyphi1[j]  * dzphi0[k];
04751 
04752             t[C1INDEX(D101,D000)] =  xphi1[i] *  yphi0[j]  *  zphi1[k];
04753             t[C1INDEX(D101,D100)] = dxphi1[i] *  yphi0[j]  *  zphi1[k];
04754             t[C1INDEX(D101,D010)] =  xphi1[i] * dyphi0[j]  *  zphi1[k];
04755             t[C1INDEX(D101,D001)] =  xphi1[i] *  yphi0[j]  * dzphi1[k];
04756             t[C1INDEX(D101,D110)] = dxphi1[i] * dyphi0[j]  *  zphi1[k];
04757             t[C1INDEX(D101,D101)] = dxphi1[i] *  yphi0[j]  * dzphi1[k];
04758             t[C1INDEX(D101,D011)] =  xphi1[i] * dyphi0[j]  * dzphi1[k];
04759             t[C1INDEX(D101,D111)] = dxphi1[i] * dyphi0[j]  * dzphi1[k];
04760 
04761             t[C1INDEX(D011,D000)] =  xphi0[i] *  yphi1[j]  *  zphi1[k];
04762             t[C1INDEX(D011,D100)] = dxphi0[i] *  yphi1[j]  *  zphi1[k];
04763             t[C1INDEX(D011,D010)] =  xphi0[i] * dyphi1[j]  *  zphi1[k];
04764             t[C1INDEX(D011,D001)] =  xphi0[i] *  yphi1[j]  * dzphi1[k];
04765             t[C1INDEX(D011,D110)] = dxphi0[i] * dyphi1[j]  *  zphi1[k];
04766             t[C1INDEX(D011,D101)] = dxphi0[i] *  yphi1[j]  * dzphi1[k];
04767             t[C1INDEX(D011,D011)] =  xphi0[i] * dyphi1[j]  * dzphi1[k];
04768             t[C1INDEX(D011,D111)] = dxphi0[i] * dyphi1[j]  * dzphi1[k];
04769 
04770             t[C1INDEX(D111,D000)] =  xphi1[i] *  yphi1[j]  *  zphi1[k];
04771             t[C1INDEX(D111,D100)] = dxphi1[i] *  yphi1[j]  *  zphi1[k];
04772             t[C1INDEX(D111,D010)] =  xphi1[i] * dyphi1[j]  *  zphi1[k];
04773             t[C1INDEX(D111,D001)] =  xphi1[i] *  yphi1[j]  * dzphi1[k];
04774             t[C1INDEX(D111,D110)] = dxphi1[i] * dyphi1[j]  *  zphi1[k];
04775             t[C1INDEX(D111,D101)] = dxphi1[i] *  yphi1[j]  * dzphi1[k];
04776             t[C1INDEX(D111,D011)] =  xphi1[i] * dyphi1[j]  * dzphi1[k];
04777             t[C1INDEX(D111,D111)] = dxphi1[i] * dyphi1[j]  * dzphi1[k];
04778           }
04779         }
04780       } // end loops over restriction stencil matrices
04781 
04782       // loop over prolongation stencil matrices
04783       // prolongation stencil matrices are the transpose of restriction
04784       for (k = -NR;  k <= NR;  k++) {
04785         for (j = -NR;  j <= NR;  j++) {
04786           for (i = -NR;  i <= NR;  i++) {
04787             Float *t = map.gres_c1hermite[level](i+NR,j+NR,k+NR).melem;
04788             Float *tt = map.gpro_c1hermite[level](i+NR,j+NR,k+NR).melem;
04789 
04790             tt[C1INDEX(D000,D000)] = t[C1INDEX(D000,D000)];
04791             tt[C1INDEX(D000,D100)] = t[C1INDEX(D100,D000)];
04792             tt[C1INDEX(D000,D010)] = t[C1INDEX(D010,D000)];
04793             tt[C1INDEX(D000,D001)] = t[C1INDEX(D001,D000)];
04794             tt[C1INDEX(D000,D110)] = t[C1INDEX(D110,D000)];
04795             tt[C1INDEX(D000,D101)] = t[C1INDEX(D101,D000)];
04796             tt[C1INDEX(D000,D011)] = t[C1INDEX(D011,D000)];
04797             tt[C1INDEX(D000,D111)] = t[C1INDEX(D111,D000)];
04798 
04799             tt[C1INDEX(D100,D000)] = t[C1INDEX(D000,D100)];
04800             tt[C1INDEX(D100,D100)] = t[C1INDEX(D100,D100)];
04801             tt[C1INDEX(D100,D010)] = t[C1INDEX(D010,D100)];
04802             tt[C1INDEX(D100,D001)] = t[C1INDEX(D001,D100)];
04803             tt[C1INDEX(D100,D110)] = t[C1INDEX(D110,D100)];
04804             tt[C1INDEX(D100,D101)] = t[C1INDEX(D101,D100)];
04805             tt[C1INDEX(D100,D011)] = t[C1INDEX(D011,D100)];
04806             tt[C1INDEX(D100,D111)] = t[C1INDEX(D111,D100)];
04807 
04808             tt[C1INDEX(D010,D000)] = t[C1INDEX(D000,D010)];
04809             tt[C1INDEX(D010,D100)] = t[C1INDEX(D100,D010)];
04810             tt[C1INDEX(D010,D010)] = t[C1INDEX(D010,D010)];
04811             tt[C1INDEX(D010,D001)] = t[C1INDEX(D001,D010)];
04812             tt[C1INDEX(D010,D110)] = t[C1INDEX(D110,D010)];
04813             tt[C1INDEX(D010,D101)] = t[C1INDEX(D101,D010)];
04814             tt[C1INDEX(D010,D011)] = t[C1INDEX(D011,D010)];
04815             tt[C1INDEX(D010,D111)] = t[C1INDEX(D111,D010)];
04816 
04817             tt[C1INDEX(D001,D000)] = t[C1INDEX(D000,D001)];
04818             tt[C1INDEX(D001,D100)] = t[C1INDEX(D100,D001)];
04819             tt[C1INDEX(D001,D010)] = t[C1INDEX(D010,D001)];
04820             tt[C1INDEX(D001,D001)] = t[C1INDEX(D001,D001)];
04821             tt[C1INDEX(D001,D110)] = t[C1INDEX(D110,D001)];
04822             tt[C1INDEX(D001,D101)] = t[C1INDEX(D101,D001)];
04823             tt[C1INDEX(D001,D011)] = t[C1INDEX(D011,D001)];
04824             tt[C1INDEX(D001,D111)] = t[C1INDEX(D111,D001)];
04825 
04826             tt[C1INDEX(D110,D000)] = t[C1INDEX(D000,D110)];
04827             tt[C1INDEX(D110,D100)] = t[C1INDEX(D100,D110)];
04828             tt[C1INDEX(D110,D010)] = t[C1INDEX(D010,D110)];
04829             tt[C1INDEX(D110,D001)] = t[C1INDEX(D001,D110)];
04830             tt[C1INDEX(D110,D110)] = t[C1INDEX(D110,D110)];
04831             tt[C1INDEX(D110,D101)] = t[C1INDEX(D101,D110)];
04832             tt[C1INDEX(D110,D011)] = t[C1INDEX(D011,D110)];
04833             tt[C1INDEX(D110,D111)] = t[C1INDEX(D111,D110)];
04834 
04835             tt[C1INDEX(D101,D000)] = t[C1INDEX(D000,D101)];
04836             tt[C1INDEX(D101,D100)] = t[C1INDEX(D100,D101)];
04837             tt[C1INDEX(D101,D010)] = t[C1INDEX(D010,D101)];
04838             tt[C1INDEX(D101,D001)] = t[C1INDEX(D001,D101)];
04839             tt[C1INDEX(D101,D110)] = t[C1INDEX(D110,D101)];
04840             tt[C1INDEX(D101,D101)] = t[C1INDEX(D101,D101)];
04841             tt[C1INDEX(D101,D011)] = t[C1INDEX(D011,D101)];
04842             tt[C1INDEX(D101,D111)] = t[C1INDEX(D111,D101)];
04843 
04844             tt[C1INDEX(D011,D000)] = t[C1INDEX(D000,D011)];
04845             tt[C1INDEX(D011,D100)] = t[C1INDEX(D100,D011)];
04846             tt[C1INDEX(D011,D010)] = t[C1INDEX(D010,D011)];
04847             tt[C1INDEX(D011,D001)] = t[C1INDEX(D001,D011)];
04848             tt[C1INDEX(D011,D110)] = t[C1INDEX(D110,D011)];
04849             tt[C1INDEX(D011,D101)] = t[C1INDEX(D101,D011)];
04850             tt[C1INDEX(D011,D011)] = t[C1INDEX(D011,D011)];
04851             tt[C1INDEX(D011,D111)] = t[C1INDEX(D111,D011)];
04852 
04853             tt[C1INDEX(D111,D000)] = t[C1INDEX(D000,D111)];
04854             tt[C1INDEX(D111,D100)] = t[C1INDEX(D100,D111)];
04855             tt[C1INDEX(D111,D010)] = t[C1INDEX(D010,D111)];
04856             tt[C1INDEX(D111,D001)] = t[C1INDEX(D001,D111)];
04857             tt[C1INDEX(D111,D110)] = t[C1INDEX(D110,D111)];
04858             tt[C1INDEX(D111,D101)] = t[C1INDEX(D101,D111)];
04859             tt[C1INDEX(D111,D011)] = t[C1INDEX(D011,D111)];
04860             tt[C1INDEX(D111,D111)] = t[C1INDEX(D111,D111)];
04861           }
04862         }
04863       } // end loops over prolongation stencil matrices
04864 
04865     } // end loop over levels
04866 
04867   } // end if C1HERMITE
04868 
04869   // calculate self energy factor for splitting
04870   BigReal gs=0, d=0;
04871   splitting(gs, d, 0, split);
04872   gzero = gs * a_p;
04873 
04874   if (CkMyPe() == 0) {
04875     iout << iINFO << "MSM finished calculating stencils\n" << endi;
04876   }
04877 
04878   // allocate map for patches
04879   PatchMap *pm = PatchMap::Object();
04880   int numpatches = pm->numPatches();
04881   map.patchList.resize(numpatches);
04882 #ifdef DEBUG_MSM_VERBOSE
04883   printf("numPatches = %d\n", numpatches);
04884 #endif
04885 
04886   // allocate map for blocks for each grid level
04887   map.blockLevel.resize(nlevels);
04888   map.bsx.resize(nlevels);
04889   map.bsy.resize(nlevels);
04890   map.bsz.resize(nlevels);
04891 #ifdef MSM_FOLD_FACTOR
04892   map.foldfactor.resize(nlevels);
04893 #endif
04894   for (level = 0;  level < nlevels;  level++) {
04895     msm::IndexRange& g = map.gridrange[level];
04896     msm::Grid<msm::BlockDiagram>& b = map.blockLevel[level];
04897     int gia = g.ia();
04898     int gni = g.ni();
04899     int gja = g.ja();
04900     int gnj = g.nj();
04901     int gka = g.ka();
04902     int gnk = g.nk();
04903     map.bsx[level] = bsx;
04904     map.bsy[level] = bsy;
04905     map.bsz[level] = bsz;
04906     if (/* map.bsx[level] < gni ||
04907         map.bsy[level] < gnj ||
04908         map.bsz[level] < gnk */ 1) {
04909       // make sure that block sizes divide evenly into periodic dimensions
04910       if (ispx) setup_periodic_blocksize(map.bsx[level], gni);
04911       if (ispy) setup_periodic_blocksize(map.bsy[level], gnj);
04912       if (ispz) setup_periodic_blocksize(map.bsz[level], gnk);
04913 #ifdef MSM_DEBUG_VERBOSE
04914       if (CkMyPe() == 0) {
04915         printf("level = %d\n  map.bs* = %d %d %d  gn* = %d %d %d\n",
04916             level, map.bsx[level], map.bsy[level], map.bsz[level],gni,gnj,gnk);
04917       }
04918 #endif
04919       // subdivide grid into multiple blocks
04920       //   == ceil(gni / bsx), etc.
04921       int bni = (gni / map.bsx[level]) + (gni % map.bsx[level] != 0);
04922       int bnj = (gnj / map.bsy[level]) + (gnj % map.bsy[level] != 0);
04923       int bnk = (gnk / map.bsz[level]) + (gnk % map.bsz[level] != 0);
04924 #ifdef MSM_FOLD_FACTOR
04925       if (/* level > 2 && */ (bni == 1 || bnj == 1 || bnk == 1)) {
04926         map.foldfactor[level].set(bsx / gni, bsy / gnj, bsz / gnk);
04927 #if 0
04928         if (CkMyPe() == 0) {
04929           printf("Setting MSM FoldFactor level %d:  %d %d %d\n",
04930               level, bsx / gni, bsy / gnj, bsz / gnk);
04931         }
04932 #endif
04933       }
04934 #endif
04935       b.set(0, bni, 0, bnj, 0, bnk);
04936       for (k = 0;  k < bnk;  k++) {
04937         for (j = 0;  j < bnj;  j++) {
04938           for (i = 0;  i < bni;  i++) {
04939             b(i,j,k).reset();
04940             int ia = gia + i*map.bsx[level];
04941             int ib = ia + map.bsx[level] - 1;
04942             int ja = gja + j*map.bsy[level];
04943             int jb = ja + map.bsy[level] - 1;
04944             int ka = gka + k*map.bsz[level];
04945             int kb = ka + map.bsz[level] - 1;
04946             if (ib >= gia + gni) ib = gia + gni - 1;
04947             if (jb >= gja + gnj) jb = gja + gnj - 1;
04948             if (kb >= gka + gnk) kb = gka + gnk - 1;
04949             b(i,j,k).nrange.setbounds(ia, ib, ja, jb, ka, kb);
04950           }
04951         }
04952       }
04953     }
04954     /*
04955     else {
04956       // make entire grid into single block
04957       b.set(0, 1, 0, 1, 0, 1);
04958       b(0,0,0).reset();
04959       b(0,0,0).nrange.set(gia, gni, gja, gnj, gka, gnk);
04960       // set correct block dimensions
04961       map.bsx[level] = gni;
04962       map.bsy[level] = gnj;
04963       map.bsz[level] = gnk;
04964     }
04965     */
04966   }
04967   //CkExit();
04968 #ifdef DEBUG_MSM_VERBOSE
04969   printf("Done allocating map for grid levels\n");
04970   printf("Grid level decomposition:\n");
04971   for (level = 0;  level < nlevels;  level++) {
04972     msm::Grid<msm::BlockDiagram>& b = map.blockLevel[level];
04973     int bia = b.ia();
04974     int bib = b.ib();
04975     int bja = b.ja();
04976     int bjb = b.jb();
04977     int bka = b.ka();
04978     int bkb = b.kb();
04979     for (k = bka;  k <= bkb;  k++) {
04980       for (j = bja;  j <= bjb;  j++) {
04981         for (i = bia;  i <= bib;  i++) {
04982           int ia = b(i,j,k).nrange.ia();
04983           int ib = b(i,j,k).nrange.ib();
04984           int ja = b(i,j,k).nrange.ja();
04985           int jb = b(i,j,k).nrange.jb();
04986           int ka = b(i,j,k).nrange.ka();
04987           int kb = b(i,j,k).nrange.kb();
04988           printf("level=%d  id=%d %d %d  [%d..%d] x [%d..%d] x [%d..%d]"
04989               " --> %d\n",
04990               level, i, j, k, ia, ib, ja, jb, ka, kb,
04991               b(i,j,k).nrange.nn());
04992         }
04993       }
04994     }
04995   }
04996 #endif
04997   if (CkMyPe() == 0) {
04998     iout << iINFO << "MSM finished creating map for grid levels\n" << endi;
04999   }
05000 
05001   initialize2();
05002 } // ComputeMsmMgr::initialize()

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, MsmGridCutoffProxyMsg::put(), MsmC1HermiteGridCutoffProxyMsg::put(), MsmBlockProxyMsg::put(), MsmC1HermiteBlockProxyMsg::put(), and msm::Array< T >::resize().

05794                                       {
05795   int i, j, k, n, level;
05796 
05797   if (CkMyPe() == 0) {
05798 
05799     // on PE 0, create 3D chare array of MsmBlock for each level;
05800     // broadcast this array of proxies to the rest of the group
05801     if (approx == C1HERMITE) {
05802       msmC1HermiteBlock.resize(nlevels);
05803     }
05804     else {
05805       msmBlock.resize(nlevels);
05806     }
05807     for (level = 0;  level < nlevels;  level++) {
05808       int ni = map.blockLevel[level].ni();
05809       int nj = map.blockLevel[level].nj();
05810       int nk = map.blockLevel[level].nk();
05811 #ifdef MSM_NODE_MAPPING
05812       CkPrintf("Using MsmBlockMap for level %d\n", level);
05813       CProxy_MsmBlockMap blockMap = CProxy_MsmBlockMap::ckNew(level);
05814       CkArrayOptions opts(ni, nj, nk);
05815       opts.setMap(blockMap);
05816       if (approx == C1HERMITE) {
05817         msmC1HermiteBlock[level] =
05818           CProxy_MsmC1HermiteBlock::ckNew(level, opts);
05819       }
05820       else {
05821         msmBlock[level] = CProxy_MsmBlock::ckNew(level, opts);
05822       }
05823 #else
05824       if (approx == C1HERMITE) {
05825         msmC1HermiteBlock[level] =
05826           CProxy_MsmC1HermiteBlock::ckNew(level, ni, nj, nk);
05827       }
05828       else {
05829         msmBlock[level] = CProxy_MsmBlock::ckNew(level, ni, nj, nk);
05830       }
05831 #endif
05832 #ifdef DEBUG_MSM_VERBOSE
05833       printf("Create MsmBlock[%d] 3D chare array ( %d x %d x %d )\n",
05834           level, ni, nj, nk);
05835 #endif
05836       char msg[128];
05837       int nijk = ni * nj * nk;
05838       sprintf(msg, "MSM grid level %d decomposed into %d block%s"
05839           " ( %d x %d x %d )\n",
05840           level, nijk, (nijk==1 ? "" : "s"), ni, nj, nk);
05841       iout << iINFO << msg;
05842     }
05843     if (approx == C1HERMITE) {
05844       MsmC1HermiteBlockProxyMsg *msg = new MsmC1HermiteBlockProxyMsg;
05845       msg->put(msmC1HermiteBlock);
05846       msmProxy.recvMsmC1HermiteBlockProxy(msg);  // broadcast
05847     }
05848     else {
05849       MsmBlockProxyMsg *msg = new MsmBlockProxyMsg;
05850       msg->put(msmBlock);
05851       msmProxy.recvMsmBlockProxy(msg);  // broadcast
05852     }
05853 
05854 #ifdef MSM_GRID_CUTOFF_DECOMP
05855     // on PE 0, create 1D chare array of MsmGridCutoff
05856     // broadcast this array proxy to the rest of the group
05857 #ifdef MSM_NODE_MAPPING
05858     CkPrintf("Using MsmGridCutoffMap\n");
05859     CProxy_MsmGridCutoffMap gcutMap = CProxy_MsmGridCutoffMap::ckNew();
05860     CkArrayOptions optsgcut(numGridCutoff);
05861     optsgcut.setMap(gcutMap);
05862     if (approx == C1HERMITE) {
05863       msmC1HermiteGridCutoff = CProxy_MsmC1HermiteGridCutoff::ckNew(optsgcut);
05864     }
05865     else {
05866       msmGridCutoff = CProxy_MsmGridCutoff::ckNew(optsgcut);
05867     }
05868 #else
05869     if (approx == C1HERMITE) {
05870       msmC1HermiteGridCutoff =
05871         CProxy_MsmC1HermiteGridCutoff::ckNew(numGridCutoff);
05872     }
05873     else {
05874       msmGridCutoff = CProxy_MsmGridCutoff::ckNew(numGridCutoff);
05875     }
05876 #endif
05877     if (approx == C1HERMITE) {
05878       MsmC1HermiteGridCutoffProxyMsg *gcmsg =
05879         new MsmC1HermiteGridCutoffProxyMsg;
05880       gcmsg->put(&msmC1HermiteGridCutoff);
05881       msmProxy.recvMsmC1HermiteGridCutoffProxy(gcmsg);
05882     }
05883     else {
05884       MsmGridCutoffProxyMsg *gcmsg = new MsmGridCutoffProxyMsg;
05885       gcmsg->put(&msmGridCutoff);
05886       msmProxy.recvMsmGridCutoffProxy(gcmsg);
05887     }
05888 
05889     // XXX PE 0 initializes each MsmGridCutoff
05890     // one-to-many
05891     // for M length chare array, better for each PE to initialize M/P?
05892     for (level = 0;  level < nlevels;  level++) { // for all levels
05893       msm::Grid<msm::BlockDiagram>& b = map.blockLevel[level];
05894       int bni = b.ni();
05895       int bnj = b.nj();
05896       int bnk = b.nk();
05897       for (k = 0;  k < bnk;  k++) { // for all blocks
05898         for (j = 0;  j < bnj;  j++) {
05899           for (i = 0;  i < bni;  i++) {
05900             // source for charges
05901             msm::BlockIndex bi = msm::BlockIndex(level, msm::Ivec(i,j,k));
05902             int numSendAcross = b(i,j,k).sendAcross.len();
05903             ASSERT( numSendAcross == b(i,j,k).indexGridCutoff.len() );
05904             // for this source, loop over destinations for potentials
05905             for (n = 0;  n < numSendAcross;  n++) {
05906               msm::BlockSend &bs = b(i,j,k).sendAcross[n];
05907               int index = b(i,j,k).indexGridCutoff[n];
05908               MsmGridCutoffInitMsg *bsmsg = new MsmGridCutoffInitMsg(bi, bs);
05909               if (approx == C1HERMITE) {
05910                 msmC1HermiteGridCutoff[index].setup(bsmsg);
05911               }
05912               else {
05913                 msmGridCutoff[index].setup(bsmsg);
05914               }
05915             } // traverse sendAcross, indexGridCutoff arrays
05916 
05917           }
05918         }
05919       } // end for all blocks
05920 
05921     } // end for all levels
05922 
05923     iout << iINFO << "MSM grid cutoff calculation decomposed into "
05924       << numGridCutoff << " work objects\n";
05925 #endif
05926     iout << endi;
05927   }
05928 
05929 #ifdef DEBUG_MSM_VERBOSE
05930   printf("end of initialization\n");
05931 #endif
05932 } // ComputeMsmMgr::initialize_create()

void ComputeMsmMgr::initVirialContrib (  )  [inline]

Definition at line 529 of file ComputeMsm.C.

References cntVirialContrib, and gvsum.

Referenced by doneVirialContrib().

00529                            {
00530     gvsum.reset(0);
00531     cntVirialContrib = 0;
00532   }

msm::Map& ComputeMsmMgr::mapData (  )  [inline]

Definition at line 447 of file ComputeMsm.C.

References map.

Referenced by ComputeMsm::doWork(), MsmBlockKernel< Vtype, Mtype >::MsmBlockKernel(), msm::PatchData::PatchData(), and ComputeMsm::saveResults().

00447 { return map; }

static void ComputeMsmMgr::ndsplitting ( BigReal  pg[],
BigReal  s,
int  n,
int  _split 
) [inline, static]

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

01266                                                                       {
01267     int k = 0;
01268     if (k == n) return;
01269     if (s <= 1) {
01270       // compute derivatives of smoothed part
01271       switch (_split) {
01272         case TAYLOR2:
01273           pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8));
01274           if (k == n) break;
01275           pg[k++] = -1./2 + (s-1)*(3./4);
01276           if (k == n) break;
01277           pg[k++] = 3./4;
01278           break;
01279         case TAYLOR3:
01280           pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16)));
01281           if (k == n) break;
01282           pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16));
01283           if (k == n) break;
01284           pg[k++] = 3./4 + (s-1)*(-15./8);
01285           if (k == n) break;
01286           pg[k++] = -15./8;
01287           break;
01288         case TAYLOR4:
01289           pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
01290                   + (s-1)*(35./128))));
01291           if (k == n) break;
01292           pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32)));
01293           if (k == n) break;
01294           pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32));
01295           if (k == n) break;
01296           pg[k++] = -15./8 + (s-1)*(105./16);
01297           if (k == n) break;
01298           pg[k++] = 105./16;
01299           break;
01300         case TAYLOR5:
01301           pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
01302                   + (s-1)*(35./128 + (s-1)*(-63./256)))));
01303           if (k == n) break;
01304           pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
01305                   + (s-1)*(-315./256))));
01306           if (k == n) break;
01307           pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64)));
01308           if (k == n) break;
01309           pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64));
01310           if (k == n) break;
01311           pg[k++] = 105./16 + (s-1)*(-945./32);
01312           if (k == n) break;
01313           pg[k++] = -945./32;
01314           break;
01315         case TAYLOR6:
01316           pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
01317                   + (s-1)*(35./128 + (s-1)*(-63./256 + (s-1)*(231./1024))))));
01318           if (k == n) break;
01319           pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
01320                   + (s-1)*(-315./256 + (s-1)*(693./512)))));
01321           if (k == n) break;
01322           pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64
01323                   + (s-1)*(3465./512))));
01324           if (k == n) break;
01325           pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64
01326                 + (s-1)*(3465./128)));
01327           if (k == n) break;
01328           pg[k++] = 105./16 + (s-1)*(-945./32 + (s-1)*(10395./128));
01329           if (k == n) break;
01330           pg[k++] = -945./32 + (s-1)*(10395./64);
01331           if (k == n) break;
01332           pg[k++] = 10395./64;
01333           break;
01334         case TAYLOR7:
01335           pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
01336                   + (s-1)*(35./128 + (s-1)*(-63./256
01337                       + (s-1)*(231./1024 + (s-1)*(-429./2048)))))));
01338           if (k == n) break;
01339           pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
01340                   + (s-1)*(-315./256 + (s-1)*(693./512
01341                       + (s-1)*(-3003./2048))))));
01342           if (k == n) break;
01343           pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64
01344                   + (s-1)*(3465./512 + (s-1)*(-9009./1024)))));
01345           if (k == n) break;
01346           pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64 + (s-1)*(3465./128
01347                   + (s-1)*(-45045./1024))));
01348           if (k == n) break;
01349           pg[k++] = 105./16 + (s-1)*(-945./32 + (s-1)*(10395./128
01350                 + (s-1)*(-45045./256)));
01351           if (k == n) break;
01352           pg[k++] = -945./32 + (s-1)*(10395./64 + (s-1)*(-135135./256));
01353           if (k == n) break;
01354           pg[k++] = 10395./64 + (s-1)*(-135135./128);
01355           if (k == n) break;
01356           pg[k++] = -135135./128;
01357           break;
01358         case TAYLOR8:
01359           pg[k++] = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
01360                   + (s-1)*(35./128 + (s-1)*(-63./256
01361                       + (s-1)*(231./1024 + (s-1)*(-429./2048
01362                           + (s-1)*(6435./32768))))))));
01363           if (k == n) break;
01364           pg[k++] = -1./2 + (s-1)*(3./4 + (s-1)*(-15./16 + (s-1)*(35./32
01365                   + (s-1)*(-315./256 + (s-1)*(693./512
01366                       + (s-1)*(-3003./2048 + (s-1)*(6435./4096)))))));
01367           if (k == n) break;
01368           pg[k++] = 3./4 + (s-1)*(-15./8 + (s-1)*(105./32 + (s-1)*(-315./64
01369                   + (s-1)*(3465./512 + (s-1)*(-9009./1024
01370                       + (s-1)*(45045./4096))))));
01371           if (k == n) break;
01372           pg[k++] = -15./8 + (s-1)*(105./16 + (s-1)*(-945./64 + (s-1)*(3465./128
01373                   + (s-1)*(-45045./1024 + (s-1)*(135135./2048)))));
01374           if (k == n) break;
01375           pg[k++] = 105./16 + (s-1)*(-945./32 + (s-1)*(10395./128
01376                 + (s-1)*(-45045./256 + (s-1)*(675675./2048))));
01377           if (k == n) break;
01378           pg[k++] = -945./32 + (s-1)*(10395./64 + (s-1)*(-135135./256
01379                 + (s-1)*(675675./512)));
01380           if (k == n) break;
01381           pg[k++] = 10395./64 + (s-1)*(-135135./128 + (s-1)*(2027025./512));
01382           if (k == n) break;
01383           pg[k++] = -135135./128 + (s-1)*(2027025./256);
01384           if (k == n) break;
01385           pg[k++] = 2027025./256;
01386           break;
01387         default:
01388           NAMD_die("Unknown MSM splitting.");
01389       }
01390     } // if (s <= 1)
01391     else { // (s > 1)
01392       // compute derivatives of s^(-1/2)
01393       const BigReal s_1 = 1./s;
01394       BigReal s_p = sqrt(s_1);
01395       BigReal p = -0.5;
01396       BigReal _c = 1;
01397       pg[k++] = _c * s_p;
01398       while (k < n) {
01399         s_p *= s_1;
01400         _c *= p;
01401         p -= 1;
01402         pg[k++] = _c * s_p;
01403       }
01404     } // else (s > 1)
01405     // higher derivatives are zero
01406     while (k < n) pg[k++] = 0;
01407   } // ndsplitting()

int ComputeMsmMgr::numLevels (  )  const [inline]

Definition at line 449 of file ComputeMsm.C.

References nlevels.

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

00449 { return nlevels; }

msm::PatchPtrArray& ComputeMsmMgr::patchPtrArray (  )  [inline]

Definition at line 445 of file ComputeMsm.C.

References patchPtr.

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

00445 { return patchPtr; }

void ComputeMsmMgr::recvMsmBlockProxy ( MsmBlockProxyMsg  ) 

Definition at line 5935 of file ComputeMsm.C.

References MsmBlockProxyMsg::get(), and msmBlock.

05936 {
05937   msg->get(msmBlock);
05938   delete(msg);
05939 }

void ComputeMsmMgr::recvMsmC1HermiteBlockProxy ( MsmC1HermiteBlockProxyMsg  ) 

Definition at line 5947 of file ComputeMsm.C.

References MsmC1HermiteBlockProxyMsg::get(), and msmC1HermiteBlock.

05950 {
05951   msg->get(msmC1HermiteBlock);
05952   delete(msg);
05953 }

void ComputeMsmMgr::recvMsmC1HermiteGridCutoffProxy ( MsmC1HermiteGridCutoffProxyMsg  ) 

Definition at line 5955 of file ComputeMsm.C.

References MsmC1HermiteGridCutoffProxyMsg::get(), and msmC1HermiteGridCutoff.

05958 {
05959   msg->get(&msmC1HermiteGridCutoff);
05960   delete(msg);
05961 }

void ComputeMsmMgr::recvMsmGridCutoffProxy ( MsmGridCutoffProxyMsg  ) 

Definition at line 5941 of file ComputeMsm.C.

References MsmGridCutoffProxyMsg::get(), and msmGridCutoff.

05942 {
05943   msg->get(&msmGridCutoff);
05944   delete(msg);
05945 }

void ComputeMsmMgr::setCompute ( ComputeMsm c  )  [inline]

Definition at line 443 of file ComputeMsm.C.

References c, and msmCompute.

00443 { msmCompute = c;  c->setMgr(this); } // local

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

03891 {
03892   ASSERT(gridspacing > 0);
03893   if (isperiodic) {
03894     const BigReal hmin = (4./5) * gridspacing;
03895     const BigReal hmax = 1.5 * hmin;
03896     hh = len;
03897     nn = 1;  // start with one grid point across length
03898     while (hh >= hmax) {
03899       hh *= 0.5;  // halve spacing and double grid points
03900       nn <<= 1;
03901     }
03902     if (hh < hmin) {
03903       if (nn < 4) {
03904         NAMD_die("Basis vector is too short or MSM grid spacing is too large");
03905       }
03906       hh *= (4./3);  // scale hh by 4/3 and nn by 3/4
03907       nn >>= 2;
03908       nn *= 3;
03909     }
03910     // now we have:  hmin <= h < hmax,
03911     // where nn is a power of 2 times no more than one power of 3
03912     ia = 0;
03913     ib = nn-1;
03914   }
03915   else {
03916     hh = gridspacing;
03917     // Instead of "nn = (int) ceil(len / hh);"
03918     // len is divisible by hh, up to roundoff error, so round to closest nn
03919     nn = (int) floor(len/hh + 0.5);
03920     ia = -s_edge;
03921     ib = nn + s_edge;
03922   }
03923 } // ComputeMsmMgr::setup_hgrid_1d()

void ComputeMsmMgr::setup_periodic_blocksize ( int &  bsize,
int  n 
)

Definition at line 3928 of file ComputeMsm.C.

References NAMD_die().

Referenced by initialize().

03929 {
03930   if (n % bsize != 0) {
03931     // n is either 2^k or 3*2^k
03932     int newbsize = 1;
03933     if (n % 3 == 0) newbsize = 3;
03934     while (newbsize < bsize && newbsize < n) newbsize *= 2;
03935     if (bsize < newbsize) newbsize /= 2;
03936     if (n % newbsize != 0) {
03937       NAMD_die("MSM grid size for periodic dimensions must be "
03938           "a power of 2 times at most one power of 3");
03939     }
03940     bsize = newbsize;
03941   }
03942   return;
03943 }

static int ComputeMsmMgr::sign ( int  n  )  [inline, static]

Definition at line 452 of file ComputeMsm.C.

00452                                 {
00453     return (n < 0 ? -1 : (n > 0 ? 1 : 0));
00454   }

static void ComputeMsmMgr::splitting ( BigReal g,
BigReal dg,
BigReal  r_a,
int  _split 
) [inline, static]

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

00667                                                                           {
00668     BigReal s = r_a * r_a;  // s = (r/a)^2, assuming 0 <= s <= 1
00669     switch (_split) {
00670       case TAYLOR2:
00671         g = 1 + (s-1)*(-1./2 + (s-1)*(3./8));
00672         dg = (2*r_a)*(-1./2 + (s-1)*(3./4));
00673         break;
00674       case TAYLOR3:
00675         g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16)));
00676         dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16)));
00677         break;
00678       case TAYLOR4:
00679         g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
00680                 + (s-1)*(35./128))));
00681         dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
00682                 + (s-1)*(35./32))));
00683         break;
00684       case TAYLOR5:
00685         g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
00686                 + (s-1)*(35./128 + (s-1)*(-63./256)))));
00687         dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
00688                 + (s-1)*(35./32 + (s-1)*(-315./256)))));
00689         break;
00690       case TAYLOR6:
00691         g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
00692                 + (s-1)*(35./128 + (s-1)*(-63./256
00693                     + (s-1)*(231./1024))))));
00694         dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
00695                 + (s-1)*(35./32 + (s-1)*(-315./256
00696                     + (s-1)*(693./512))))));
00697         break;
00698       case TAYLOR7:
00699         g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
00700             + (s-1)*(35./128 + (s-1)*(-63./256
00701                 + (s-1)*(231./1024 + (s-1)*(-429./2048)))))));
00702         dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
00703                 + (s-1)*(35./32 + (s-1)*(-315./256
00704                     + (s-1)*(693./512 + (s-1)*(-3003./2048)))))));
00705         break;
00706       case TAYLOR8:
00707         g = 1 + (s-1)*(-1./2 + (s-1)*(3./8 + (s-1)*(-5./16
00708                 + (s-1)*(35./128 + (s-1)*(-63./256
00709                     + (s-1)*(231./1024 + (s-1)*(-429./2048
00710                         + (s-1)*(6435./32768))))))));
00711         dg = (2*r_a)*(-1./2 + (s-1)*(3./4 + (s-1)*(-15./16
00712                 + (s-1)*(35./32 + (s-1)*(-315./256
00713                     + (s-1)*(693./512 + (s-1)*(-3003./2048
00714                         + (s-1)*(6435./4096))))))));
00715         break;
00716       case TAYLOR2_DISP:
00717         g = 1 + (s-1)*(-3 + (s-1)*(6));
00718         dg = (2*r_a)*(-3 + (s-1)*(12));
00719         break;
00720       case TAYLOR3_DISP:
00721         g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10)));
00722         dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30)));
00723         break;
00724       case TAYLOR4_DISP:
00725         g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10 + (s-1)*(15))));
00726         dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30 + (s-1)*(60))));
00727         break;
00728       case TAYLOR5_DISP:
00729         g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
00730                 + (s-1)*(15 + (s-1)*(-21)))));
00731         dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
00732                 + (s-1)*(60 + (s-1)*(-105)))));
00733         break;
00734       case TAYLOR6_DISP:
00735         g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
00736                 + (s-1)*(15 + (s-1)*(-21 + (s-1)*(28))))));
00737         dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
00738                 + (s-1)*(60 + (s-1)*(-105 + (s-1)*(168))))));
00739         break;
00740       case TAYLOR7_DISP:
00741         g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
00742                 + (s-1)*(15 + (s-1)*(-21 + (s-1)*(28
00743                       + (s-1)*(-36)))))));
00744         dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
00745                 + (s-1)*(60 + (s-1)*(-105 + (s-1)*(168
00746                       + (s-1)*(-252)))))));
00747         break;
00748       case TAYLOR8_DISP:
00749         g = 1 + (s-1)*(-3 + (s-1)*(6 + (s-1)*(-10
00750                 + (s-1)*(15 + (s-1)*(-21 + (s-1)*(28
00751                       + (s-1)*(-36 + (s-1)*(45))))))));
00752         dg = (2*r_a)*(-3 + (s-1)*(12 + (s-1)*(-30
00753                 + (s-1)*(60 + (s-1)*(-105 + (s-1)*(168
00754                       + (s-1)*(-252 + (s-1)*(360))))))));
00755         break;
00756       default:
00757         NAMD_die("Unknown MSM splitting.");
00758     } // switch
00759   } // splitting()

void ComputeMsmMgr::stencil_1d ( Float  phi[],
Float  t 
) [inline]

Definition at line 761 of file ComputeMsm.C.

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

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

00761                                         {
00762     switch (approx) {
00763       case CUBIC:
00764         phi[0] = 0.5f * (1 - t) * (2 - t) * (2 - t);
00765         t--;
00766         phi[1] = (1 - t) * (1 + t - 1.5f * t * t);
00767         t--;
00768         phi[2] = (1 + t) * (1 - t - 1.5f * t * t);
00769         t--;
00770         phi[3] = 0.5f * (1 + t) * (2 + t) * (2 + t);
00771         break;
00772       case QUINTIC:
00773         phi[0] = (1.f/24) * (1-t) * (2-t) * (3-t) * (3-t) * (4-t);
00774         t--;
00775         phi[1] = (1-t) * (2-t) * (3-t) * ((1.f/6)
00776             + t * (0.375f - (5.f/24)*t));
00777         t--;
00778         phi[2] = (1-t*t) * (2-t) * (0.5f + t * (0.25f - (5.f/12)*t));
00779         t--;
00780         phi[3] = (1-t*t) * (2+t) * (0.5f - t * (0.25f + (5.f/12)*t));
00781         t--;
00782         phi[4] = (1+t) * (2+t) * (3+t) * ((1.f/6)
00783             - t * (0.375f + (5.f/24)*t));
00784         t--;
00785         phi[5] = (1.f/24) * (1+t) * (2+t) * (3+t) * (3+t) * (4+t);
00786         break;
00787       case QUINTIC2:
00788         phi[0] = (1.f/24) * (3-t) * (3-t) * (3-t) * (t-2) * (5*t-8);
00789         t--;
00790         phi[1] = (-1.f/24) * (2-t) * (t-1) * (-48+t*(153+t*(-114+t*25)));
00791         t--;
00792         phi[2] = (1.f/12) * (1-t) * (12+t*(12+t*(-3+t*(-38+t*25))));
00793         t--;
00794         phi[3] = (1.f/12) * (1+t) * (12+t*(-12+t*(-3+t*(38+t*25))));
00795         t--;
00796         phi[4] = (-1.f/24) * (2+t) * (t+1) * (48+t*(153+t*(114+t*25)));
00797         t--;
00798         phi[5] = (1.f/24) * (3+t) * (3+t) * (3+t) * (t+2) * (5*t+8);
00799         break;
00800       case SEPTIC:
00801         phi[0] = (-1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-4)*(t-5)*(t-6);
00802         t--;
00803         phi[1] = (1.f/720)*(t-1)*(t-2)*(t-3)*(t-4)*(t-5)*(-6+t*(-20+7*t));
00804         t--;
00805         phi[2] = (-1.f/240)*(t*t-1)*(t-2)*(t-3)*(t-4)*(-10+t*(-12+7*t));
00806         t--;
00807         phi[3] = (1.f/144)*(t*t-1)*(t*t-4)*(t-3)*(-12+t*(-4+7*t));
00808         t--;
00809         phi[4] = (-1.f/144)*(t*t-1)*(t*t-4)*(t+3)*(-12+t*(4+7*t));
00810         t--;
00811         phi[5] = (1.f/240)*(t*t-1)*(t+2)*(t+3)*(t+4)*(-10+t*(12+7*t));
00812         t--;
00813         phi[6] = (-1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(-6+t*(20+7*t));
00814         t--;
00815         phi[7] = (1.f/720)*(t+1)*(t+2)*(t+3)*(t+4)*(t+4)*(t+5)*(t+6);
00816         break;
00817       case SEPTIC3:
00818         phi[0] = (3632.f/5) + t*((-7456.f/5) + t*((58786.f/45) + t*(-633
00819                 + t*((26383.f/144) + t*((-22807.f/720) + t*((727.f/240)
00820                       + t*(-89.f/720)))))));
00821         t--;
00822         phi[1] = -440 + t*((25949.f/20) + t*((-117131.f/72) + t*((2247.f/2)
00823                 + t*((-66437.f/144) + t*((81109.f/720) + t*((-727.f/48)
00824                       + t*(623.f/720)))))));
00825         t--;
00826         phi[2] = (138.f/5) + t*((-8617.f/60) + t*((12873.f/40) + t*((-791.f/2)
00827                 + t*((4557.f/16) + t*((-9583.f/80) + t*((2181.f/80)
00828                       + t*(-623.f/240)))))));
00829         t--;
00830         phi[3] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((2569.f/144)
00831                 + t*((-727.f/48) + t*(623.f/144)))));
00832         t--;
00833         phi[4] = 1 + t*t*((-49.f/36) + t*t*((-959.f/144) + t*((-2569.f/144)
00834                 + t*((-727.f/48) + t*(-623.f/144)))));
00835         t--;
00836         phi[5] = (138.f/5) + t*((8617.f/60) + t*((12873.f/40) + t*((791.f/2)
00837                 + t*((4557.f/16) + t*((9583.f/80) + t*((2181.f/80)
00838                       + t*(623.f/240)))))));
00839         t--;
00840         phi[6] = -440 + t*((-25949.f/20) + t*((-117131.f/72) + t*((-2247.f/2)
00841                 + t*((-66437.f/144) + t*((-81109.f/720) + t*((-727.f/48)
00842                       + t*(-623.f/720)))))));
00843         t--;
00844         phi[7] = (3632.f/5) + t*((7456.f/5) + t*((58786.f/45) + t*(633
00845                 + t*((26383.f/144) + t*((22807.f/720) + t*((727.f/240)
00846                       + t*(89.f/720)))))));
00847         break;
00848       case NONIC:
00849         phi[0] = (-1.f/40320)*(t-8)*(t-7)*(t-6)*(t-5)*(t-5)*(t-4)*(t-3)*
00850           (t-2)*(t-1);
00851         t--;
00852         phi[1] = (1.f/40320)*(t-7)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*
00853           (-8+t*(-35+9*t));
00854         t--;
00855         phi[2] = (-1.f/10080)*(t-6)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*
00856           (-14+t*(-25+9*t));
00857         t--;
00858         phi[3] = (1.f/1440)*(t-5)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*
00859           (-6+t*(-5+3*t));
00860         t--;
00861         phi[4] = (-1.f/2880)*(t-4)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*
00862           (-20+t*(-5+9*t));
00863         t--;
00864         phi[5] = (1.f/2880)*(t-3)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*
00865           (-20+t*(5+9*t));
00866         t--;
00867         phi[6] = (-1.f/1440)*(t-2)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*
00868           (-6+t*(5+3*t));
00869         t--;
00870         phi[7] = (1.f/10080)*(t-1)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*
00871           (-14+t*(25+9*t));
00872         t--;
00873         phi[8] = (-1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+6)*(t+7)*
00874           (-8+t*(35+9*t));
00875         t--;
00876         phi[9] = (1.f/40320)*(t+1)*(t+2)*(t+3)*(t+4)*(t+5)*(t+5)*(t+6)*
00877           (t+7)*(t+8);
00878         break;
00879       case NONIC4:
00880       { // begin grouping to define local variables
00881         double Tphi[10];
00882         double T=t;
00883         Tphi[0] = 439375./7+T*(-64188125./504+T*(231125375./2016
00884               +T*(-17306975./288+T*(7761805./384+T*(-2895587./640
00885                     +T*(129391./192+T*(-259715./4032+T*(28909./8064
00886                           +T*(-3569./40320)))))))));
00887         T--;
00888         Tphi[1] = -56375+T*(8314091./56+T*(-49901303./288+T*(3763529./32
00889                 +T*(-19648027./384+T*(9469163./640+T*(-545977./192
00890                       +T*(156927./448+T*(-28909./1152
00891                           +T*(3569./4480)))))))));
00892         T--;
00893         Tphi[2] = 68776./7+T*(-1038011./28+T*(31157515./504+T*(-956669./16
00894                 +T*(3548009./96+T*(-2422263./160+T*(197255./48
00895                       +T*(-19959./28+T*(144545./2016
00896                           +T*(-3569./1120)))))))));
00897         T--;
00898         Tphi[3] = -154+T*(12757./12+T*(-230123./72+T*(264481./48
00899                 +T*(-576499./96+T*(686147./160+T*(-96277./48
00900                       +T*(14221./24+T*(-28909./288+T*(3569./480)))))))));
00901         T--;
00902         Tphi[4] = 1+T*T*(-205./144+T*T*(91./192+T*(-6181./320
00903                 +T*(6337./96+T*(-2745./32+T*(28909./576
00904                       +T*(-3569./320)))))));
00905         T--;
00906         Tphi[5] = 1+T*T*(-205./144+T*T*(91./192+T*(6181./320
00907                 +T*(6337./96+T*(2745./32+T*(28909./576
00908                       +T*(3569./320)))))));
00909         T--;
00910         Tphi[6] = -154+T*(-12757./12+T*(-230123./72+T*(-264481./48
00911                 +T*(-576499./96+T*(-686147./160+T*(-96277./48
00912                       +T*(-14221./24+T*(-28909./288+T*(-3569./480)))))))));
00913         T--;
00914         Tphi[7] = 68776./7+T*(1038011./28+T*(31157515./504+T*(956669./16
00915                 +T*(3548009./96+T*(2422263./160+T*(197255./48
00916                       +T*(19959./28+T*(144545./2016+T*(3569./1120)))))))));
00917         T--;
00918         Tphi[8] = -56375+T*(-8314091./56+T*(-49901303./288+T*(-3763529./32
00919                 +T*(-19648027./384+T*(-9469163./640+T*(-545977./192
00920                       +T*(-156927./448+T*(-28909./1152
00921                           +T*(-3569./4480)))))))));
00922         T--;
00923         Tphi[9] = 439375./7+T*(64188125./504+T*(231125375./2016
00924               +T*(17306975./288+T*(7761805./384+T*(2895587./640
00925                     +T*(129391./192+T*(259715./4032+T*(28909./8064
00926                           +T*(3569./40320)))))))));
00927         for (int i=0;  i < 10;  i++) {
00928           phi[i] = Float(Tphi[i]);
00929         }
00930       } // end grouping to define local variables
00931         break;
00932       default:
00933         NAMD_die("Unknown MSM approximation.");
00934     } // switch
00935   } // stencil_1d()

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

01244                                                                         {
01245     phi[0] = (1 - t) * (1 - t) * (1 + 2*t);
01246     psi[0] = h * t * (1 - t) * (1 - t);
01247     t--;
01248     phi[1] = (1 + t) * (1 + t) * (1 - 2*t);
01249     psi[1] = h * t * (1 + t) * (1 + t);
01250   }

void ComputeMsmMgr::subtractVirialContrib (  )  [inline]

Definition at line 536 of file ComputeMsm.C.

References numVirialContrib.

00536                                {
00537     numVirialContrib--;
00538   }

void ComputeMsmMgr::update ( CkQdMsg *   ) 

Definition at line 5963 of file ComputeMsm.C.

References approx, C1HERMITE, msmBlock, and msmC1HermiteBlock.

05964 {
05965 #ifdef DEBUG_MSM_VERBOSE
05966   printf("ComputeMsmMgr:  update() PE %d\n", CkMyPe());
05967 #endif
05968   delete msg;
05969 
05970   // have to setup sections AFTER initialization is finished
05971   if (CkMyPe() == 0) {
05972     for (int level = 0;  level < nlevels;  level++) {
05973       if (approx == C1HERMITE) {
05974         msmC1HermiteBlock[level].setupSections();
05975       }
05976       else {
05977         msmBlock[level].setupSections();
05978       }
05979     }
05980   }
05981 
05982   // XXX how do update for constant pressure simulation?
05983 }


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

Definition at line 604 of file ComputeMsm.C.

Referenced by addPotential(), msm::PatchData::anterpolation(), calcBlockWork(), compute(), d_stencil_1d(), msm::PatchData::init(), initialize(), initialize_create(), msm::PatchData::interpolation(), msm::PatchData::PatchData(), MsmBlockKernel< Vtype, Mtype >::prolongationKernel(), MsmBlockKernel< Vtype, Mtype >::restrictionKernel(), stencil_1d(), and update().

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

Definition at line 608 of file ComputeMsm.C.

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

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

Definition at line 599 of file ComputeMsm.C.

Referenced by msm::PatchData::anterpolationC1Hermite(), initialize(), and msm::PatchData::interpolationC1Hermite().

BigReal ComputeMsmMgr::hxlen_1

Definition at line 600 of file ComputeMsm.C.

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

BigReal ComputeMsmMgr::hylen

Definition at line 599 of file ComputeMsm.C.

Referenced by msm::PatchData::anterpolationC1Hermite(), initialize(), and msm::PatchData::interpolationC1Hermite().

BigReal ComputeMsmMgr::hylen_1

Definition at line 600 of file ComputeMsm.C.

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

BigReal ComputeMsmMgr::hzlen

Definition at line 599 of file ComputeMsm.C.

Referenced by msm::PatchData::anterpolationC1Hermite(), initialize(), and msm::PatchData::interpolationC1Hermite().

BigReal ComputeMsmMgr::hzlen_1

Definition at line 600 of file ComputeMsm.C.

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

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

Definition at line 592 of file ComputeMsm.C.

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

msm::Map ComputeMsmMgr::map

Definition at line 471 of file ComputeMsm.C.

Referenced by blockFlatIndex(), calcBlockWork(), calcGcutWork(), MsmGridCutoffKernel< C1Vector, C1Matrix >::init(), initialize(), initialize_create(), mapData(), MsmC1HermiteGridCutoff::setup(), MsmGridCutoff::setup(), and MsmGridCutoffKernel< C1Vector, C1Matrix >::setup().

msm::Array<CProxy_MsmBlock> ComputeMsmMgr::msmBlock

Definition at line 464 of file ComputeMsm.C.

Referenced by initialize_create(), recvMsmBlockProxy(), msm::PatchData::sendCharge(), MsmBlock::sendDownPotential(), MsmBlock::sendUpCharge(), and update().

msm::Array<CProxy_MsmC1HermiteBlock> ComputeMsmMgr::msmC1HermiteBlock

Definition at line 465 of file ComputeMsm.C.

Referenced by initialize_create(), recvMsmC1HermiteBlockProxy(), msm::PatchData::sendChargeC1Hermite(), MsmC1HermiteBlock::sendDownPotential(), MsmC1HermiteBlock::sendUpCharge(), and update().

CProxy_MsmC1HermiteGridCutoff ComputeMsmMgr::msmC1HermiteGridCutoff

Definition at line 468 of file ComputeMsm.C.

Referenced by initialize_create(), recvMsmC1HermiteGridCutoffProxy(), and MsmC1HermiteBlock::setupSections().

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 MsmC1HermiteBlock::gridCutoff(), MsmBlock::gridCutoff(), initialize_create(), recvMsmGridCutoffProxy(), and MsmBlock::setupSections().

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

Definition at line 606 of file ComputeMsm.C.

Referenced by MsmC1HermiteBlock::gridCutoff(), MsmBlock::gridCutoff(), initialize(), numLevels(), MsmC1HermiteBlock::sendDownPotential(), MsmBlock::sendDownPotential(), MsmC1HermiteBlock::sendPatch(), MsmBlock::sendPatch(), MsmC1HermiteBlock::sumReducedPotential(), and MsmBlock::sumReducedPotential().

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

Definition at line 622 of file ComputeMsm.C.

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

Vector ComputeMsmMgr::sglower

Definition at line 610 of file ComputeMsm.C.

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

BigReal ComputeMsmMgr::shx

Definition at line 613 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::shx_1

Definition at line 614 of file ComputeMsm.C.

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

BigReal ComputeMsmMgr::shy

Definition at line 613 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::shy_1

Definition at line 614 of file ComputeMsm.C.

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

BigReal ComputeMsmMgr::shz

Definition at line 613 of file ComputeMsm.C.

Referenced by initialize().

BigReal ComputeMsmMgr::shz_1

Definition at line 614 of file ComputeMsm.C.

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

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

Definition at line 618 of file ComputeMsm.C.

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

Float ComputeMsmMgr::srx_y

Definition at line 618 of file ComputeMsm.C.

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

Float ComputeMsmMgr::srx_z

Definition at line 618 of file ComputeMsm.C.

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

Float ComputeMsmMgr::sry_x

Definition at line 619 of file ComputeMsm.C.

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

Float ComputeMsmMgr::sry_y

Definition at line 619 of file ComputeMsm.C.

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

Float ComputeMsmMgr::sry_z

Definition at line 619 of file ComputeMsm.C.

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

Float ComputeMsmMgr::srz_x

Definition at line 620 of file ComputeMsm.C.

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

Float ComputeMsmMgr::srz_y

Definition at line 620 of file ComputeMsm.C.

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

Float ComputeMsmMgr::srz_z

Definition at line 620 of file ComputeMsm.C.

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

msm::Grid<Float> ComputeMsmMgr::subgrid

Definition at line 480 of file ComputeMsm.C.

Referenced by addPotential().

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:
Generated on Mon Sep 25 01:17:17 2017 for NAMD by  doxygen 1.4.7