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 "phastaIO.h" 19 20 //Provided by phastaIO 21 void Gather_Headers( int* fileDescriptor, std::vector< std::string >& headers ); 22 23 char read_solution(double** solutiono, int* size, int* nshgo, int* ndofo, 24 int nump, int rank, int timestep, char* casedir); 25 26 std::set<int>* find_timesteps(char* casedir, int nump); 27 double compare_solution(char* lpath, char* rpath, int timestep, int nump); 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 assert(argc>2); 39 char* lpath = argv[1]; 40 char* rpath = argv[2]; 41 42 int ndof; 43 int nshg; 44 int solsize; 45 double* solution; 46 47 std::set<int>* l_timesteps = find_timesteps(lpath, size); 48 std::set<int>* r_timesteps = find_timesteps(rpath, size); 49 std::set<int>* timesteps_to_check = new std::set<int>; 50 std::set_intersection(l_timesteps->begin(), l_timesteps->end(), 51 r_timesteps->begin(), r_timesteps->end(), 52 std::inserter(*timesteps_to_check, timesteps_to_check->begin())); 53 delete l_timesteps; 54 delete r_timesteps; 55 if(rank == 0) 56 printf("Found %d common timesteps\n", 57 timesteps_to_check->size()); 58 #ifdef DBGONLY 59 read_solution(&solution, &solsize, &nshg, &ndof, size, rank, 0, "./"); 60 printf("nshg: %d, ndof: %d\n", nshg, ndof); 61 assert(solsize == ndof*nshg); 62 #endif 63 double maxerror = 0.0; 64 double error; 65 double gblmaxerror; 66 for(std::set<int>::iterator i = timesteps_to_check->begin(); 67 i!=timesteps_to_check->end();i++) 68 { 69 error = compare_solution(lpath, rpath, *i, size); 70 if(error>maxerror) maxerror = error; 71 } 72 delete timesteps_to_check; 73 MPI_Barrier(MPI_COMM_WORLD); 74 MPI_Reduce(&maxerror, &gblmaxerror, 1, MPI_DOUBLE, MPI_MAX, 0, 75 MPI_COMM_WORLD); 76 if(rank == 0) printf("Maximum difference across all timesteps: %e\n", 77 gblmaxerror); 78 MPI_Finalize(); 79 } 80 double compare_solution(char* lpath, char* rpath, int timestep, int nump) 81 { 82 int rank; 83 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 84 double* lsol; 85 double* rsol; 86 int lsize; 87 int rsize; 88 89 read_solution(&lsol, &lsize, NULL, NULL, nump, rank, timestep, lpath); 90 read_solution(&rsol, &rsize, NULL, NULL, nump, rank, timestep, rpath); 91 92 double maxdiff=0.0; 93 double gblmaxdiff; 94 if(lsize != rsize) 95 { 96 printf("Error: Solution sizes different: %d, %d\n", 97 lsize, rsize); 98 assert(lsize == rsize); 99 } 100 for(int i=0;i<lsize;i++) 101 { 102 double diff = fabs(lsol[i]-rsol[i]); 103 if(diff > maxdiff) maxdiff = diff; 104 } 105 free(lsol); 106 free(rsol); 107 MPI_Reduce(&maxdiff, &gblmaxdiff, 1, MPI_DOUBLE, MPI_MAX, 0, 108 MPI_COMM_WORLD); 109 MPI_Barrier(MPI_COMM_WORLD); //TODO: debugging only 110 if(rank == 0) 111 printf("Timestep: %d, maximum difference: %e\n", timestep, gblmaxdiff); 112 return gblmaxdiff; 113 114 } 115 char read_solution(double** solutiono, int* size, int* nshgo, int* ndofo, 116 int nump, int rank, int timestep, char* casedir) 117 { 118 int iarray[10]; 119 const char* iformat = "binary"; 120 int ithree=3; 121 int igeombc; 122 char* fn; 123 int nshg; 124 int ndof; 125 double* solution; 126 if(nump == 1) 127 asprintf(&fn, "%s/restart.%d.%d", 128 casedir,timestep,rank+1); 129 else 130 asprintf(&fn, "%s/%d-procs_case/restart.%d.%d", 131 casedir, nump,timestep,rank+1); 132 openfile(fn, "read", &igeombc); 133 //TODO: error handle 134 readheader(&igeombc, "solution", (void*) iarray, &ithree, "integer", iformat); 135 nshg = iarray[0]; 136 ndof = iarray[1]; 137 if(size != NULL) 138 *size = nshg*ndof; 139 solution = (double*) malloc(sizeof(double)*nshg*ndof); 140 readdatablock(&igeombc, "solution", solution, size, "double", iformat); 141 closefile(&igeombc, "read"); 142 if(solutiono != NULL) 143 *solutiono = solution; 144 if(nshgo != NULL) 145 *nshgo = nshg; 146 if(ndofo != NULL) 147 *ndofo = ndof; 148 free(fn); 149 return(0); 150 } 151 152 std::set<int>* find_timesteps(char* casedir, int nump) 153 { 154 char* path; 155 char* fullpath; 156 DIR* dir; 157 struct dirent* d; 158 int part, ts; 159 std::set<int>* step_list = new std::set<int>; 160 161 if(nump == 1) 162 asprintf(&path, "%s", casedir); 163 else 164 asprintf(&path, "%s/%d-procs_case", casedir, nump); 165 dir = opendir(path); 166 if(!dir) 167 { 168 perror("Error opening case: "); 169 MPI_Abort(MPI_COMM_WORLD,1); 170 } 171 while((d=readdir(dir))) 172 { 173 asprintf(&fullpath, "%s/%s", path, d->d_name); 174 if(sscanf(d->d_name, "restart.%d.%d", &ts, &part)==2) 175 { 176 step_list->insert(ts); 177 } 178 } 179 return(step_list); 180 free(path); 181 } 182 183 void check_ilwork(int* ilwork, int nlwork) 184 { 185 int numtask = ilwork[0]; 186 int itkbeg = 0; //task offset 187 printf("%d tasks\n", numtask); 188 189 for(int i=0;i<numtask;i++) 190 { 191 int itag = ilwork[itkbeg+1]; //mpi tag 192 int iacc = ilwork[itkbeg+2]; //0 for slave, 1 for master 193 int iother = ilwork[itkbeg+3]-1; //other rank (see ctypes.f for off by one) 194 int numseg = ilwork[itkbeg+4]; //number of segments 195 printf("Comm with rank: %d\n", iother); 196 for(int j=0;j<numseg;j++) 197 { 198 int isgbeg = ilwork[itkbeg+5+(j*2)]; //first idx of seg 199 int lenseg = ilwork[itkbeg+6+(j*2)]; //length of seg 200 201 printf("isgbeg: %d, len: %d\n", isgbeg, lenseg); 202 assert(itkbeg+6+(j*2) < nlwork); 203 } 204 itkbeg+= 4+2*numseg; 205 } 206 } 207 208 void read_ilwork() 209 { 210 int rank; 211 int size; 212 213 MPI_Comm_rank(MPI_COMM_WORLD, &rank); 214 MPI_Comm_size(MPI_COMM_WORLD, &size); 215 216 int iarray[10]; 217 const char* iformat = "binary"; 218 int ione=1; 219 int itwo=2; 220 int igeombc; 221 int nilwork=0; 222 int tmp; 223 int* ilwork; 224 char* fn; 225 asprintf(&fn, "%d-procs_case/geombc.dat.%d",size,rank+1); 226 openfile(fn, "read", &igeombc); 227 228 std::vector<std::string> headers; 229 Gather_Headers(&igeombc, headers); 230 231 readheader(&igeombc, "size of ilwork array", (void*) &nilwork, &ione, 232 "integer", iformat); 233 assert(nilwork > 0); 234 readheader(&igeombc, "ilwork", (void*) &tmp, &ione, 235 "integer", iformat); 236 assert(tmp == nilwork); 237 238 ilwork = (int*) malloc(sizeof(int)*nilwork); 239 240 readdatablock(&igeombc, "ilwork", ilwork, &nilwork, "integer", iformat); 241 closefile(&igeombc, "read"); 242 } 243