OpenAtom  Version1.5a
CP_Rho_RHartExt.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 //////////////////////////////////////////////////////////////////////////////
3 //////////////////////////////////////////////////////////////////////////////
4 /** \file CP_Rho_RHartExt.C
5  *
6  * This is a description of the "life" of a CP_Rho_RHartExt object
7  *
8  * Fill in details here
9  */
10 //////////////////////////////////////////////////////////////////////////////
11 
12 #include "charm++.h"
13 #include <iostream>
14 #include <fstream>
15 #include <cmath>
16 
17 #include "debug_flags.h"
18 #include "utility/util.h"
19 #include "main/cpaimd.h"
20 #include "main/AtomsCache.h"
21 #include "fft_slab_ctrl/fftCacheSlab.h"
22 #include "main/eesCache.h"
23 #include "cp_state_ctrl/CP_State_Plane.h"
24 
25 #include "src_piny_physics_v1.0/include/class_defs/CP_OPERATIONS/class_cplocal.h"
26 #include "src_piny_physics_v1.0/include/class_defs/CP_OPERATIONS/class_cpxcfnctls.h"
27 
28 //////////////////////////////////////////////////////////////////////////////
29 
30 extern CkVec <CProxy_CP_Rho_RealSpacePlane> UrhoRealProxy;
31 extern CkVec <CProxy_AtomsCache> UatomsCacheProxy;
32 extern CkVec <CProxy_CP_Rho_RHartExt> UrhoRHartExtProxy;
33 extern CkVec <CProxy_CP_Rho_GHartExt> UrhoGHartExtProxy;
34 extern CkVec <CProxy_eesCache> UeesCacheProxy;
35 extern CkVec <CProxy_FFTcache> UfftCacheProxy;
36 
37 extern ComlibInstanceHandle commRHartGHartIns;
38 
39 extern int sizeX;
40 extern Config config;
41 
42 
43 //#define _CP_RHART_VERBOSE_
44 
45 /** \file CP_Rho_RHartExt.C
46  * @addtogroup Density
47  @{
48 */
49 //////////////////////////////////////////////////////////////////////////////
50 //////////////////////////////////////////////////////////////////////////////
51 //////////////////////////////////////////////////////////////////////////////
52 /**
53  * This object just gets a rho message, computes GHartExt, and sends
54  * vks.
55  */
56 //////////////////////////////////////////////////////////////////////////////
57 CP_Rho_RHartExt::CP_Rho_RHartExt(int _ngrida, int _ngridb, int _ngridc,
58  int _ees_eext_on, int _natmTyp,
59  UberCollection _instance)
60  : thisInstance(_instance)
61 {
62 //////////////////////////////////////////////////////////////////////////////
63 /// Set some pareameters
64 
65  CPcharmParaInfo *sim = CPcharmParaInfo::get();
66  nplane_rho_x = sim->nplane_rho_x;
67  listSubFlag = sim->listSubFlag;
68  rhoRsubplanes = config.rhoRsubplanes;
69 
70  CkAssert(nplane_rho_x >= rhoRsubplanes); // safety : should already be checked.
71 
72  ngrida = _ngrida;
73  ngridb = _ngridb;
74  ngridc = _ngridc;
75  ees_eext_on = _ees_eext_on;
76  natmTypTot = _natmTyp;
77 
78  nchareHartAtmT = config.nchareHartAtmT;
79  int div = (natmTypTot/nchareHartAtmT);
80  int rem = (natmTypTot % nchareHartAtmT);
81  int max = (thisIndex.z < rem ? thisIndex.z : rem);
82  natmTyp = (thisIndex.z<rem ? div+1 : div);
83  atmTypoff = div*thisIndex.z + max;
84 
85  iplane_ind = thisIndex.y*ngridc + thisIndex.x;
86 
87 //////////////////////////////////////////////////////////////////////////////
88 /// Compute messages sizes and zero message counters
89 
90  int nchareG = sim->nchareRhoGEext;
91  if(rhoRsubplanes>1){
92  recvCountFromGHartExt = 0;
93  for(int i=0;i<nchareG;i++){
94  if(sim->nline_send_eext_y[i][thisIndex.y]>0){recvCountFromGHartExt++;}
95  }//endfor
96  }else{
97  recvCountFromGHartExt=nchareG;
98  }//endif
99 
100  countIntRtoG = 0;
101  countIntGtoR[0] = 0;
102  countIntGtoR[1] = 0;
103  countFFT[0] = 0;
104  countFFT[1] = 0;
105  nAtmTypRecv = 0;
106  registrationFlag = 0;
107  launchFlag = 0;
108  iteration = 0;
109  iterAtmTyp = 0;
110  countDebug = 0;
111 
112 //////////////////////////////////////////////////////////////////////////////
113 /// Parallelization
114 
115  //------------------------------------------------------
116  // rho(x,y,z) by (y,z)
117  div = (ngridb/rhoRsubplanes);
118  rem = (ngridb % rhoRsubplanes);
119  max = (thisIndex.y < rem ? thisIndex.y : rem);
120  myNgridb = (thisIndex.y<rem ? div+1 : div); // number of y values/lines of x
121  myBoff = div*thisIndex.y + max; // offset into y
122  nptsB = ngrida*myNgridb; // size of plane without extra room
123  nptsExpndB = (ngrida+2)*myNgridb; // extra memory for RealToComplex FFT
124 
125  //------------------------------------------------------
126  // Parallelization after transpose : rho(gx,y,z) : parallelize gx and z
127  if(rhoRsubplanes>1){
128  myNplane_rho = sim->numSubGx[thisIndex.y];
129  }else{
130  myNplane_rho = nplane_rho_x;
131  }//endif
132  nptsA = 2*myNplane_rho*ngridb; // memory size for fft in doubles
133  nptsExpndA = 2*myNplane_rho*ngridb; // memory size for fft in doubles
134 
135 
136 //---------------------------------------------------------------------------
137  }//end routine
138 //////////////////////////////////////////////////////////////////////////////
139 
140 //////////////////////////////////////////////////////////////////////////////
141 //////////////////////////////////////////////////////////////////////////////
142 //////////////////////////////////////////////////////////////////////////////
143 void CP_Rho_RHartExt::init(){
144 //////////////////////////////////////////////////////////////////////////////
145 /// Malloc data sets : Register in the cache
146 
147  csize = 0;
148  csizeInt = 0;
149  if(ees_eext_on==1){
150  csize = (ngrida/2+1)*myNgridb; // complex variable size
151 
152  eesCache *eesData = UeesCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
153  eesData->registerCacheRHart(thisIndex.x);
154 
155  atmSFC = (complex*) fftw_malloc(csize*sizeof(complex));
156  atmSFR = reinterpret_cast<double*> (atmSFC);
157  atmForcC = atmSFC;
158  atmForcR = atmSFR;
159 
160  atmEwdSFC = (complex*) fftw_malloc(csize*sizeof(complex));
161  atmEwdSFR = reinterpret_cast<double*> (atmEwdSFC);
162  atmEwdForcC = atmEwdSFC;
163  atmEwdForcR = atmEwdSFR;
164 
165  if(rhoRsubplanes>1){
166  csizeInt = nptsA/2;
167  atmSFCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
168  atmSFRint = reinterpret_cast<double*> (atmSFC);
169  atmForcCint = atmSFCint;
170  atmForcRint = atmSFRint;
171 
172  atmEwdSFCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
173  atmEwdSFRint = reinterpret_cast<double*> (atmEwdSFC);
174  atmEwdForcCint = atmEwdSFCint;
175  atmEwdForcRint = atmEwdSFRint;
176  }//endif
177 
178  int i=1;
179  // CkCallback cb(CkIndex_CP_Rho_RHartExt::registrationDone(NULL),UrhoRHartExtProxy[thisInstance.proxyOffset]);
180  CkCallback cb(CkIndex_CP_Rho_RHartExt::registrationDone(NULL),thisProxy);
181  contribute(sizeof(int),&i,CkReduction::sum_int,cb);
182  }//endif
183 
184 //////////////////////////////////////////////////////////////////////////////
185 /// Set up proxies
186 
187  rhoGHartProxy_com = UrhoGHartExtProxy[thisInstance.proxyOffset];;
188 #ifdef USE_COMLIB
189  if (config.useRHartInsGHart){
190  ComlibAssociateProxy(commRHartGHartIns,rhoGHartProxy_com);
191  }//endif
192 #endif
193 
194 //---------------------------------------------------------------------------
195  }//end routine
196 //////////////////////////////////////////////////////////////////////////////
197 
198 
199 //////////////////////////////////////////////////////////////////////////////
200 //////////////////////////////////////////////////////////////////////////////
201 //////////////////////////////////////////////////////////////////////////////
202 CP_Rho_RHartExt::~CP_Rho_RHartExt(){
203 }
204 //////////////////////////////////////////////////////////////////////////////
205 
206 
207 //////////////////////////////////////////////////////////////////////////////
208 //////////////////////////////////////////////////////////////////////////////
209 //////////////////////////////////////////////////////////////////////////////
210 void CP_Rho_RHartExt::pup(PUP::er &p){
211 //////////////////////////////////////////////////////////////////////////////
212 
213  ArrayElement3D::pup(p);
214  p|listSubFlag;
215  p|countDebug;
216  p|nplane_rho_x;
217  p|rhoRsubplanes;
218  p|registrationFlag;
219  p|launchFlag;
220  p|ngrida;
221  p|ngridb;
222  p|ngridc;
223  p|iplane_ind;
224  p|ees_eext_on;
225  p|natmTypTot;
226  p|atmTypoff;
227  p|natmTyp;
228  p|nchareHartAtmT;
229  p|nAtmTypRecv;
230  p|countIntRtoG;
231  p(countIntGtoR,2);
232  p(countFFT,2);
233  p|iteration;
234  p|iterAtmTyp;
235  p|csize;
236  p|csizeInt;
237  p|myNgridb;
238  p|myBoff;
239  p|nptsB;
240  p|nptsExpndB;
241  p|myNplane_rho;
242  p|myAoff;
243  p|nptsA;
244  p|nptsExpndA;
245 
246  if (p.isUnpacking() && ees_eext_on==1) {
247  atmSFC = (complex*) fftw_malloc(csize*sizeof(complex));
248  atmSFR = reinterpret_cast<double*> (atmSFC);
249  atmForcC = atmSFC;
250  atmForcR = atmSFR;
251  atmEwdSFC = (complex*) fftw_malloc(csize*sizeof(complex));
252  atmEwdSFR = reinterpret_cast<double*> (atmEwdSFC);
253  atmEwdForcC = atmEwdSFC;
254  atmEwdForcR = atmEwdSFR;
255  if(rhoRsubplanes>1){
256  atmSFCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
257  atmSFRint = reinterpret_cast<double*> (atmSFCint);
258  atmForcCint = atmSFC;
259  atmForcRint = atmSFR;
260  atmEwdSFCint = (complex*) fftw_malloc(csizeInt*sizeof(complex));
261  atmEwdSFRint = reinterpret_cast<double*> (atmEwdSFCint);
262  atmEwdForcCint = atmEwdSFCint;
263  atmEwdForcRint = atmEwdSFRint;
264  }//endif
265  }//endif
266 
267  if(ees_eext_on==1) {
268  p((char*)atmSFC,csize*sizeof(complex));
269  p((char*)atmEwdSFC,csize*sizeof(complex));
270  if(rhoRsubplanes>1){
271  p((char*)atmSFCint,csizeInt*sizeof(complex));
272  p((char*)atmEwdSFCint,csizeInt*sizeof(complex));
273  }//endif
274  }//endif
275 
276 //////////////////////////////////////////////////////////////////////////////
277  }//end routine
278 //////////////////////////////////////////////////////////////////////////////
279 
280 
281 //////////////////////////////////////////////////////////////////////////////
282 /// Invoke by Rspace-density : Density has arrived in r-space and will soon arrive
283 /// in g-space. Get moving RhartExt
284 //////////////////////////////////////////////////////////////////////////////
285 //////////////////////////////////////////////////////////////////////////////
286 //////////////////////////////////////////////////////////////////////////////
287 void CP_Rho_RHartExt::startEextIter(){
288 //////////////////////////////////////////////////////////////////////////////
289 /// Check for error
290 
291  CPcharmParaInfo *sim = CPcharmParaInfo::get();
292  int cp_min_opt = sim->cp_min_opt;
293 
294  iteration ++;
295  // the atoms haven't moved yet
296  if(cp_min_opt==0){
297  if(UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch()->iteration != iteration-1){
298  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
299  CkPrintf("Flow of Control Error in GHartExtVks : atoms slow %d %d\n",
300  UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch()->iteration,iteration);
301  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
302  CkExit();
303  }//endif
304  }//endif
305 
306 //////////////////////////////////////////////////////////////////////////////
307 /// This is a new time step : Increment time step counter, zero atm typ counter
308 /// Launch if we are ready
309 
310 #ifdef _CP_RHART_VERBOSE_
311  CkPrintf("HI, I am rhart chare %d %d in start : %d on %d\n",
312  thisIndex.x,thisIndex.y,registrationFlag,CkMyPe());
313 #endif
314 
315  if(iterAtmTyp==natmTyp){CkPrintf("%d signing off\n",thisIndex.x);CkExit();}
316 
317  iterAtmTyp = 0;
318  nAtmTypRecv = 0;
319 
320  launchFlag=1;
321  if(registrationFlag==1){computeAtmSF();}
322 
323 //////////////////////////////////////////////////////////////////////////////
324  }//end routine
325 //////////////////////////////////////////////////////////////////////////////
326 
327 
328 ///////////////////////////////////////////////////////////////////////////=
329 /// Make sure everyone is registered on the 1st time step
330 ///////////////////////////////////////////////////////////////////////////=
331 ///////////////////////////////////////////////////////////////////////////c
332 ///////////////////////////////////////////////////////////////////////////=
333  void CP_Rho_RHartExt::registrationDone(CkReductionMsg *msg) {
334 ///////////////////////////////////////////////////////////////////////////=
335 
336  int sum = ((int *)msg->getData())[0];
337  delete msg;
338 
339 #ifdef _CP_RHART_VERBOSE_
340  CkPrintf("HI, I am rhart chare %d %d in reg : %d on %d\n",
341  thisIndex.x,thisIndex.y,sum,CkMyPe());
342 #endif
343 
344  // use barrier as a debug tool
345  if(iterAtmTyp==natmTyp){
346  CkPrintf("HI, I am rhart chare %d exiting at %d\n",thisIndex.x,iterAtmTyp+1);
347  CkExit();
348  }//endif
349 
350  registrationFlag=1;
351  if(launchFlag==1){computeAtmSF();}
352 
353  }//end routine
354 ///////////////////////////////////////////////////////////////////////////=
355 
356 
357 
358 //////////////////////////////////////////////////////////////////////////////
359 /// Start the real space part of the EES interpolation for atmSF(iatmTyp)
360 //////////////////////////////////////////////////////////////////////////////
361 //////////////////////////////////////////////////////////////////////////////
362 //////////////////////////////////////////////////////////////////////////////
363 void CP_Rho_RHartExt::computeAtmSF(){
364 //////////////////////////////////////////////////////////////////////////////
365 
366  iterAtmTyp++;
367  if(iterAtmTyp>natmTyp){
368  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
369  CkPrintf("Too many iterations RHartExtVks\n");
370  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
371  CkExit();
372  }//endif
373 
374 #ifdef _CP_RHART_VERBOSE_
375  CkPrintf("HI, I am rhart chare %d %d %d in computeatmsf at iter %d on %d\n",
376  thisIndex.x,thisIndex.y,thisIndex.z, iterAtmTyp,CkMyPe());
377 #endif
378 
379 //////////////////////////////////////////////////////////////////////////////
380 /// Look up some constants and fill the r-space grid
381 
382  int myPlane = thisIndex.x;
383  int mySubplane = thisIndex.y;
384  int itime = iteration;
385  int iterAtmTypFull = iterAtmTyp+atmTypoff;
386 
387  eesCache *eesData= UeesCacheProxy[thisInstance.proxyOffset].ckLocalBranch ();
388  if(iterAtmTyp==1){eesData->queryCacheRHart(myPlane,itime,iterAtmTyp);}
389 
390  int *plane_index = eesData->RhoRHartData[myPlane]->plane_index;
391  int **nSub = eesData->RhoRHartData[myPlane]->nSub;
392  int ***igrid = eesData->RhoRHartData[myPlane]->igrid;
393  double ***mn = eesData->RhoRHartData[myPlane]->mn;
394 
395  AtomsCache *ag = UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch(); // find me the local copy
396  int natm = ag->natm;
397 
398  CPLOCAL::eesPackGridRchare(natm,iterAtmTypFull,atmSFR,myPlane,
399  mySubplane,igrid,mn,
400  plane_index,nSub,myNgridb);
401 
402 //////////////////////////////////////////////////////////////////////////////
403 /// FFT the result to G-space
404 
405  fftAtmSfRtoG();
406 
407 //////////////////////////////////////////////////////////////////////////////
408  }//end routine
409 //////////////////////////////////////////////////////////////////////////////
410 
411 
412 //////////////////////////////////////////////////////////////////////////////
413 //////////////////////////////////////////////////////////////////////////////
414 //////////////////////////////////////////////////////////////////////////////
415 /// FFT of SF(x,y,z,iatmTyp) -> SF(gx,gy,z,iatmTyp)
416 //////////////////////////////////////////////////////////////////////////////
417 void CP_Rho_RHartExt::fftAtmSfRtoG(){
418 //////////////////////////////////////////////////////////////////////////////
419 #ifdef _CP_RHART_VERBOSE_
420  CkPrintf("HI, I am rhart chare %d %d %d in fftsback at %d on %d\n",
421  thisIndex.x,thisIndex.y,thisIndex.z, iterAtmTyp,CkMyPe());
422 #endif
423 
424 ///////////////////////////////////////////////////////////////////////////=
425 /// Do the FFT
426 #if CMK_TRACE_ENABLED
427  double StartTime=CmiWallTimer();
428 #endif
429 
430  FFTcache *fftcache = UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
431 
432  if(rhoRsubplanes>1){
433  fftcache->doEextFFTRxToGx_Rchare(atmSFC,atmSFR,nplane_rho_x,ngrida,myNgridb,iplane_ind);
434  }else{
435  fftcache->doEextFFTRtoG_Rchare(atmSFC,atmSFR,nplane_rho_x,ngrida,ngridb,iplane_ind);
436  }//endif
437 
438 #if CMK_TRACE_ENABLED
439  traceUserBracketEvent(doEextFFTRtoG_, StartTime, CmiWallTimer());
440 #endif
441 
442 ///////////////////////////////////////////////////////////////////////////=
443 /// Do the communication if you are not debugging
444 
445  if(rhoRsubplanes>1){
446  sendAtmSfRyToGy(); // double transpose method (yz ---> gx,z)
447  }else{
448 #ifdef DEBUG_INT_TRANS_FWD
449  char name[100];
450  sprintf(name,"partFFTGxGyZ%d.out.%d.%d",rhoRsubplanes,thisIndex.x,thisIndex.y);
451  FILE *fp = fopen(name,"w");
452  for(int ix =0;ix<nplane_rho_x;ix++){
453  for(int iy =0;iy<ngridb;iy++){
454  int i = iy*(ngrida/2+1) + ix;
455  fprintf(fp,"%d %d : %g %g\n",iy,ix,atmSFC[i].re,atmSFC[i].im);
456  }//endfor
457  }//endof
458  fclose(fp);
459  UrhoRHartExtProxy[thisInstance.proxyOffset](0,0,0).exitForDebugging();
460 #else
461 #ifdef _GLENN_DBG_KPT_
462  char name[100];
463  sprintf(name,"sfAtmTypGxGyZ_inR.p%d.t%d.out",thisIndex.x,iterAtmTyp);
464  FILE *fp = fopen(name,"w");
465  for(int ix =0;ix<nplane_rho_x;ix++){
466  for(int iy =0;iy<ngridb;iy++){
467  int i = iy*(ngrida/2+1) + ix;
468  fprintf(fp,"%d %d : %g %g\n",iy,ix,atmSFC[i].re,atmSFC[i].im);
469  }//endfor
470  }//endof
471  fclose(fp);
472 #endif
473  sendAtmSfRhoGHart(); // single transpose method (z ---> gx,gy)
474 #endif
475  }//endif
476 
477 //////////////////////////////////////////////////////////////////////////////
478  }// end routine
479 //////////////////////////////////////////////////////////////////////////////
480 
481 
482 
483 //////////////////////////////////////////////////////////////////////////////
484 //////////////////////////////////////////////////////////////////////////////
485 //////////////////////////////////////////////////////////////////////////////
486 ///
487 /// Double Transpose Fwd Send : A(gx,y,z) on the way to A(gx,gy,z)
488 /// Send so that (y,z) parallelism is
489 /// switched to (gx,z)
490 ///
491 //////////////////////////////////////////////////////////////////////////////
492 void CP_Rho_RHartExt::sendAtmSfRyToGy(){
493 //////////////////////////////////////////////////////////////////////////////
494 /// Launch the communication
495 
496  CkAssert(rhoRsubplanes>1);
497  CPcharmParaInfo *sim = CPcharmParaInfo::get();
498  int **listSubGx = sim->listSubGx;
499  int *numSubGx = sim->numSubGx;
500 
501  //-----------------------------------------------------------------------------
502  // Commlib launch :
503 
504 #ifdef _ERIC_SETS_UP_COMMLIB_
505  if(config.useRInsRhoRP) commRealInstanceRx.beginIteration();
506 #endif
507 
508  //-----------------------------------------------------------------------------
509  // Send the data : I have myNgridB values of y (gx,y) y=1...myNgridB and all gx
510  // Send all the `y' I have for the gx range desired after transpose
511 
512  int stride = ngrida/2+1;
513  int ix = thisIndex.x;
514  for(int ic = 0; ic < rhoRsubplanes; ic ++) { // chare arrays to which we will send
515  int num = numSubGx[ic]; // number of gx values chare ic wants
516  int size = num*myNgridb; // num*(all the y's i have)
517 
518  int sendFFTDataSize = size;
519  RhoGHartMsg *msg = new (sendFFTDataSize, 8 * sizeof(int)) RhoGHartMsg;
520  msg->size = size;
521  msg->offset = myBoff; // where the myNgridB y-lines start.
522  msg->num = myNgridb; // number of y-lines I have.
523  msg->iter = iterAtmTyp;
524  complex *data = msg->data; // data
525 
526  if(config.prioFFTMsg){
527  CkSetQueueing(msg, CK_QUEUEING_IFIFO);
528  *(int*)CkPriorityPtr(msg) = config.rhogpriority+thisIndex.y;
529  }//endif
530 
531  if(listSubFlag==1){
532  for(int i=0,koff=0;i<num;i++,koff+=myNgridb){
533  for(int k=koff,ii=listSubGx[ic][i];k<myNgridb+koff;k++,ii+=stride){
534  data[k] = atmSFC[ii]; // all y's of this gx
535  }//endfor
536  }//endfor
537  }else{
538  int nst=listSubGx[ic][0];
539  for(int i=0,ist=nst,koff=0;i<num;i++,koff+=myNgridb,ist++){
540  for(int k=koff,ii=ist;k<myNgridb+koff;k++,ii+=stride){
541  data[k] = atmSFC[ii]; // all y's of this gx
542  }//endfor
543  }//endfor
544  }//endif
545 
546  UrhoRHartExtProxy[thisInstance.proxyOffset](ix,ic,thisIndex.z).recvAtmSfRyToGy(msg);
547 
548 #ifdef _ERIC_SETS_UP_COMMLIB_
549  UrhoRHartExtProxy[thisInstance.proxyOffset](ix,ic,thisIndex.z).recvAtmSfRyToGy(msg);
550 #endif
551 
552  }//end for : chare sending
553 
554  //-----------------------------------------------------------------------------
555  // Commlib stop
556 
557 #ifdef _ERIC_SETS_UP_COMMLIB_
558  if(config.useRInsRhoRP) commRealInstanceRx.endIteration();
559 #endif
560 
561 //---------------------------------------------------------------------------
562  }//end routine
563 //////////////////////////////////////////////////////////////////////////////
564 
565 
566 
567 //////////////////////////////////////////////////////////////////////////////
568 //////////////////////////////////////////////////////////////////////////////
569 //////////////////////////////////////////////////////////////////////////////
570 ///
571 /// Double Transpose Fwd Recv : A(gx,y,z) on the way to A(gx,gy,z)
572 /// Recv so that (y,z) parallel switched to (gx,z)
573 ///
574 /// Invoked natm_typ times per algorithm step :
575 ///
576 //////////////////////////////////////////////////////////////////////////////
577 void CP_Rho_RHartExt::recvAtmSfRyToGy(RhoGHartMsg *msg){
578 //////////////////////////////////////////////////////////////////////////////
579 
580  int size = msg->size; // msg size
581  int iter = msg->iter; //
582  int num = msg->num;
583  int offset = msg->offset;
584  complex *msgData = msg->data;
585 
586  CkAssert(size==myNplane_rho*num);
587  CkAssert(rhoRsubplanes>1);
588 
589 //////////////////////////////////////////////////////////////////////////////
590 
591  for(int js=0,j=offset;js<size;js+=num,j+=ngridb){
592  for(int is=js,i=j;is<num+js;is++,i++){
593  atmSFCint[i] = msgData[is];
594  }//endfor
595  }//endfor
596 
597  delete msg;
598 
599 //////////////////////////////////////////////////////////////////////////////
600 /// Do the Y fft, invoke communication if not debugging
601 
602  countIntRtoG++;
603  if(countIntRtoG==rhoRsubplanes){
604 
605  countIntRtoG = 0;
606  FFTcache *fftcache = UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
607 #if CMK_TRACE_ENABLED
608  double StartTime=CmiWallTimer();
609 #endif
610  fftcache->doEextFFTRyToGy_Rchare(atmSFCint,atmSFRint,myNplane_rho,
611  ngrida,ngridb,iplane_ind);
612 #if CMK_TRACE_ENABLED
613  traceUserBracketEvent(doEextFFTRytoGy_, StartTime, CmiWallTimer());
614 #endif
615 
616 #ifndef DEBUG_INT_TRANS_FWD
617  sendAtmSfRhoGHart();
618 #endif
619  }//endif
620 
621 //////////////////////////////////////////////////////////////////////////////
622 /// Debugging
623 
624 #ifdef DEBUG_INT_TRANS_FWD
625  if(countIntRtoG==rhoRsubplanes){
626  CPcharmParaInfo *sim = CPcharmParaInfo::get();
627  int **listSubGx = sim->listSubGx;
628  int ic = thisIndex.y;
629  char name[100];
630  sprintf(name,"partFFTGxGyZ%d.out.%d.%d",rhoRsubplanes,thisIndex.x,thisIndex.y);
631  FILE *fp = fopen(name,"w");
632  for(int ix =0;ix<myNplane_rho;ix++){
633  for(int iy =0;iy<ngridb;iy++){
634  int i = ix*ngridb + iy;
635  fprintf(fp,"%d %d : %g %g\n",iy,listSubGx[ic][ix],atmSFCint[i].re,atmSFCint[i].im);
636  }//endfor
637  }//endof
638  fclose(fp);
639  UrhoRHartExtProxy[thisInstance.proxyOffset](0,0,0).exitForDebugging();
640  }//endif
641 #endif
642 
643 //////////////////////////////////////////////////////////////////////////////
644  }//end routine
645 //////////////////////////////////////////////////////////////////////////////
646 
647 
648 
649 //////////////////////////////////////////////////////////////////////////////
650 //////////////////////////////////////////////////////////////////////////////
651 //////////////////////////////////////////////////////////////////////////////
652 /// Send SF(gx,gy,z,iatmTYP) to g-space whence the FFT will be completed
653 //////////////////////////////////////////////////////////////////////////////
654 void CP_Rho_RHartExt::sendAtmSfRhoGHart(){
655 //////////////////////////////////////////////////////////////////////////////
656 
657 #ifdef _CP_RHART_VERBOSE_
658  CkPrintf("HI, I am rhart chare %d %d %d in sendsf at %d on %d\n",
659  thisIndex.x,thisIndex.y,thisIndex.z,iterAtmTyp,CkMyPe());
660 #endif
661 
662  CPcharmParaInfo *sim = CPcharmParaInfo::get();
663  int nchareG = sim->nchareRhoGEext;
664 
665  int **tranpack = sim->index_tran_upack_eext;
666  int *nlines_per_chareG = sim->nlines_per_chareRhoGEext;
667 
668  int ***tranupack_Y = sim->index_tran_upack_eext_y;
669  int **nlines_per_chareG_Y = sim->nline_send_eext_y;
670 
671 ///////////////////////////////////////////////////////////////////==
672 /// Perform the transpose and then the blast off the final 1D-FFT
673 
674  //----------------------------------------------------------------
675  // start commlib
676 #ifdef USE_COMLIB
677  if(rhoRsubplanes==1){
678 #ifdef OLD_COMMLIB
679  if(config.useRHartInsGHart){commRHartGHartIns.beginIteration();}
680 #else
681  if(config.useRHartInsGHart){ComlibBegin(rhoGHartProxy_com,0);}
682 #endif
683  }//endif
684 #endif
685 
686  //----------------------------------------------------------------
687  // do the send
688 
689  int iy = thisIndex.y;
690  for (int ic=0;ic<nchareG;ic++) { // chare arrays to which we will send
691 
692  //---------------------
693  // malloc the message
694  int sendFFTDataSize = nlines_per_chareG[ic];
695  if(rhoRsubplanes!=1){sendFFTDataSize = nlines_per_chareG_Y[ic][iy];}
696  if(sendFFTDataSize>0)
697  {
698  RhoGHartMsg *msg = new (sendFFTDataSize, 8 * sizeof(int)) RhoGHartMsg;
699 
700  //----------------------
701  // pack the message
702  if(config.prioEextFFTMsg){
703  CkSetQueueing(msg, CK_QUEUEING_IFIFO);
704  *(int*)CkPriorityPtr(msg) = config.rhogHartpriority+thisIndex.x*10;
705  }//endif
706  msg->iter = iterAtmTyp;
707  msg->size = sendFFTDataSize;
708  msg->offset = thisIndex.x; // c-plane-index
709  msg->offsetGx = thisIndex.y; // gx parallelization index
710  complex *data = msg->data;
711 
712  if(rhoRsubplanes==1){
713  for(int i=0;i<sendFFTDataSize;i++){data[i] = atmSFC[tranpack[ic][i]];}
714  }else{
715  for(int i=0;i<sendFFTDataSize;i++){
716  data[i] = atmSFCint[tranupack_Y[ic][iy][i]];
717  }//endfor
718  }//endif
719 
720  //-----------------
721  // Send the message
722  if(rhoRsubplanes==1){
723  rhoGHartProxy_com(ic,thisIndex.z).recvAtmSFFromRhoRHart(msg); // send the message
724  }else{
725  UrhoGHartExtProxy[thisInstance.proxyOffset](ic,thisIndex.z).recvAtmSFFromRhoRHart(msg); // send the message
726  }//endif
727  }
728  }//end for : chare sending
729 
730  //----------------------------------------------------------------
731  // stop commlib
732 #ifdef USE_COMLIB
733  if(rhoRsubplanes==1){
734 #ifdef OLD_COMMLIB
735  if(config.useRHartInsGHart){commRHartGHartIns.endIteration();}
736 #else
737  if(config.useRHartInsGHart){ComlibBegin(rhoGHartProxy_com,1);}
738 #endif
739  }//endif
740 #endif
741 
742 //////////////////////////////////////////////////////////////////////////////
743  }// end routine
744 //////////////////////////////////////////////////////////////////////////////
745 
746 
747 //////////////////////////////////////////////////////////////////////////////
748 /// Hartree sends back atom forces from e-atm interation
749 /// Depending on the flag, it is Ewald or e-atm interation
750 //////////////////////////////////////////////////////////////////////////////
751 //////////////////////////////////////////////////////////////////////////////
752 //////////////////////////////////////////////////////////////////////////////
753 void CP_Rho_RHartExt::recvAtmForcFromRhoGHart(RhoRHartMsg *msg){
754 //////////////////////////////////////////////////////////////////////////////
755 
756  CPcharmParaInfo *sim = CPcharmParaInfo::get();
757  int nchareG = sim->nchareRhoGEext;
758 
759  int **tranUnpack = sim->index_tran_upack_eext;
760  int *nlines_per_chareG = sim->nlines_per_chareRhoGEext;
761 
762  int ***tranupack_Y = sim->index_tran_upack_eext_y;
763  int **nlines_per_chareG_Y = sim->nline_send_eext_y;
764  int iy = thisIndex.y;
765 
766  int size = msg->size;
767  int Index = msg->senderIndex;
768  int iopt = msg->iopt;
769  int iter = msg->iter;
770  complex *partiallyFFTd = msg->data;
771 
772  int mySize;
773  int csizenow;
774  if(rhoRsubplanes==1){
775  mySize = nlines_per_chareG[Index];
776  csizenow = csize;
777  }else{
778  mySize = nlines_per_chareG_Y[Index][iy];
779  csizenow = csizeInt;
780  }//endif
781 
782 #ifdef _NAN_CHECK_
783  for(int i=0;i<msg->size ;i++){
784  CkAssert(isnan(msg->data[i].re)==0);
785  CkAssert(isnan(msg->data[i].im)==0);
786  }
787 #endif
788 
789 //////////////////////////////////////////////////////////////////////////////
790 /// Perform some error checking
791 
792  countFFT[iopt]++;
793  if (countFFT[iopt] > recvCountFromGHartExt) {
794  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
795  CkPrintf("Mismatch in allowed rho_gspace chare arrays : %d %d %d %d\n",
796  countFFT[iopt],nchareG,thisIndex.x,thisIndex.y);
797  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
798  CkExit();
799  }//endif
800 
801  if(size!=mySize){
802  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
803  CkPrintf("Dude.1, %d != %d for rho rhart chare %d %d\n",size,mySize,
804  thisIndex.x,Index);
805  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
806  CkExit();
807  }//endif
808 
809  if(iter!=iterAtmTyp){
810  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
811  CkPrintf("Dude.2, %d != %d for rho rhart chare %d %d %d\n",iter,iterAtmTyp,
812  thisIndex.x,Index,iopt);
813  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
814  CkExit();
815  }//endif
816 
817  if((iopt==1) && (iterAtmTyp!=natmTyp) ){
818  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
819  CkPrintf("Dude.3, %d != %d for rho rhart chare %d %d : %d %d\n",iter,natmTyp,
820  thisIndex.x,iy,Index,iopt);
821  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
822  CkExit();
823  }//endif
824 
825 //////////////////////////////////////////////////////////////////////////////
826 /// unpack the data and delete the message
827 
828  complex *data;
829  if(rhoRsubplanes==1){
830  switch(iopt){
831  case 0 : data = atmSFC; break;
832  case 1 : data = atmEwdSFC; break;
833  default: CkAbort("impossible iopt"); break;
834  }//end switch
835  }else{
836  switch(iopt){
837  case 0 : data = atmSFCint; break;
838  case 1 : data = atmEwdSFCint; break;
839  default: CkAbort("impossible iopt"); break;
840  }//end switch
841  }//endif
842 
843  if(countFFT[iopt]==1){memset(data,0,sizeof(complex)*csizenow);}
844 
845  if(rhoRsubplanes==1){
846  for(int i=0;i<size;i++){data[tranUnpack[Index][i]] = partiallyFFTd[i];}
847  }else{
848  for(int i=0;i<size;i++){
849  data[tranupack_Y[Index][iy][i]] = partiallyFFTd[i];
850  }//endfor
851  }//endif
852 
853  delete msg;
854 
855 //////////////////////////////////////////////////////////////////////////////
856 /// When you have all the data : finish the FFT back to real space
857 
858  if (countFFT[iopt] == recvCountFromGHartExt){
859 #ifdef _CP_RHART_VERBOSE_
860  CkPrintf("HI, I am rhart chare %d %d %d in recvsf with opt %d at %d on %d\n",
861  thisIndex.x,thisIndex.y, thisIndex.z,iopt,iterAtmTyp,CkMyPe());
862 #endif
863  countFFT[iopt]=0;
864  if(rhoRsubplanes==1){nAtmTypRecv++;}
865  fftAtmForcGtoR(iopt);
866  }//endif
867 
868 //////////////////////////////////////////////////////////////////////////////
869  }//end routine
870 //////////////////////////////////////////////////////////////////////////////
871 
872 
873 //////////////////////////////////////////////////////////////////////////////
874 /// Complete the FFT : sf(gx,gy,z) ->sf(x,y,z)
875 //////////////////////////////////////////////////////////////////////////////
876 //////////////////////////////////////////////////////////////////////////////
877 //////////////////////////////////////////////////////////////////////////////
878 void CP_Rho_RHartExt::fftAtmForcGtoR(int flagEwd){
879 //////////////////////////////////////////////////////////////////////////////
880 
881 #ifdef _CP_RHART_VERBOSE_
882  CkPrintf("HI, I am rhart chare %d %d %d in FFTAtmForc w %d at %d on %d\n",
883  thisIndex.x,thisIndex.y,thisIndex.z,flagEwd,iterAtmTyp,CkMyPe());
884 #endif
885 
886  double *dataR;
887  complex *dataC;
888  if(rhoRsubplanes==1){
889  switch(flagEwd){
890  case 0: dataR=atmSFR; dataC=atmSFC; break;
891  case 1: dataR=atmEwdSFR; dataC=atmEwdSFC; break;
892  default: CkAbort("impossible iopt"); break;
893  }//endif
894  }else{
895  switch(flagEwd){
896  case 0: dataR=atmSFRint; dataC=atmSFCint; break;
897  case 1: dataR=atmEwdSFRint; dataC=atmEwdSFCint; break;
898  default: CkAbort("impossible iopt"); break;
899  }//endif
900  }//endif
901 
902 //////////////////////////////////////////////////////////////////////////////
903 
904 #if CMK_TRACE_ENABLED
905  double StartTime=CmiWallTimer();
906 #endif
907 
908  FFTcache *fftcache = UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
909  if(rhoRsubplanes==1){
910  fftcache->doEextFFTGtoR_Rchare(dataC,dataR,nplane_rho_x,ngrida,ngridb,iplane_ind);
911  computeAtmForc(flagEwd);
912  }else{
913  fftcache->doEextFFTGyToRy_Rchare(dataC,dataR,myNplane_rho,ngrida,ngridb,iplane_ind);
914  sendAtmForcGxToRx(flagEwd);
915  }//endif
916 
917 #if CMK_TRACE_ENABLED
918  traceUserBracketEvent(doEextFFTGtoR_, StartTime, CmiWallTimer());
919 #endif
920 
921 //////////////////////////////////////////////////////////////////////////////
922  }//end routine
923 //////////////////////////////////////////////////////////////////////////////
924 
925 
926 
927 //////////////////////////////////////////////////////////////////////////////
928 //////////////////////////////////////////////////////////////////////////////
929 //////////////////////////////////////////////////////////////////////////////
930 void CP_Rho_RHartExt::sendAtmForcGxToRx(int iopt){
931 //////////////////////////////////////////////////////////////////////////////
932  CkAssert(rhoRsubplanes>1);
933 
934  complex *FFTresult;
935  switch(iopt){
936  case 0 : FFTresult = atmSFCint; break;
937  case 1 : FFTresult = atmEwdSFCint; break;
938  default: CkAbort("impossible iopt"); break;
939  }//end switch
940 
941 //////////////////////////////////////////////////////////////////////////////
942 /// Launch the communication
943 
944  //-----------------------------------------------------------------------------
945  // Commlib launch :
946 
947 #ifdef _ERIC_SETS_UP_COMMLIB_
948  switch(iopt){
949  case 0 : if(config.useRInsRhoRP) commRealInstanceRx.beginIteration(); break;
950  case 1 : if(config.useRInsIGXRhoRP) commRealIGXInstanceRx.beginIteration(); break;
951  default: CkAbort("impossible iopt");break;
952  }//end switch
953 #endif
954 
955  //-----------------------------------------------------------------------------
956  // Send the data : I have myNgridB values of y (gx,y) y=1...myNgridB and all gx
957  // Send all the `y' I have for the gx range desired after transpose
958 
959  int ix = thisIndex.x;
960  for(int ic = 0; ic < rhoRsubplanes; ic ++) { // chare arrays to which we will send
961  int div = (ngridb/rhoRsubplanes); //parallelize y
962  int rem = (ngridb % rhoRsubplanes);
963  int add = (ic < rem ? 1 : 0);
964  int max = (ic < rem ? ic : rem);
965  int ist = div*ic + max; // start of y desired by chare ic
966  int iend = ist + div + add; // end of y desired by chare ic
967  int size = (iend-ist)*myNplane_rho; // data size : send all the gx I have
968 
969  int sendFFTDataSize = size;
970  RhoGHartMsg *msg = new (sendFFTDataSize, 8 * sizeof(int)) RhoGHartMsg;
971  msg->size = size;
972  msg->iopt = iopt;
973  msg->iter = iterAtmTyp;
974  msg->offset = thisIndex.y; // my chare index
975  msg->num = myNplane_rho; // number of gx-lines I have.
976  complex *data = msg->data; // data
977 
978  if(config.prioFFTMsg){
979  CkSetQueueing(msg, CK_QUEUEING_IFIFO);
980  *(int*)CkPriorityPtr(msg) = config.rhogpriority+thisIndex.y;
981  }//endif
982 
983  for(int i=ist,koff=0;i<iend;i++,koff+=myNplane_rho){
984  for(int k=koff,ii=i;k<myNplane_rho+koff;k++,ii+=ngridb){
985  data[k] = FFTresult[ii];
986  }//endfor
987  }//endfor
988 
989  //HEY is thisIndex.z right here?
990  UrhoRHartExtProxy[thisInstance.proxyOffset](ix,ic,thisIndex.z).recvAtmForcGxToRx(msg);
991 
992 #ifdef _ERIC_SETS_UP_COMMLIB_
993  switch(iopt){
994  case 0 : UrhoGProxy[thisInstance.proxyOffset]_com(ic,0).acceptRhoGradVksGxToRx(msg);break;
995  case 1 : UrhoGProxy[thisInstance.proxyOffset]IGX_com(ic,0).acceptRhoGradVksGxToRx(msg); break;
996  default: CkAbort("impossible iopt");break;
997  }//end switch
998 #endif
999 
1000  }//end for : chare sending
1001 
1002  //-----------------------------------------------------------------------------
1003  // Commlib stop
1004 
1005 #ifdef _ERIC_SETS_UP_COMMLIB_
1006  switch(iopt){
1007  case 0 : if(config.useRInsRhoRP) commRealInstanceRx.endIteration(); break;
1008  case 1 : if(config.useRInsIGXRhoRP) commRealIGXInstanceRx.endIteration(); break;
1009  default: CkAbort("impossible iopt");break;
1010  }//end switch
1011 #endif
1012 
1013 //---------------------------------------------------------------------------
1014  }//end routine
1015 //////////////////////////////////////////////////////////////////////////////
1016 
1017 
1018 //////////////////////////////////////////////////////////////////////////////
1019 //////////////////////////////////////////////////////////////////////////////
1020 //////////////////////////////////////////////////////////////////////////////
1021 void CP_Rho_RHartExt::recvAtmForcGxToRx(RhoGHartMsg *msg){
1022 //////////////////////////////////////////////////////////////////////////////
1023 
1024  CPcharmParaInfo *sim = CPcharmParaInfo::get();
1025  int **listSubGx = sim->listSubGx;
1026  int *numSubGx = sim->numSubGx;
1027 
1028  int size = msg->size; // msg size
1029  int iopt = msg->iopt;
1030  int iter = msg->iter;
1031  int num = msg->num; // number of lines along `a' sent
1032  int offset = msg->offset; // offset into lines along `a'
1033  complex *msgData = msg->data;
1034 
1035  CkAssert(rhoRsubplanes>1);
1036  CkAssert(size==myNgridb*num);
1037  if(iterAtmTyp!=iter){
1038  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1039  CkPrintf("iter!=iterAtmTyp recvatmforcgxtorx : %d %d %d : %d %d %d\n",
1040  thisIndex.x,thisIndex.y,thisIndex.z,iterAtmTyp,iter,iopt);
1041  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1042  CkExit();
1043  }//endif
1044 
1045  complex *dataC;
1046  double *dataR;
1047  switch(iopt){
1048  case 0: dataR=atmSFR; dataC=atmSFC; break;
1049  case 1: dataR=atmEwdSFR; dataC=atmEwdSFC; break;
1050  default: CkAbort("impossible iopt"); break;
1051  }//end switch
1052 
1053 //////////////////////////////////////////////////////////////////////////////
1054 /// Unpack the message
1055 
1056  countIntGtoR[iopt]++;
1057  if(countIntGtoR[iopt]==1){bzero(dataR,sizeof(double)*nptsExpndB);}
1058 
1059  int stride = ngrida/2+1;
1060  if(listSubFlag==1){
1061  for(int js=0,j=0;js<size;js+=num,j++){
1062  int jj = j*stride;
1063  for(int is=js,i=0;is<(num+js);is++,i++){
1064  dataC[(listSubGx[offset][i]+jj)] = msgData[is];
1065  }//endfor
1066  }//endfor
1067  }else{
1068  int nst = listSubGx[offset][0];
1069  for(int js=0,j=0;js<size;js+=num,j++){
1070  int jj = j*stride+nst;
1071  for(int is=js,i=jj;is<(num+js);is++,i++){
1072  dataC[i] = msgData[is];
1073  }//endfor
1074  }//endfor
1075  }//endif
1076  delete msg;
1077 
1078 //////////////////////////////////////////////////////////////////////////////
1079 /// Do the FFT when you have all the parts
1080 
1081  if(countIntGtoR[iopt]==rhoRsubplanes){
1082  countIntGtoR[iopt]=0;
1083  FFTcache *fftcache = UfftCacheProxy[thisInstance.proxyOffset].ckLocalBranch();
1084 #if CMK_TRACE_ENABLED
1085  double StartTime=CmiWallTimer();
1086 #endif
1087  fftcache->doEextFFTGxToRx_Rchare(dataC,dataR,nplane_rho_x,ngrida,myNgridb,iplane_ind);
1088 #if CMK_TRACE_ENABLED
1089  traceUserBracketEvent(doEextFFTGxtoRx_, StartTime, CmiWallTimer());
1090 #endif
1091 
1092  nAtmTypRecv++;
1093  computeAtmForc(iopt);
1094  }//endif
1095 
1096 //////////////////////////////////////////////////////////////////////////////
1097  }//end routine
1098 //////////////////////////////////////////////////////////////////////////////
1099 
1100 
1101 
1102 //////////////////////////////////////////////////////////////////////////////
1103 /// Get forces on atoms from Ewald or electrons
1104 //////////////////////////////////////////////////////////////////////////////
1105 //////////////////////////////////////////////////////////////////////////////
1106 //////////////////////////////////////////////////////////////////////////////
1107 void CP_Rho_RHartExt::computeAtmForc(int flagEwd){
1108 //////////////////////////////////////////////////////////////////////////////
1109 #ifdef _CP_RHART_VERBOSE_
1110  CkPrintf("HI, I am rhart chare %d %d %d in computeAtmForc.1 w %d at %d %d %d on %d\n",
1111  thisIndex.x,thisIndex.y, thisIndex.z, flagEwd,iterAtmTyp,nAtmTypRecv,natmTyp,CkMyPe());
1112 #endif
1113 //////////////////////////////////////////////////////////////////////////////
1114 /// Check for errors
1115 
1116  // Ewald total SF arrives on the last iteration only .
1117  // In parallel it can arrive before the last atmtyp SF
1118  if( (flagEwd==1) && (iterAtmTyp!=natmTyp)){
1119  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1120  CkPrintf("You can't have flagEwd==1 unless you are on the last step\n");
1121  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1122  CkExit();
1123  }//endif
1124 
1125  // If we are not on the last iteration when ewald can show up
1126  // the iteration # must match the number of SF received
1127  if( (flagEwd==0) && (iterAtmTyp!=nAtmTypRecv) && (iterAtmTyp!=natmTyp)){
1128  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1129  CkPrintf("Atm SF Recv and atm SF calc out of sync %d\n",thisIndex.x);
1130  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1131  CkExit();
1132  }//endif
1133 
1134  // the number of iteration <=natmTyp and you recieve natmtyp SF plus Ewald total
1135  if( iterAtmTyp>natmTyp || nAtmTypRecv>natmTyp+1 ){
1136  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1137  CkPrintf("Too much action for rhoRhart\n");
1138  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
1139  CkExit();
1140  }//endif
1141 
1142 //////////////////////////////////////////////////////////////////////////////
1143 /// set the variables to evaluate atom forces
1144 
1145  int myPlane = thisIndex.x;
1146  int mySubplane = thisIndex.y;
1147  int iterAtmTypFull = iterAtmTyp+atmTypoff;
1148 
1149  double *data;
1150  if(flagEwd==0){data=atmSFR;}else{data=atmEwdSFR;}
1151 
1152  // you have already queried for this step:
1153  eesCache *eesData = UeesCacheProxy[thisInstance.proxyOffset].ckLocalBranch ();
1154  AtomsCache *ag = UatomsCacheProxy[thisInstance.proxyOffset].ckLocalBranch(); // find me the local copy
1155 
1156  CkAssert(eesData->allowedRhoRHartChares[myPlane]!=0);
1157  int *plane_index = eesData->RhoRHartData[myPlane]->plane_index;
1158  int **nSub = eesData->RhoRHartData[myPlane]->nSub;
1159  int ***igrid = eesData->RhoRHartData[myPlane]->igrid;
1160  double ***dmn_x = eesData->RhoRHartData[myPlane]->dmn_x;
1161  double ***dmn_y = eesData->RhoRHartData[myPlane]->dmn_y;
1162  double ***dmn_z = eesData->RhoRHartData[myPlane]->dmn_z;
1163 
1164  FastAtoms *fastAtoms = &(ag->fastAtoms);
1165  int natm = ag->natm;
1166 
1167 #if CMK_TRACE_ENABLED
1168  double StartTime=CmiWallTimer();
1169 #endif
1170 
1171  CPLOCAL::eesAtmForceRchare(natm,fastAtoms,iterAtmTypFull,igrid,
1172  dmn_x,dmn_y,dmn_z, plane_index,nSub,data,
1173  myPlane,mySubplane,flagEwd);
1174 #if CMK_TRACE_ENABLED
1175  traceUserBracketEvent(eesAtmForcR_, StartTime, CmiWallTimer());
1176 #endif
1177 
1178 //////////////////////////////////////////////////////////////////////////////
1179 /// If you aren't done, start another round.
1180 
1181  if(iterAtmTyp<natmTyp){computeAtmSF();}
1182 
1183 //////////////////////////////////////////////////////////////////////////////
1184 /// If your are done, Tell rho real and reset yourself for the next time step
1185 
1186  int nrecvMax=natmTyp;
1187  if(thisIndex.z==0){nrecvMax++;}
1188 
1189 #define _CP_DEBUG_STOP_RHART_OFF
1190 
1191 #ifndef _CP_DEBUG_STOP_RHART_
1192  if(iterAtmTyp==natmTyp && nAtmTypRecv==nrecvMax){
1193 #ifdef _CP_RHART_VERBOSE_
1194  CkPrintf("HI, I am rhart chare %d %d %d in computeAtmForc.2 w %d at %d done on %d\n",
1195  thisIndex.x,thisIndex.y,thisIndex.z,flagEwd,iterAtmTyp,CkMyPe());
1196 #endif
1197  CPcharmParaInfo *sim = CPcharmParaInfo::get();
1198  int index = (thisIndex.x % sim->sizeZ);
1199  UrhoRealProxy[thisInstance.proxyOffset](index,thisIndex.y).RHartReport();
1200  iterAtmTyp = 0;
1201  nAtmTypRecv = 0;
1202  }//endif
1203 #endif
1204 
1205 //////////////////////////////////////////////////////////////////////////////
1206 /// If your are done and debugging, exit
1207 
1208 #ifdef _CP_DEBUG_STOP_RHART_
1209  if(iterAtmTyp==natmTyp && nAtmTypRecv==nrecvMax){
1210  CkPrintf("HI, I am rhart chare %d in computeAtmForc w %d at %d contrib\n",
1211  thisIndex.x,flagEwd,iterAtmTyp);
1212  int i=1;
1213  CkCallback cb(CkIndex_CP_Rho_RHartExt::registrationDone(NULL),UrhoRHartExtProxy[thisInstance.proxyOffset]);
1214  contribute(sizeof(int),&i,CkReduction::sum_int,cb);
1215  }//endif
1216 #endif
1217 
1218 //////////////////////////////////////////////////////////////////////////////
1219  }// end routine
1220 //////////////////////////////////////////////////////////////////////////////
1221 
1222 
1223 //////////////////////////////////////////////////////////////////////////////
1224 ///////////////////////////////////////////////////////////////////////////cc
1225 //////////////////////////////////////////////////////////////////////////////
1226 /// Glenn's Rhart exit
1227 //////////////////////////////////////////////////////////////////////////////
1228 void CP_Rho_RHartExt::exitForDebugging(){
1229  countDebug++;
1230  if(countDebug==(rhoRsubplanes*ngridc*nchareHartAtmT)){
1231  countDebug=0;
1232  CkPrintf("I am in the exitfordebuging rhorhartext puppy. Bye-bye\n");
1233  CkExit();
1234  }//endif
1235 }
1236 //////////////////////////////////////////////////////////////////////////////
1237 /*@}*/
eesCache::registerCacheRHart
void registerCacheRHart(int)
= RhoRhart Cache Management tool
Definition: eesCache.C:158
UberCollection
holds the UberIndex and the offset for proxies
Definition: Uber.h:68
CP_Rho_RHartExt::startEextIter
void startEextIter()
Invoke by Rspace-density : Density has arrived in r-space and will soon arrive in g-space...
Definition: CP_Rho_RHartExt.C:287
FFTcache::doEextFFTGxToRx_Rchare
void doEextFFTGxToRx_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward //////////////////////////////////////////////...
Definition: fftCache.C:1075
eesCache::queryCacheRHart
void queryCacheRHart(int, int, int)
= RhoGhart Cache Management tool
Definition: eesCache.C:438
CP_Rho_RHartExt::sendAtmSfRyToGy
void sendAtmSfRyToGy()
Double Transpose Fwd Send : A(gx,y,z) on the way to A(gx,gy,z) Send so that (y,z) parallelism is swit...
Definition: CP_Rho_RHartExt.C:492
CP_Rho_RHartExt::exitForDebugging
void exitForDebugging()
Glenn's Rhart exit.
Definition: CP_Rho_RHartExt.C:1228
FFTcache::doEextFFTGtoR_Rchare
void doEextFFTGtoR_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward //////////////////////////////////////////////...
Definition: fftCache.C:1028
CP_Rho_RHartExt::computeAtmForc
void computeAtmForc(int)
Get forces on atoms from Ewald or electrons /////////////////////////////////////////////////////////...
Definition: CP_Rho_RHartExt.C:1107
FFTcache::doEextFFTRtoG_Rchare
void doEextFFTRtoG_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward /////////////////////////////////////////////...
Definition: fftCache.C:923
CP_Rho_RHartExt::fftAtmSfRtoG
void fftAtmSfRtoG()
FFT of SF(x,y,z,iatmTyp) -> SF(gx,gy,z,iatmTyp)
Definition: CP_Rho_RHartExt.C:417
fftCacheSlab.h
Add type declarations for simulationConstants class (readonly vars) and once class for each type of o...
CP_Rho_RHartExt::init
void init()
Definition: CP_Rho_RHartExt.C:143
CP_Rho_RHartExt::sendAtmSfRhoGHart
void sendAtmSfRhoGHart()
Send SF(gx,gy,z,iatmTYP) to g-space whence the FFT will be completed.
Definition: CP_Rho_RHartExt.C:654
CP_Rho_RHartExt::computeAtmSF
void computeAtmSF()
Start the real space part of the EES interpolation for atmSF(iatmTyp) ///////////////////////////////...
Definition: CP_Rho_RHartExt.C:363
cpaimd.h
Some basic data structures and the array map classes are defined here.
CP_Rho_RHartExt::fftAtmForcGtoR
void fftAtmForcGtoR(int flagEwd)
Complete the FFT : sf(gx,gy,z) ->sf(x,y,z) //////////////////////////////////////////////////////////...
Definition: CP_Rho_RHartExt.C:878
FFTcache::doEextFFTRxToGx_Rchare
void doEextFFTRxToGx_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward /////////////////////////////////////////////...
Definition: fftCache.C:969
FFTcache::doEextFFTGyToRy_Rchare
void doEextFFTGyToRy_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(gx,gy,z) -> data(x,y,z) : forward //////////////////////////////////////////////...
Definition: fftCache.C:1106
CP_Rho_RHartExt::recvAtmForcGxToRx
void recvAtmForcGxToRx(RhoGHartMsg *msg)
Definition: CP_Rho_RHartExt.C:1021
CP_State_Plane.h
eesCache
Group Container class : Only allowed chare data classes have data.
Definition: eesCache.h:18
debug_flags.h
Useful debugging flags.
CP_Rho_RHartExt::registrationDone
void registrationDone(CkReductionMsg *msg)
= Make sure everyone is registered on the 1st time step
Definition: CP_Rho_RHartExt.C:333
CP_Rho_RHartExt::recvAtmSfRyToGy
void recvAtmSfRyToGy(RhoGHartMsg *msg)
Double Transpose Fwd Recv : A(gx,y,z) on the way to A(gx,gy,z) Recv so that (y,z) parallel switched t...
Definition: CP_Rho_RHartExt.C:577
CP_Rho_RHartExt::recvAtmForcFromRhoGHart
void recvAtmForcFromRhoGHart(RhoRHartMsg *msg)
Hartree sends back atom forces from e-atm interation Depending on the flag, it is Ewald or e-atm inte...
Definition: CP_Rho_RHartExt.C:753
FFTcache::doEextFFTRyToGy_Rchare
void doEextFFTRyToGy_Rchare(complex *, double *, int, int, int, int)
Eext : Rchare : data(x,y,z) -> data(gx,gy,z) : Backward /////////////////////////////////////////////...
Definition: fftCache.C:1001
CP_Rho_RHartExt::sendAtmForcGxToRx
void sendAtmForcGxToRx(int iopt)
Definition: CP_Rho_RHartExt.C:930