OpenAtom  Version1.5a
stateSlab.C
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////=
2 ///////////////////////////////////////////////////////////////////////////////c
3 ///////////////////////////////////////////////////////////////////////////////=
4 /** \file stateSlab.C
5  * Add functions to allow application programmers to initialize these and
6  * the corresponding functions in Charm++ to call these functions with
7  * appropriate parameters
8  */
9 ///////////////////////////////////////////////////////////////////////////////=
10 
11 #include "utility/util.h"
12 #include "debug_flags.h"
13 #include <cmath>
14 #include "main/cpaimd.h"
15 #include "fft_slab_ctrl/fftCacheSlab.h"
16 #include "cp_state_ctrl/CP_State_Plane.h"
17 #include "main/eesCache.h"
18 
19 ///////////////////////////////////////////////////////////////////////////////=
20 
21 extern CkVec <CProxy_FFTcache> UfftCacheProxy;
22 extern Config config;
23 extern CkVec <CProxy_eesCache> UeesCacheProxy;
24 extern int nstates;
25 extern int sizeX;
26 
27 ///////////////////////////////////////////////////////////////////////////////=
28 
29 
30 ///////////////////////////////////////////////////////////////////////////////=
31 ///////////////////////////////////////////////////////////////////////////////c
32 ///////////////////////////////////////////////////////////////////////////////=
33 /* This gets called at the end of the GStatePlane constructor */
34 ///////////////////////////////////////////////////////////////////////////////=
35 void initGStateSlab(GStateSlab *gs, int sizeX, int sizeY, int sizeZ, int gSpaceUnits,
36  int realSpaceUnits, int s_grain, int iplane_ind,int istate_ind,
37  int len_nhc_cp, int num_nhc_cp, int nck_nhc_cp)
38 ///////////////////////////////////////////////////////////////////////////////=
39  {//begin routine
40 ///////////////////////////////////////////////////////////////////////////////=
41 /// Explanation of organization of data: A point is psi(kx,ky,kz)
42 ///
43 /// Each gStateSlab is a collection of pts in lines of constant kx,ky.
44 /// The maximum number of pts in any line is size[1] = nfftz
45 /// Some lines are longer/shorter to spherical trunction (|k|<k_cut).
46 /// Also, kx>=0 when at the Gamma point or doublePack=1, the only
47 /// thing that has been recently tested. Collections of lines are parallelized.
48 /// Each collection does not correspond to a unique plane of kx.
49 /// The parameter gSpaceUnits is obsolete and must be unity.
50 /// In order to create the state in real space, psi(x,y,z), an fft
51 /// of all kz lines is performed followed by a transpose. The result of the transpose
52 /// is parallelized by planes of z. This is different than it used to be.
53 /// The number of points AFTER the fft is nlines*nfftz will be different
54 /// for each chare array
55 
56  if(gSpaceUnits!=1 || realSpaceUnits!=1){
57  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
58  CkPrintf("gspacePPC==1 obsolete. We no longer parallelize g-space by plane\n");
59  CkPrintf("realSpacePPC==1 although real space is parallelized by plane\n");
60  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
61  CkExit();
62  }//endif
63  gs->numNonZeroPlanes=1;
64  gs->mysizeX = sizeX;
65  gs->planeSize = size2d(sizeY,sizeZ); // fftsizes (sizeY,sizeZ) sizeX is a global
66  gs->numPoints = 0; // number of packed points (data size before fft)
67  gs->numRuns = 0; // 2*(number of lines) in the collection
68  gs->numLines = 0; // (number of lines) in the collection
69  gs->fftReqd = false; // false if this chare has no packed pts
70  gs->numFull = 0; // number of pts : numLines*nfftz
71 
72  gs->S_grainSize = s_grain; // PC grainsize
73  gs->xdim = 1; // may need some love
74  gs->ydim = gs->planeSize[0];
75  gs->zdim = gs->planeSize[1];
76  gs->iplane_ind = iplane_ind;
77  gs->istate_ind = istate_ind;
78  gs->initNHC(len_nhc_cp,num_nhc_cp,nck_nhc_cp);
79 
80 ///////////////////////////////////////////////////////////////////////////////=
81  }//end routine
82 ///////////////////////////////////////////////////////////////////////////////=
83 
84 
85 ///////////////////////////////////////////////////////////////////////////////=
86 ///////////////////////////////////////////////////////////////////////////////c
87 ///////////////////////////////////////////////////////////////////////////////=
88 GStateSlab::~GStateSlab() {
89 
90  if(packedPlaneData !=NULL) fftw_free( packedPlaneData);
91  if(packedForceData !=NULL) fftw_free( packedForceData);
92  if(packedVelData !=NULL) fftw_free( packedVelData);
93  if(cp_min_opt==0){
94  if(packedPlaneDataScr !=NULL) fftw_free( packedPlaneDataScr);
95  }//endif
96 #ifdef _CP_DEBUG_UPDATE_OFF_
97  if(cp_min_opt==1){
98  if(packedPlaneDataTemp!=NULL) fftw_free( packedPlaneDataTemp);
99  }//endif
100 #endif
101  destroyNHC();
102 
103  packedPlaneData = NULL;
104  packedForceData = NULL;
105  packedVelData = NULL;
106  packedPlaneDataScr = NULL;
107  packedPlaneDataTemp = NULL;
108 
109 }
110 ///////////////////////////////////////////////////////////////////////////////=
111 
112 ///////////////////////////////////////////////////////////////////////////////=
113 ///////////////////////////////////////////////////////////////////////////////c
114 ///////////////////////////////////////////////////////////////////////////////=
115 void GStateSlab::pup(PUP::er &p) {
116 ///////////////////////////////////////////////////////////////////////////////=
117 /// Dont have to pup fftw plans - they live in the fft cache group
118 
119 /// CkPrintf("gs pup\n:");
120  p|cp_min_opt;
121  p|ees_nonlocal;
122  p|numNonZeroPlanes;
123  p|numRuns;
124  p|numLines;
125  p|numPoints;
126  p|numFull;
127  p|mysizeX;
128  p|planeSize;
129  p|fftReqd;
130  p|S_grainSize;
131  p|xdim;
132  p|ydim;
133  p|zdim;
134  p|ngridaNL;
135  p|ngridbNL;
136  p|ngridcNL;
137  p|iplane_ind;
138  p|istate_ind;
139  p|ihave_kx0;
140  p|ihave_g000;
141  p|ind_g000;
142  p|kx0_strt;
143  p|kx0_end;
144  p|nkx0; p|nkx0_uni; p|nkx0_red; p|nkx0_zero;
145  p|eke_ret;
146  p|fictEke_ret;
147  p|ekeNhc_ret;
148  p|potNHC_ret;
149  p|degfree;
150  p|degfreeNHC;
151 
152  if (p.isUnpacking()) {
153  packedPlaneData = (complex *)fftw_malloc(numPoints*sizeof(complex));
154  packedForceData = (complex *)fftw_malloc(numFull*sizeof(complex));
155  packedVelData = (complex *)fftw_malloc(numPoints*sizeof(complex));
156  packedRedPsi = (complex *)fftw_malloc(nkx0*sizeof(complex));
157  if(cp_min_opt==0){
158  packedPlaneDataScr = (complex *)fftw_malloc(numPoints*sizeof(complex));
159  }//endif
160 #ifdef _CP_DEBUG_UPDATE_OFF_
161  if(cp_min_opt==1){
162  packedPlaneDataTemp = (complex *)fftw_malloc(numPoints*sizeof(complex));
163  }//endif
164 #endif
165  }//endif
166  p((char *) packedPlaneData, numPoints*sizeof(complex));
167  p((char *) packedForceData, numFull*sizeof(complex));
168  p((char *) packedVelData, numPoints*sizeof(complex)); ///g under min
169  p((char *) packedRedPsi, nkx0*sizeof(complex));
170  if(cp_min_opt==0){
171  p((char *) packedPlaneDataScr, numPoints*sizeof(complex));
172  }//endif
173 #ifdef _CP_DEBUG_UPDATE_OFF_
174  if(cp_min_opt==1){
175  p((char *) packedPlaneDataTemp, numPoints*sizeof(complex));
176  }//endif
177 #endif
178 
179  p|nck_nhc_cp;
180  p|len_nhc_cp;
181  p|num_nhc_cp;
182  p|kTCP;
183  p|tauNHCCP;
184  if (p.isUnpacking()) {
185  initNHC(len_nhc_cp,num_nhc_cp,nck_nhc_cp);
186  }//endif recving
187  int nsize = num_nhc_cp*len_nhc_cp*nck_nhc_cp;
188  double *xt = new double[nsize];
189  double *xtp = new double[nsize];
190  double *vt = new double[nsize];
191  double *ft = new double[nsize];
192  if(p.isPacking()){
193  int iii=0;
194  for(int k =0;k<num_nhc_cp;k++){
195  for(int i =0;i<num_nhc_cp;i++){
196  for(int j =0;j<len_nhc_cp;j++){
197  xt[iii] = xNHC[k][i][j];
198  xtp[iii] = xNHCP[k][i][j];
199  vt[iii] = vNHC[k][i][j];
200  ft[iii] = fNHC[k][i][j];
201  iii++;
202  }}}
203  }//endif sending
204  p(xt,nsize);
205  p(xtp,nsize);
206  p(vt,nsize);
207  p(ft,nsize);
208  p(degFreeSplt,nck_nhc_cp);
209  p(istrNHC,nck_nhc_cp);
210  p(iendNHC,nck_nhc_cp);
211  p(mNHC,len_nhc_cp);
212  p(v0NHC,num_nhc_cp);
213  p(a2NHC,num_nhc_cp);
214  p(a4NHC,num_nhc_cp);
215  if (p.isUnpacking()) {
216  int iii=0;
217  for(int k =0;k<nck_nhc_cp;k++){
218  for(int i =0;i<num_nhc_cp;i++){
219  for(int j =0;j<len_nhc_cp;j++){
220  xNHC[k][i][j] = xt[iii];
221  xNHCP[k][i][j] = xtp[iii];
222  vNHC[k][i][j] = vt[iii];
223  fNHC[k][i][j] = ft[iii];
224  iii++;
225  }}}
226  }//endif receiving
227  delete []xt;
228  delete []xtp;
229  delete []vt;
230  delete []ft;
231 
232 ///////////////////////////////////////////////////////////////////////////////=
233  }//end routine
234 ///////////////////////////////////////////////////////////////////////////////=
235 
236 
237 
238 
239 ///////////////////////////////////////////////////////////////////////////////=
240 ///////////////////////////////////////////////////////////////////////////////c
241 ///////////////////////////////////////////////////////////////////////////////=
242 void GStateSlab::addForces(complex *points, const int *k_x){
243 ///////////////////////////////////////////////////////////////////////////////=
244  int i;
245  double wght;
246 
247  if(!config.doublePack){
248 
249  for(i = 0; i < numPoints; i++){
250  packedForceData[i] += points[i];
251  }//endfor
252 
253  }else{
254 
255  int nfreq=1000;
256  for(i = 0; i < nkx0; i++){
257 #ifdef _CP_DEBUG_VKS_OFF_ // only non-local, no vks forces
258  packedForceData[i].re = 0.0; packedForceData[i].im = 0.0;
259 #endif
260  packedForceData[i] += points[i];
261  }//endfor
262 
263  for(i = nkx0; i < numPoints; i++){
264 #ifdef _CP_DEBUG_VKS_OFF_ // only non-local, no vks forces
265  packedForceData[i].re = 0.0; packedForceData[i].im = 0.0;
266 #endif
267  packedForceData[i] *= 2.0;
268  packedForceData[i] += points[i];
269  }//endfor
270 
271  }//endif : doublePack
272 
273 ///////////////////////////////////////////////////////////////////////////////=
274  }//end routine
275 ///////////////////////////////////////////////////////////////////////////////=
276 
277 
278 ///////////////////////////////////////////////////////////////////////////////=
279 ///////////////////////////////////////////////////////////////////////////////c
280 ///////////////////////////////////////////////////////////////////////////////=
281 /*
282  * In G-space our representation uses wrapping of coordinates.
283  * When reading from the files, if x < 0, x += sizeX;
284  */
285 ///////////////////////////////////////////////////////////////////////////////=
286 
287 void GStateSlab::setKRange(int n, int *k_x, int *k_y, int *k_z){
288 
289 ///////////////////////////////////////////////////////////////////////=
290 /// Construct the k-vectors
291 
292  CkAssert(n == numPoints);
293 
294 ///////////////////////////////////////////////////////////////////////=
295 /// Find pts with k_x==0:
296 
297  int i;
298  ihave_g000 = 0;
299  ind_g000 = -1;
300  ihave_kx0 = 0;
301  nkx0 = 0;
302  nkx0_uni = 0;
303  nkx0_red = 0;
304  nkx0_zero = 0;
305  kx0_strt = 0;
306  kx0_end = 0;
307 
308  if(config.doublePack){
309  for(i=0;i<numPoints;i++){
310  if(k_x[i]==0 && k_y[i]>0){nkx0_uni++;}
311  if(k_x[i]==0 && k_y[i]<0){nkx0_red++;}
312  if(k_x[i]==0 && k_y[i]==0 && k_z[i]>=0){nkx0_uni++;}
313  if(k_x[i]==0 && k_y[i]==0 && k_z[i]<0){nkx0_red++;}
314  if(k_x[i]==0 && k_y[i]==0 && k_z[i]==0){nkx0_zero++;ihave_g000=1;ind_g000=i;}
315  if(k_x[i]==0){
316  if(ihave_kx0==0){kx0_strt=i;}
317  ihave_kx0=1;
318  nkx0++;
319  }//endif
320  }//endfor
321  kx0_end = kx0_strt + nkx0;
322  }else{
323  for(int i=0;i<numPoints;i++){
324  if(k_x[i]==0 && k_y[i]==0 && k_z[i]==0){ihave_g000=1;ind_g000=i;}
325  }//endfor
326  }//endif
327 
328 ///////////////////////////////////////////////////////////////////////=
329 /// Check the layout for doublePack case: kx=0 first
330 
331  if(config.doublePack){
332 
333  if(kx0_strt!=0){
334  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
335  CkPrintf("kx=0 should be stored first | kx_srt !=0\n");
336  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
337  CkAbort("kx=0 should be stored first | kx_srt !=0\n");
338  }//endif
339 
340  if(nkx0!=nkx0_uni+nkx0_red){
341  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
342  CkPrintf("Incorrect count of redundant guys\n");
343  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
344  CkAbort("Incorrect count of redundant guys\n");
345  }//endif
346 
347  for(i=0;i<nkx0;i++){
348  if(k_x[i]!=0){
349  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
350  CkPrintf("kx should be stored consecutively and first\n");
351  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
352  CkAbort("kx should be stored consecutively and first\n");
353  }//endif
354  }//endif
355 
356  for(i=0;i<nkx0_red;i++){
357  if(k_y[i]>0 || (k_y[i]==0 && k_z[i]>=0)){
358  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
359  CkPrintf("ky <0 should be stored first\n");
360  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
361  CkAbort("ky <0 should be stored first\n");
362  }//endif
363  }//endfor
364 
365  for(i=nkx0_red;i<nkx0_uni;i++){
366  if(k_y[i]<0 || (k_y[i]==0 && k_z[i]<0)){
367  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
368  CkPrintf("ky <0 should be stored first\n");
369  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
370  CkAbort("ky <0 should be stored first\n");
371  }//endif
372  }//endfor
373 
374  }//endif
375 
376  packedRedPsi = (complex *)fftw_malloc(nkx0*sizeof(complex));
377  packedRedPsiV = (complex *)fftw_malloc(nkx0*sizeof(complex));
378  memset(packedRedPsi, 0, sizeof(complex)*nkx0);
379  memset(packedRedPsiV, 0, sizeof(complex)*nkx0);
380 
381 ///////////////////////////////////////////////////////////////////////////////=
382  }//end routine
383 ///////////////////////////////////////////////////////////////////////////////=
384 
385 
386 
387 ///////////////////////////////////////////////////////////////////////////////=
388 ///////////////////////////////////////////////////////////////////////////////c
389 ///////////////////////////////////////////////////////////////////////////////=
390 /* This gets called at the end of the CP_State_RealSpacePlane constructor */
391 ///////////////////////////////////////////////////////////////////////////////=
392 void initRealStateSlab(RealStateSlab *rs, int ngrid_a,int ngrid_b, int ngrid_c,
393  int gSpaceUnits, int realSpaceUnits, int stateIndex, int planeIndex)
394 ///////////////////////////////////////////////////////////////////////////////=
395  {//begin routine
396 ///////////////////////////////////////////////////////////////////////////////=
397 /// Explanation of the organization:
398 /// Each CP_State_RealSpacePlane instance is actually a planes of psi_I(x,y,z).
399 /// We index into planeArr using the z-coordinate
400 /// and each chare is an x-y plane.
401 /// Data arrives from the gspaceplane of total size : nlines_tot*nfftz
402 
403  if(gSpaceUnits!=1 || realSpaceUnits!=1){
404  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
405  CkPrintf("gspacePPC==1 obsolete. We no longer parallelize g-space by plane\n");
406  CkPrintf("realSpacePPC==1 although real space is parallelized by plane\n");
407  CkPrintf("@@@@@@@@@@@@@@@@@@@@_error_@@@@@@@@@@@@@@@@@@@@\n");
408  CkExit();
409  }//endif
410 
411  rs->thisState = stateIndex; // my state (I)
412  rs->thisPlane = planeIndex; // my plane (z)
413  CPcharmParaInfo *sim = CPcharmParaInfo::get();
414  rs->numPlanesToExpect = sim->nchareG;
415  if(config.doublePack){
416  int rsize_a = ngrid_a*(ngrid_b/2+1);
417  int rsize_b = ngrid_b*(ngrid_a/2+1);
418  rs->rsize = (rsize_a > rsize_b ? rsize_a : rsize_b);
419  rs->nsize = ngrid_a*ngrid_b; // when fft is completed
420  rs->size = rs->rsize;
421  }else{
422  rs->rsize = 2*ngrid_a*ngrid_b;
423  rs->nsize = ngrid_a*ngrid_b; // when fft is complete
424  rs->size = rs->nsize;
425  }//endif
426 
427  rs->planeArr = (complex *) fftw_malloc(rs->size * sizeof(complex));
428  rs->planeArrR = reinterpret_cast<double*> (rs->planeArr);
429 
430  rs->ngrid_a = ngrid_a;
431  rs->ngrid_b = ngrid_b;
432 
433 ///////////////////////////////////////////////////////////////////////////////=
434  }//end routine
435 ///////////////////////////////////////////////////////////////////////////////=
436 
437 
438 ///////////////////////////////////////////////////////////////////////////////=
439 ///////////////////////////////////////////////////////////////////////////////c
440 ///////////////////////////////////////////////////////////////////////////////=
441 RealStateSlab::~RealStateSlab() {
442 ///////////////////////////////////////////////////////////////////////////////=
443 
444  if(planeArr != NULL) {fftw_free(planeArr);planeArr = NULL; planeArrR=NULL;}
445 
446 }
447 ///////////////////////////////////////////////////////////////////////////////=
448 
449 
450 ///////////////////////////////////////////////////////////////////////////////=
451 ///////////////////////////////////////////////////////////////////////////////c
452 ///////////////////////////////////////////////////////////////////////////////=
453 void RealStateSlab::pup(PUP::er &p) {
454 ///////////////////////////////////////////////////////////////////////////////=
455  p|ngrid_a;
456  p|ngrid_b;
457  p|size;
458  if (p.isUnpacking()) {
459  planeArr = (complex *) fftw_malloc(size*sizeof(complex));
460  planeArrR = reinterpret_cast<double*> (planeArr);
461  }
462  PUParray(p, planeArr, size);
463  p|thisState;
464  p|thisPlane;
465  p|numPlanesToExpect;
466  p|nsize;
467  p|rsize;
468  p|size;
469  p|e_gga;
470 }
471 ///////////////////////////////////////////////////////////////////////////////=
472 
473 
474 ///////////////////////////////////////////////////////////////////////////////=
475 ///////////////////////////////////////////////////////////////////////////////c
476 ///////////////////////////////////////////////////////////////////////////////=
477 void RealStateSlab::allocate() {
478  if(planeArr == NULL) {
479  planeArr = (complex *) fftw_malloc(size*sizeof(complex));
480  planeArrR = reinterpret_cast<double*> (planeArr);
481  }//endif
482 }
483 ///////////////////////////////////////////////////////////////////////////////=
484 
485 
486 ///////////////////////////////////////////////////////////////////////////////=
487 ///////////////////////////////////////////////////////////////////////////////c
488 ///////////////////////////////////////////////////////////////////////////////=
489 void RealStateSlab::destroy() {
490  if(planeArr != NULL) {fftw_free(planeArr); planeArr=NULL; planeArrR = NULL;}
491 }
492 ///////////////////////////////////////////////////////////////////////////////=
493 
494 
GStateSlab::pup
void pup(PUP::er &)
= data
Definition: stateSlab.C:115
size2d
== Size or location in a regular 2D array
Definition: util.h:59
GStateSlab::setKRange
void setKRange(int, int *, int *, int *)
Definition: stateSlab.C:287
initGStateSlab
void initGStateSlab(GStateSlab *gs, int sizeX, int sizeY, int sizeZ, int gSpaceUnits, int realSpaceUnits, int s_grain, int iplane_ind, int istate_ind, int len_nhc_cp, int num_nhc_cp, int nck_nhc_cp)
= slab initialization helpers
Definition: stateSlab.C:35
fftCacheSlab.h
Add type declarations for simulationConstants class (readonly vars) and once class for each type of o...
cpaimd.h
Some basic data structures and the array map classes are defined here.
CP_State_Plane.h
debug_flags.h
Useful debugging flags.
initRealStateSlab
void initRealStateSlab(RealStateSlab *rs, int ngrid_a, int ngrid_b, int ngrid_c, int gSpaceUnits, int realSpaceUnits, int stateIndex, int planeIndex)
Definition: stateSlab.C:392