1 /* ----------------------------------------------------------------------
2 SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer
3 http://sparta.sandia.gov
4 Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov
5 Sandia National Laboratories
6
7 Copyright (2014) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level SPARTA directory.
13 ------------------------------------------------------------------------- */
14
15 #include "mpi.h"
16 #include "spatype.h"
17 #include "string.h"
18 #include "write_isurf.h"
19 #include "update.h"
20 #include "domain.h"
21 #include "grid.h"
22 #include "surf.h"
23 #include "comm.h"
24 #include "modify.h"
25 #include "fix.h"
26 #include "fix_ablate.h"
27 #include "input.h"
28 #include "memory.h"
29 #include "error.h"
30
31 using namespace SPARTA_NS;
32
33 enum{INT,DOUBLE};
34
35 /* ---------------------------------------------------------------------- */
36
WriteISurf(SPARTA * sparta)37 WriteISurf::WriteISurf(SPARTA *sparta) : Pointers(sparta) {}
38
39 /* ---------------------------------------------------------------------- */
40
command(int narg,char ** arg)41 void WriteISurf::command(int narg, char **arg)
42 {
43 MPI_Comm_rank(world,&me);
44 dim = domain->dimension;
45
46 if (!surf->exist || surf->implicit == 0)
47 error->all(FLERR,"Cannot write_isurf when implicit surfs do not exist");
48
49 if (narg < 6) error->all(FLERR,"Illegal write_isurf command");
50
51 ggroup = grid->find_group(arg[0]);
52 if (ggroup < 0) error->all(FLERR,"Write_isurf grid group ID does not exist");
53 groupbit = grid->bitmask[ggroup];
54
55 nx = input->inumeric(FLERR,arg[1]);
56 ny = input->inumeric(FLERR,arg[2]);
57 nz = input->inumeric(FLERR,arg[3]);
58
59 if (dim == 2 && nz != 1) error->all(FLERR,"Invalid write_isurf command");
60
61 // if filename contains a "*", replace with current timestep
62
63 char *ptr;
64 int n = strlen(arg[4]) + 16;
65 char *file = new char[n];
66
67 if ((ptr = strchr(arg[4],'*'))) {
68 *ptr = '\0';
69 sprintf(file,"%s" BIGINT_FORMAT "%s",arg[4],update->ntimestep,ptr+1);
70 } else strcpy(file,arg[4]);
71
72 // ablation fix ID
73
74 char *ablateID = arg[5];
75 int ifix = modify->find_fix(ablateID);
76 if (ifix < 0)
77 error->all(FLERR,"Fix ID for write_isurf does not exist");
78 if (strcmp(modify->fix[ifix]->style,"ablate") != 0)
79 error->all(FLERR,"Fix for write_surf is not a fix ablate");
80 ablate = (FixAblate *) modify->fix[ifix];
81
82 // check that group and Nx,Ny,Nz match FixAblate
83
84 if (ggroup != ablate->igroup)
85 error->all(FLERR,"Write_isurf group does not match fix ablate group");
86
87 if (nx != ablate->nx || ny != ablate->ny || nz != ablate->nz)
88 error->all(FLERR,"Write_isurf Nxyz does not match fix ablate Nxyz");
89
90 // process optional command line args
91
92 precision = DOUBLE;
93
94 int iarg = 6;
95 while (iarg < narg) {
96 if (strcmp(arg[iarg],"precision") == 0) {
97 if (iarg+2 > narg) error->all(FLERR,"Invalid write_isurf command");
98 if (strcmp(arg[iarg+1],"int") == 0) precision = INT;
99 else if (strcmp(arg[iarg+1],"double") == 0) precision = DOUBLE;
100 else error->all(FLERR,"Invalid write_isurf command");
101 iarg += 2;
102 } else error->all(FLERR,"Invalid write_isurf command");
103 }
104
105 // collect all corner point data into one big vector
106
107 if (me == 0 && screen) fprintf(screen,"Writing isurf file ...\n");
108
109 MPI_Barrier(world);
110 double time1 = MPI_Wtime();
111
112 dbuf = dbufall = NULL;
113
114 collect_values();
115
116 MPI_Barrier(world);
117 double time2 = MPI_Wtime();
118
119 // write file
120
121 FILE *fp;
122
123 if (me == 0) {
124 fp = fopen(file,"wb");
125 if (!fp) {
126 char str[128];
127 sprintf(str,"Cannot open grid corner point file %s",file);
128 error->one(FLERR,str);
129 }
130 }
131
132 if (me == 0) write_file(fp);
133
134 // close file
135
136 if (me == 0) fclose(fp);
137
138 memory->destroy(dbuf);
139 memory->destroy(dbufall);
140
141 MPI_Barrier(world);
142 double time3 = MPI_Wtime();
143
144 // stats
145
146 double time_total = time3-time1;
147
148 if (comm->me == 0) {
149 if (screen) {
150 fprintf(screen," corner points = %d\n",ncorner);
151 fprintf(screen," CPU time = %g secs\n",time_total);
152 fprintf(screen," collect/write percent = %g %g\n",
153 100.0*(time2-time1)/time_total,100.0*(time3-time2)/time_total);
154 }
155
156 if (logfile) {
157 fprintf(logfile," corner points = %d\n",ncorner);
158 fprintf(logfile," CPU time = %g secs\n",time_total);
159 fprintf(logfile," collect/write percent = %g %g\n",
160 100.0*(time2-time1)/time_total,100.0*(time3-time2)/time_total);
161 }
162 }
163 }
164
165 /* ----------------------------------------------------------------------
166 collect corner point values into one big array
167 each proc owns copy, use MPI_Allreduce to sum
168 ------------------------------------------------------------------------- */
169
collect_values()170 void WriteISurf::collect_values()
171 {
172 double *cornerlo = ablate->cornerlo;
173 double *xyzsize = ablate->xyzsize;
174
175 bigint bncorner = bigint (nx+1) * (ny+1);
176 if (dim == 3) bncorner *= (nz+1);
177 if (bncorner > MAXSMALLINT) error->all(FLERR,"Write_isurf grid is too large");
178 ncorner = bncorner;
179
180 memory->create(dbuf,ncorner,"write_isurf:dbuf");
181 memory->create(dbufall,ncorner,"write_isurf:dbufall");
182 for (int i = 0; i < ncorner; i++) dbuf[i] = 0.0;
183
184 Grid::ChildCell *cells = grid->cells;
185 Grid::ChildInfo *cinfo = grid->cinfo;
186 int nglocal = grid->nlocal;
187
188 int ix,iy,iz,index;
189 double **array_grid = ablate->array_grid;
190
191 for (int icell = 0; icell < nglocal; icell++) {
192 if (!(cinfo[icell].mask & groupbit)) continue;
193 if (cells[icell].nsplit <= 0) continue;
194
195 // ix,iy,iz = cell index (1 to Nxyz) within array of grid cells
196
197 ix =
198 static_cast<int> ((cells[icell].lo[0]-cornerlo[0]) / xyzsize[0] + 0.5) + 1;
199 iy =
200 static_cast<int> ((cells[icell].lo[1]-cornerlo[1]) / xyzsize[1] + 0.5) + 1;
201 iz =
202 static_cast<int> ((cells[icell].lo[2]-cornerlo[2]) / xyzsize[2] + 0.5) + 1;
203
204 // index = corner point index, 0 to (Nx+1)*(Ny+1)*(Nz+1) - 1
205 // x varies fastest, then y, z slowest
206 // this is for lower/left/bottom corner point of icell's 4/8 points
207
208 index = (iz-1)*(ny+1)*(nx+1) + (iy-1)*(nx+1) + ix-1;
209
210 // always copy lower/left/bottom corner point
211 // copy other corner points if on upper boundary of Nx,Ny,Nz
212
213 dbuf[index] = array_grid[icell][0];
214 if (ix == nx) dbuf[index+1] = array_grid[icell][1];
215 if (iy == ny) dbuf[index+nx+1] = array_grid[icell][2];
216 if (ix == nx && iy == ny) dbuf[index+nx+2] = array_grid[icell][3];
217
218 if (iz == nz && dim == 3) {
219 index += (ny+1)*(nx+1);
220 dbuf[index] = array_grid[icell][4];
221 if (ix == nx) dbuf[index+1] = array_grid[icell][5];
222 if (iy == ny) dbuf[index+nx+1] = array_grid[icell][6];
223 if (ix == nx && iy == ny) dbuf[index+nx+2] = array_grid[icell][7];
224 }
225 }
226
227 // MPI_Allreduce to sum dbuf across all procs
228 // so that proc 0 has a copy to write to file
229
230 MPI_Allreduce(dbuf,dbufall,ncorner,MPI_DOUBLE,MPI_SUM,world);
231 }
232
233 /* ----------------------------------------------------------------------
234 write grid corner point file
235 only called by proc 0
236 ------------------------------------------------------------------------- */
237
write_file(FILE * fp)238 void WriteISurf::write_file(FILE *fp)
239 {
240 int nxyz[3];
241 nxyz[0] = nx + 1;
242 nxyz[1] = ny + 1;
243 nxyz[2] = nz + 1;
244 fwrite(nxyz,sizeof(int),dim,fp);
245
246 if (precision == INT) {
247 uint8_t *ibuf8;
248 memory->create(ibuf8,ncorner,"write_isurf:ibuf8");
249 for (int i = 0; i < ncorner; i++)
250 ibuf8[i] = static_cast<int> (dbufall[i]);
251 fwrite(ibuf8,sizeof(uint8_t),ncorner,fp);
252 memory->destroy(ibuf8);
253 }
254
255 if (precision == DOUBLE) fwrite(dbufall,sizeof(double),ncorner,fp);
256 }
257