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