1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 This file is from LAMMPS
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include "lmptype.h"
47 #include <mpi.h>
48 #include <string.h>
49 #include "write_data.h"
50 #include "atom.h"
51 #include "atom_vec.h"
52 #include "group.h"
53 #include "force.h"
54 #include "pair.h"
55 #include "bond.h"
56 #include "angle.h"
57 #include "dihedral.h"
58 #include "improper.h"
59 #include "update.h"
60 #include "modify.h"
61 #include "fix.h"
62 #include "domain.h"
63 #include "universe.h"
64 #include "comm.h"
65 #include "output.h"
66 #include "thermo.h"
67 #include "memory.h"
68 #include "error.h"
69
70 using namespace LAMMPS_NS;
71
72 enum{IGNORE,WARN,ERROR}; // same as thermo.cpp
73 enum{II,IJ};
74
75 /* ---------------------------------------------------------------------- */
76
WriteData(LAMMPS * lmp)77 WriteData::WriteData(LAMMPS *lmp) : Pointers(lmp)
78 {
79 MPI_Comm_rank(world,&me);
80 MPI_Comm_size(world,&nprocs);
81 }
82
83 /* ----------------------------------------------------------------------
84 called as write_data command in input script
85 ------------------------------------------------------------------------- */
86
command(int narg,char ** arg)87 void WriteData::command(int narg, char **arg)
88 {
89 if (domain->box_exist == 0)
90 error->all(FLERR,"Write_data command before simulation box is defined");
91
92 if (narg < 1) error->all(FLERR,"Illegal write_data command");
93
94 // if filename contains a "*", replace with current timestep
95
96 char *ptr;
97 int n = strlen(arg[0]) + 16;
98 char *file = new char[n];
99
100 if ((ptr = strchr(arg[0],'*'))) {
101 *ptr = '\0';
102 sprintf(file,"%s" BIGINT_FORMAT "%s",arg[0],update->ntimestep,ptr+1);
103 } else strcpy(file,arg[0]);
104
105 // read optional args
106
107 pairflag = II;
108
109 tag_offset = tag_max = 0;
110 bool tag_flag = false;
111
112 int iarg = 1;
113 while (iarg < narg) {
114 if (strcmp(arg[iarg],"pair") == 0) {
115 if (iarg+2 > narg) error->all(FLERR,"Illegal write_data command");
116 if (strcmp(arg[iarg+1],"ii") == 0) pairflag = II;
117 else if (strcmp(arg[iarg+1],"ij") == 0) pairflag = IJ;
118 else error->all(FLERR,"Illegal write_data command");
119 iarg += 2;
120 } else if (strcmp(arg[iarg],"tag_offset") == 0) {
121 if (iarg+2 > narg) error->all(FLERR,"Illegal write_data command");
122 tag_offset = atoi(arg[iarg+1]);
123 printf("Applying a tag offset of %d to atom data\n",tag_offset);
124 tag_flag = true;
125 iarg += 2;
126 } else error->all(FLERR,"Illegal write_data command");
127 }
128
129 // init entire system since comm->exchange is done
130 // comm::init needs neighbor::init needs pair::init needs kspace::init, etc
131
132 if (comm->me == 0 && screen)
133 fprintf(screen,"System init for write_data ...\n");
134 lmp->init();
135
136 // move atoms to new processors before writing file
137 // do setup_pre_exchange to force update of per-atom info if needed
138 // enforce PBC in case atoms are outside box
139 // call borders() to rebuild atom map since exchange() destroys map
140
141 modify->setup_pre_exchange();
142 if (domain->triclinic) domain->x2lamda(atom->nlocal);
143 domain->pbc();
144 domain->reset_box();
145 comm->setup();
146 comm->exchange();
147 comm->borders();
148 if (domain->triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
149 modify->forceMeshExchange(); // call explicit exchange of mesh data
150
151 write(file);
152
153 delete [] file;
154
155 if(tag_flag && 0 == comm->me)
156 {
157 FILE *fpt = fopen("max_tag","w");
158 fprintf(fpt,"%d\n",tag_max);
159 fclose(fpt);
160 }
161 }
162
163 /* ----------------------------------------------------------------------
164 called from command()
165 might later let it be directly called within run/minimize loop
166 ------------------------------------------------------------------------- */
167
write(char * file)168 void WriteData::write(char *file)
169 {
170 // special case where reneighboring is not done in integrator
171 // on timestep data file is written (due to build_once being set)
172 // if box is changing, must be reset, else data file will have
173 // wrong box size and atoms will be lost when data file is read
174 // other calls to pbc and domain and comm are not made,
175 // b/c they only make sense if reneighboring is actually performed
176
177 //if (neighbor->build_once) domain->reset_box();
178
179 // natoms = sum of nlocal = value to write into data file
180 // if unequal and thermo lostflag is "error", don't write data file
181
182 bigint nblocal = atom->nlocal;
183 bigint natoms;
184 MPI_Allreduce(&nblocal,&natoms,1,MPI_LMP_BIGINT,MPI_SUM,world);
185 if (natoms != atom->natoms && output->thermo->lostflag == ERROR)
186 error->all(FLERR,"Atom count is inconsistent, cannot write data file");
187
188 // sum up bond,angle counts
189 // may be different than atom->nbonds,nangles if broken/turned-off
190
191 if (atom->nbonds || atom->nbondtypes) {
192 nbonds_local = atom->avec->pack_bond(NULL);
193 MPI_Allreduce(&nbonds_local,&nbonds,1,MPI_LMP_BIGINT,MPI_SUM,world);
194 }
195 if (atom->nangles || atom->nangletypes) {
196 nangles_local = atom->avec->pack_angle(NULL);
197 MPI_Allreduce(&nangles_local,&nangles,1,MPI_LMP_BIGINT,MPI_SUM,world);
198 }
199
200 // open data file
201
202 if (me == 0) {
203 fp = fopen(file,"w");
204 if (fp == NULL) {
205 char str[512];
206 sprintf(str,"Cannot open data file %s",file);
207 error->one(FLERR,str);
208 }
209 }
210
211 // proc 0 writes header, ntype-length arrays, force fields
212
213 if (me == 0) {
214 header();
215 type_arrays();
216 force_fields();
217 }
218
219 // per atom info
220
221 if (natoms) atoms();
222 if (natoms) velocities();
223 if (atom->nbonds && nbonds) bonds();
224 if (atom->nangles && nangles) angles();
225 if (atom->ndihedrals) dihedrals();
226 if (atom->nimpropers) impropers();
227
228 // extra sections managed by fixes
229
230 for (int i = 0; i < modify->nfix; i++)
231 if (modify->fix[i]->wd_section)
232 for (int m = 0; m < modify->fix[i]->wd_section; m++) fix(i,m);
233
234 tag_max = static_cast<bigint>(atom->tag_max()) + tag_offset;
235
236 // close data file
237
238 if (me == 0) fclose(fp);
239 }
240
241 /* ----------------------------------------------------------------------
242 proc 0 writes out data file header
243 ------------------------------------------------------------------------- */
244
header()245 void WriteData::header()
246 {
247 fprintf(fp,"LAMMPS data file via write_data, version %s, "
248 "timestep = " BIGINT_FORMAT "\n",
249 universe->version,update->ntimestep);
250
251 fprintf(fp,"\n");
252
253 fprintf(fp,BIGINT_FORMAT " atoms\n",atom->natoms);
254 fprintf(fp,"%d atom types\n",atom->ntypes);
255
256 if (atom->nbonds || atom->nbondtypes) {
257 fprintf(fp,BIGINT_FORMAT " bonds\n",nbonds);
258 fprintf(fp,"%d bond types\n",atom->nbondtypes);
259 }
260 if (atom->nangles || atom->nangletypes) {
261 fprintf(fp,BIGINT_FORMAT " angles\n",nangles);
262 fprintf(fp,"%d angle types\n",atom->nangletypes);
263 }
264 if (atom->ndihedrals || atom->ndihedraltypes) {
265 fprintf(fp,BIGINT_FORMAT " dihedrals\n",atom->ndihedrals);
266 fprintf(fp,"%d dihedral types\n",atom->ndihedraltypes);
267 }
268 if (atom->nimpropers || atom->nimpropertypes) {
269 fprintf(fp,BIGINT_FORMAT " impropers\n",atom->nimpropers);
270 fprintf(fp,"%d improper types\n",atom->nimpropertypes);
271 }
272
273 for (int i = 0; i < modify->nfix; i++)
274 if (modify->fix[i]->wd_header)
275 for (int m = 0; m < modify->fix[i]->wd_header; m++)
276 modify->fix[i]->write_data_header(fp,m);
277
278 fprintf(fp,"\n");
279
280 fprintf(fp,"%-1.16e %-1.16e xlo xhi\n",domain->boxlo[0],domain->boxhi[0]);
281 fprintf(fp,"%-1.16e %-1.16e ylo yhi\n",domain->boxlo[1],domain->boxhi[1]);
282 fprintf(fp,"%-1.16e %-1.16e zlo zhi\n",domain->boxlo[2],domain->boxhi[2]);
283
284 if (domain->triclinic)
285 fprintf(fp,"%-1.16e %-1.16e %-1.16e xy xz yz\n",
286 domain->xy,domain->xz,domain->yz);
287 }
288
289 /* ----------------------------------------------------------------------
290 proc 0 writes out any type-based arrays that are defined
291 ------------------------------------------------------------------------- */
292
type_arrays()293 void WriteData::type_arrays()
294 {
295 if (atom->mass) {
296 double *mass = atom->mass;
297 fprintf(fp,"\nMasses\n\n");
298 for (int i = 1; i <= atom->ntypes; i++) fprintf(fp,"%d %g\n",i,mass[i]);
299 }
300 }
301
302 /* ----------------------------------------------------------------------
303 proc 0 writes out force field info
304 ------------------------------------------------------------------------- */
305
force_fields()306 void WriteData::force_fields()
307 {
308 if (force->pair && force->pair->writedata) {
309 if (pairflag == II) {
310 fprintf(fp,"\nPair Coeffs\n\n");
311 force->pair->write_data(fp);
312 } else if (pairflag == IJ) {
313 fprintf(fp,"\nPairIJ Coeffs\n\n");
314 force->pair->write_data_all(fp);
315 }
316 }
317 if (atom->avec->bonds_allow && force->bond && force->bond->writedata) {
318 fprintf(fp,"\nBond Coeffs\n\n");
319 force->bond->write_data(fp);
320 }
321 if (atom->avec->angles_allow && force->angle && force->angle->writedata) {
322 fprintf(fp,"\nAngle Coeffs\n\n");
323 force->angle->write_data(fp);
324 }
325 if (atom->avec->dihedrals_allow && force->dihedral &&
326 force->dihedral->writedata) {
327 fprintf(fp,"\nDihedral Coeffs\n\n");
328 force->dihedral->write_data(fp);
329 }
330 if (atom->avec->impropers_allow && force->improper &&
331 force->improper->writedata) {
332 fprintf(fp,"\nImproper Coeffs\n\n");
333 force->improper->write_data(fp);
334 }
335 }
336
337 /* ----------------------------------------------------------------------
338 write out Atoms section of data file
339 ------------------------------------------------------------------------- */
340
atoms()341 void WriteData::atoms()
342 {
343 // communication buffer for all my Atom info
344 // max_size = largest buffer needed by any proc
345
346 int ncol = atom->avec->size_data_atom + 3;
347
348 int sendrow = atom->nlocal;
349 int maxrow;
350 MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
351
352 double **buf;
353 if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
354 else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
355
356 // pack my atom data into buf
357
358 atom->avec->pack_data(buf,tag_offset);
359
360 // write one chunk of atoms per proc to file
361 // proc 0 pings each proc, receives its chunk, writes to file
362 // all other procs wait for ping, send their chunk to proc 0
363
364 int tmp,recvrow;
365 MPI_Status status;
366 MPI_Request request;
367
368 if (me == 0) {
369 fprintf(fp,"\nAtoms\n\n");
370 for (int iproc = 0; iproc < nprocs; iproc++) {
371 if (iproc) {
372 MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request);
373 MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
374 MPI_Wait(&request,&status);
375 MPI_Get_count(&status,MPI_DOUBLE,&recvrow);
376 recvrow /= ncol;
377 } else recvrow = sendrow;
378
379 atom->avec->write_data(fp,recvrow,buf);
380 }
381
382 } else {
383 MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
384 MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world);
385 }
386
387 memory->destroy(buf);
388 }
389
390 /* ----------------------------------------------------------------------
391 write out Velocities section of data file
392 ------------------------------------------------------------------------- */
393
velocities()394 void WriteData::velocities()
395 {
396 // communication buffer for all my Atom info
397 // max_size = largest buffer needed by any proc
398
399 int ncol = atom->avec->size_velocity + 1;
400
401 int sendrow = atom->nlocal;
402 int maxrow;
403 MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
404
405 double **buf;
406 if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
407 else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
408
409 // pack my velocity data into buf
410
411 atom->avec->pack_vel(buf,tag_offset);
412
413 // write one chunk of velocities per proc to file
414 // proc 0 pings each proc, receives its chunk, writes to file
415 // all other procs wait for ping, send their chunk to proc 0
416
417 int tmp,recvrow;
418 MPI_Status status;
419 MPI_Request request;
420
421 if (me == 0) {
422 fprintf(fp,"\nVelocities\n\n");
423 for (int iproc = 0; iproc < nprocs; iproc++) {
424 if (iproc) {
425 MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request);
426 MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
427 MPI_Wait(&request,&status);
428 MPI_Get_count(&status,MPI_DOUBLE,&recvrow);
429 recvrow /= ncol;
430 } else recvrow = sendrow;
431
432 atom->avec->write_vel(fp,recvrow,buf);
433 }
434
435 } else {
436 MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
437 MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world);
438 }
439
440 memory->destroy(buf);
441 }
442
443 /* ----------------------------------------------------------------------
444 write out Bonds section of data file
445 ------------------------------------------------------------------------- */
446
bonds()447 void WriteData::bonds()
448 {
449 // communication buffer for all my Bond info
450
451 int ncol = 3;
452 int sendrow = static_cast<int> (nbonds_local);
453 int maxrow;
454 MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
455
456 int **buf;
457 if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
458 else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
459
460 // pack my bond data into buf
461
462 atom->avec->pack_bond(buf);
463
464 // write one chunk of info per proc to file
465 // proc 0 pings each proc, receives its chunk, writes to file
466 // all other procs wait for ping, send their chunk to proc 0
467
468 int tmp,recvrow;
469 MPI_Status status;
470 MPI_Request request;
471
472 int index = 1;
473 if (me == 0) {
474 fprintf(fp,"\nBonds\n\n");
475 for (int iproc = 0; iproc < nprocs; iproc++) {
476 if (iproc) {
477 MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request);
478 MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
479 MPI_Wait(&request,&status);
480 MPI_Get_count(&status,MPI_INT,&recvrow);
481 recvrow /= ncol;
482 } else recvrow = sendrow;
483
484 atom->avec->write_bond(fp,recvrow,buf,index);
485 index += recvrow;
486 }
487
488 } else {
489 MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
490 MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world);
491 }
492
493 memory->destroy(buf);
494 }
495
496 /* ----------------------------------------------------------------------
497 write out Angles section of data file
498 ------------------------------------------------------------------------- */
499
angles()500 void WriteData::angles()
501 {
502 // communication buffer for all my Angle info
503
504 int ncol = 4;
505 int sendrow = static_cast<int> (nangles_local);
506 int maxrow;
507 MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
508
509 int **buf;
510 if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
511 else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
512
513 // pack my angle data into buf
514
515 atom->avec->pack_angle(buf);
516
517 // write one chunk of info per proc to file
518 // proc 0 pings each proc, receives its chunk, writes to file
519 // all other procs wait for ping, send their chunk to proc 0
520
521 int tmp,recvrow;
522 MPI_Status status;
523 MPI_Request request;
524
525 int index = 1;
526 if (me == 0) {
527 fprintf(fp,"\nAngles\n\n");
528 for (int iproc = 0; iproc < nprocs; iproc++) {
529 if (iproc) {
530 MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request);
531 MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
532 MPI_Wait(&request,&status);
533 MPI_Get_count(&status,MPI_INT,&recvrow);
534 recvrow /= ncol;
535 } else recvrow = sendrow;
536
537 atom->avec->write_angle(fp,recvrow,buf,index);
538 index += recvrow;
539 }
540
541 } else {
542 MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
543 MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world);
544 }
545
546 memory->destroy(buf);
547 }
548
549 /* ----------------------------------------------------------------------
550 write out Dihedrals section of data file
551 ------------------------------------------------------------------------- */
552
dihedrals()553 void WriteData::dihedrals()
554 {
555 // communication buffer for all my Dihedral info
556 // max_size = largest buffer needed by any proc
557
558 int ncol = 5;
559
560 int *tag = atom->tag;
561 int *num_dihedral = atom->num_dihedral;
562 int **dihedral_atom2 = atom->dihedral_atom2;
563 int nlocal = atom->nlocal;
564 int newton_bond = force->newton_bond;
565
566 int i,j;
567 int sendrow = 0;
568 if (newton_bond) {
569 for (i = 0; i < nlocal; i++)
570 sendrow += num_dihedral[i];
571 } else {
572 for (i = 0; i < nlocal; i++)
573 for (j = 0; j < num_dihedral[i]; j++)
574 if (tag[i] == dihedral_atom2[i][j]) sendrow++;
575 }
576
577 int maxrow;
578 MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
579
580 int **buf;
581 if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
582 else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
583
584 // pack my dihedral data into buf
585
586 atom->avec->pack_dihedral(buf);
587
588 // write one chunk of info per proc to file
589 // proc 0 pings each proc, receives its chunk, writes to file
590 // all other procs wait for ping, send their chunk to proc 0
591
592 int tmp,recvrow;
593 MPI_Status status;
594 MPI_Request request;
595
596 int index = 1;
597 if (me == 0) {
598 fprintf(fp,"\nDihedrals\n\n");
599 for (int iproc = 0; iproc < nprocs; iproc++) {
600 if (iproc) {
601 MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request);
602 MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
603 MPI_Wait(&request,&status);
604 MPI_Get_count(&status,MPI_INT,&recvrow);
605 recvrow /= ncol;
606 } else recvrow = sendrow;
607
608 atom->avec->write_dihedral(fp,recvrow,buf,index);
609 index += recvrow;
610 }
611
612 } else {
613 MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
614 MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world);
615 }
616
617 memory->destroy(buf);
618 }
619
620 /* ----------------------------------------------------------------------
621 write out Impropers section of data file
622 ------------------------------------------------------------------------- */
623
impropers()624 void WriteData::impropers()
625 {
626 // communication buffer for all my Improper info
627 // max_size = largest buffer needed by any proc
628
629 int ncol = 5;
630
631 int *tag = atom->tag;
632 int *num_improper = atom->num_improper;
633 int **improper_atom2 = atom->improper_atom2;
634 int nlocal = atom->nlocal;
635 int newton_bond = force->newton_bond;
636
637 int i,j;
638 int sendrow = 0;
639 if (newton_bond) {
640 for (i = 0; i < nlocal; i++)
641 sendrow += num_improper[i];
642 } else {
643 for (i = 0; i < nlocal; i++)
644 for (j = 0; j < num_improper[i]; j++)
645 if (tag[i] == improper_atom2[i][j]) sendrow++;
646 }
647
648 int maxrow;
649 MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
650
651 int **buf;
652 if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
653 else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
654
655 // pack my improper data into buf
656
657 atom->avec->pack_improper(buf);
658
659 // write one chunk of info per proc to file
660 // proc 0 pings each proc, receives its chunk, writes to file
661 // all other procs wait for ping, send their chunk to proc 0
662
663 int tmp,recvrow;
664 MPI_Status status;
665 MPI_Request request;
666
667 int index = 1;
668 if (me == 0) {
669 fprintf(fp,"\nImpropers\n\n");
670 for (int iproc = 0; iproc < nprocs; iproc++) {
671 if (iproc) {
672 MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_INT,iproc,0,world,&request);
673 MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
674 MPI_Wait(&request,&status);
675 MPI_Get_count(&status,MPI_INT,&recvrow);
676 recvrow /= ncol;
677 } else recvrow = sendrow;
678
679 atom->avec->write_improper(fp,recvrow,buf,index);
680 index += recvrow;
681 }
682
683 } else {
684 MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
685 MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_INT,0,0,world);
686 }
687
688 memory->destroy(buf);
689 }
690
691 /* ----------------------------------------------------------------------
692 write out Mth section of data file owned by Fix ifix
693 ------------------------------------------------------------------------- */
694
fix(int ifix,int mth)695 void WriteData::fix(int ifix, int mth)
696 {
697 // communication buffer for Fix info
698
699 int sendrow,ncol;
700 modify->fix[ifix]->write_data_section_size(mth,sendrow,ncol);
701 int maxrow;
702 MPI_Allreduce(&sendrow,&maxrow,1,MPI_INT,MPI_MAX,world);
703
704 double **buf;
705 if (me == 0) memory->create(buf,MAX(1,maxrow),ncol,"write_data:buf");
706 else memory->create(buf,MAX(1,sendrow),ncol,"write_data:buf");
707
708 // pack my fix data into buf
709
710 modify->fix[ifix]->write_data_section_pack(mth,buf);
711
712 // write one chunk of info per proc to file
713 // proc 0 pings each proc, receives its chunk, writes to file
714 // all other procs wait for ping, send their chunk to proc 0
715
716 int tmp,recvrow;
717 MPI_Status status;
718 MPI_Request request;
719
720 int index = 1;
721 if (me == 0) {
722 modify->fix[ifix]->write_data_section_keyword(mth,fp);
723 for (int iproc = 0; iproc < nprocs; iproc++) {
724 if (iproc) {
725 MPI_Irecv(&buf[0][0],maxrow*ncol,MPI_DOUBLE,iproc,0,world,&request);
726 MPI_Send(&tmp,0,MPI_INT,iproc,0,world);
727 MPI_Wait(&request,&status);
728 MPI_Get_count(&status,MPI_DOUBLE,&recvrow);
729 recvrow /= ncol;
730 } else recvrow = sendrow;
731
732 modify->fix[ifix]->write_data_section(mth,fp,recvrow,buf,index);
733 index += recvrow;
734 }
735
736 } else {
737 MPI_Recv(&tmp,0,MPI_INT,0,0,world,&status);
738 MPI_Rsend(&buf[0][0],sendrow*ncol,MPI_DOUBLE,0,0,world);
739 }
740
741 memory->destroy(buf);
742 }
743