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