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 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 } 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 } 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 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 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