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