1 #ifndef __xlC__
2 #define _GNU_SOURCE
3 #endif
4 #define _WITH_DPRINTF
5 #include <cstdio>
6 #include <cstdlib>
7 #include <cassert>
8 #include <cstring>
9 #include <cmath>
10 #include <string>
11 #include <vector>
12 #include <algorithm>
13 #include <set>
14
15 #include <sys/types.h>
16 #include <sys/stat.h>
17 #include <dirent.h>
18
19 #define OMPI_SKIP_MPICXX 1
20 #include <mpi.h>
21
22 #include "syncio.h"
23 #include "posixio.h"
24 #include "phIO.h"
25
26 char read_solution(double** solutiono, int* size, int* nshgo, int* ndofo,
27 int nump, int rank, int timestep, int nSyncFiles, char* casedir);
28 std::set<int>* find_timesteps(char* casedir, int nSyncFiles);
29 double compare_solution(char* lpath, char* rpath,
30 int timestep, int nump, int nSyncFiles);
31 char* getRestartName(int nSyncFiles);
32
main(int argc,char ** argv)33 int main(int argc, char** argv)
34 {
35 int rank;
36 int size;
37 MPI_Init(&argc, &argv);
38
39 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
40 MPI_Comm_size(MPI_COMM_WORLD, &size);
41
42 if(argc != 5) {
43 fprintf(stderr, "argc %d\n", argc);
44 fprintf(stderr,
45 "Usage: %s <left> <right> <numSyncFiles> <tolerance>\n"
46 "where <left> and <right> are different"
47 "N-procs_case directories\n", argv[0]);
48 MPI_Finalize();
49 return 1;
50 }
51 char* lpath = argv[1];
52 char* rpath = argv[2];
53 int nSyncFiles = atoi(argv[3]);
54 double tolerance = atof(argv[4]);
55
56 int ndof;
57 int nshg;
58 int solsize;
59 double* solution;
60
61 std::set<int>* l_timesteps = find_timesteps(lpath, nSyncFiles);
62 std::set<int>* r_timesteps = find_timesteps(rpath, nSyncFiles);
63 std::set<int>* timesteps_to_check = new std::set<int>;
64 std::set_intersection(l_timesteps->begin(), l_timesteps->end(),
65 r_timesteps->begin(), r_timesteps->end(),
66 std::inserter(*timesteps_to_check, timesteps_to_check->begin()));
67 delete l_timesteps;
68 delete r_timesteps;
69 if(rank == 0)
70 printf("Found %d common timesteps\n",
71 timesteps_to_check->size());
72 #ifdef DBGONLY
73 read_solution(&solution, &solsize, &nshg, &ndof,
74 size, rank, 0, numSyncFiles, "./");
75 printf("nshg: %d, ndof: %d\n", nshg, ndof);
76 assert(solsize == ndof*nshg);
77 #endif
78 double maxerror = 0.0;
79 double error;
80 double gblmaxerror;
81 for(std::set<int>::iterator i = timesteps_to_check->begin();
82 i!=timesteps_to_check->end();i++)
83 {
84 error = compare_solution(lpath, rpath, *i, size, nSyncFiles);
85 if(error>maxerror) maxerror = error;
86 }
87 delete timesteps_to_check;
88 MPI_Barrier(MPI_COMM_WORLD);
89 MPI_Reduce(&maxerror, &gblmaxerror, 1, MPI_DOUBLE, MPI_MAX, 0,
90 MPI_COMM_WORLD);
91 if(rank == 0) printf("Maximum difference across all timesteps: %e\n",
92 gblmaxerror);
93 MPI_Finalize();
94 return (gblmaxerror > tolerance);
95 }
compare_solution(char * lpath,char * rpath,int timestep,int nump,int nSyncFiles)96 double compare_solution(char* lpath, char* rpath, int timestep, int nump, int nSyncFiles)
97 {
98 int rank;
99 MPI_Comm_rank(MPI_COMM_WORLD, &rank);
100 double* lsol;
101 double* rsol;
102 int lsize;
103 int rsize;
104
105 read_solution(&lsol, &lsize, NULL, NULL,
106 nump, rank, timestep, nSyncFiles, lpath);
107 read_solution(&rsol, &rsize, NULL, NULL,
108 nump, rank, timestep, nSyncFiles, rpath);
109
110 double maxdiff=0.0;
111 double gblmaxdiff;
112 if(lsize != rsize)
113 {
114 printf("Error: Solution sizes different: %d, %d\n",
115 lsize, rsize);
116 assert(lsize == rsize);
117 }
118 for(int i=0;i<lsize;i++)
119 {
120 double diff = fabs(lsol[i]-rsol[i]);
121 if(diff > maxdiff) maxdiff = diff;
122 }
123 free(lsol);
124 free(rsol);
125 MPI_Reduce(&maxdiff, &gblmaxdiff, 1, MPI_DOUBLE, MPI_MAX, 0,
126 MPI_COMM_WORLD);
127 MPI_Barrier(MPI_COMM_WORLD); //TODO: debugging only
128 if(rank == 0)
129 printf("Timestep: %d, maximum difference: %e\n", timestep, gblmaxdiff);
130 return gblmaxdiff;
131
132 }
read_solution(double ** solutiono,int * size,int * nshgo,int * ndofo,int nump,int rank,int timestep,int nSyncFiles,char * casedir)133 char read_solution(double** solutiono, int* size, int* nshgo, int* ndofo,
134 int nump, int rank, int timestep, int nSyncFiles, char* casedir)
135 {
136 int iarray[10];
137 const char* iformat = "binary";
138 int ithree=3;
139 int igeombc;
140 char* fn;
141 int nshg;
142 int ndof;
143 double* solution;
144 phio_fp fp;
145 if( nSyncFiles == 0 )
146 posixio_setup(&fp, 'r');
147 else if( nSyncFiles > 0 )
148 syncio_setup_read(nSyncFiles, &fp);
149 char rname[1024];
150 phio_constructName(fp,"restart",rname);
151 asprintf(&fn,"%s/%s%d.",casedir,rname,timestep);
152 phio_openfile(fn, fp);
153
154 phio_readheader(fp, "solution", (void*) iarray, &ithree, "integer", iformat);
155 nshg = iarray[0];
156 ndof = iarray[1];
157 if(size != NULL)
158 *size = nshg*ndof;
159 solution = (double*) malloc(sizeof(double)*nshg*ndof);
160 phio_readdatablock(fp, "solution", solution, size, "double", iformat);
161 phio_closefile(fp);
162 if(solutiono != NULL)
163 *solutiono = solution;
164 if(nshgo != NULL)
165 *nshgo = nshg;
166 if(ndofo != NULL)
167 *ndofo = ndof;
168 free(fn);
169 return(0);
170 }
171
find_timesteps(char * casedir,int nSyncFiles)172 std::set<int>* find_timesteps(char* casedir, int nSyncFiles)
173 {
174 char* path;
175 DIR* dir;
176 struct dirent* d;
177 int part, ts;
178 std::set<int>* step_list = new std::set<int>;
179
180 asprintf(&path, "%s", casedir);
181 dir = opendir(path);
182 if(!dir)
183 {
184 perror("Error opening case: ");
185 MPI_Abort(MPI_COMM_WORLD,1);
186 }
187 char* rname = getRestartName(nSyncFiles);
188 char* fmt;
189 asprintf(&fmt, "%s.%%d.%%d", rname);
190 while((d=readdir(dir)))
191 {
192 if(sscanf(d->d_name, fmt, &ts, &part)==2)
193 {
194 step_list->insert(ts);
195 }
196 }
197 free(rname);
198 free(fmt);
199 free(path);
200 closedir(dir);
201 return(step_list);
202 }
203
getRestartName(int nSyncFiles)204 char* getRestartName(int nSyncFiles) {
205 char* f;
206 if(0 == nSyncFiles)
207 asprintf(&f, "restart");
208 else if(nSyncFiles > 0)
209 asprintf(&f, "restart-dat");
210 else {
211 fprintf(stderr,
212 "ERROR: the number of sync-io files must be"
213 "greater than or equal to zero\n");
214 MPI_Abort(MPI_COMM_WORLD, 1);
215 return NULL;
216 }
217 return f;
218 }
219