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