[694] | 1 | #include "mpi.hpp" |
---|
[688] | 2 | #include <netcdf.h> |
---|
| 3 | #include "errhandle.hpp" |
---|
| 4 | #include "libmapper.hpp" |
---|
| 5 | #include "mapper.hpp" |
---|
| 6 | |
---|
| 7 | #include "node.hpp" |
---|
| 8 | |
---|
| 9 | using namespace sphereRemap ; |
---|
| 10 | |
---|
| 11 | int ncread(double *lon, double *lat, double* val, int nElt, const char* filename) |
---|
| 12 | { |
---|
| 13 | int ncid, blonid, blatid, valid, neltid; |
---|
| 14 | size_t nEltCheck; |
---|
| 15 | exit_on_failure(nc_open(filename, NC_NOWRITE, &ncid), std::string("Cannot open netCDF file ") + filename); |
---|
| 16 | exit_on_failure(nc_inq_dimid(ncid, "elt", &neltid), std::string("No dimension elt in file ") + filename); |
---|
| 17 | nc_inq_dimlen(ncid, neltid, &nEltCheck); |
---|
| 18 | exit_on_failure(nElt != nEltCheck, std::string("Array sizes do not match!")); |
---|
| 19 | |
---|
| 20 | nc_inq_varid(ncid, "bounds_lon", &blonid); |
---|
| 21 | nc_inq_varid(ncid, "bounds_lat", &blatid); |
---|
| 22 | nc_inq_varid(ncid, "val", &valid); |
---|
| 23 | nc_get_var_double(ncid, blonid, lon); |
---|
| 24 | nc_get_var_double(ncid, blatid, lat); |
---|
| 25 | nc_get_var_double(ncid, valid, val); |
---|
| 26 | nc_close(ncid); |
---|
| 27 | } |
---|
| 28 | |
---|
| 29 | int ncwriteValue(double *val, const char* filename) |
---|
| 30 | { |
---|
| 31 | int ncid, valid ; |
---|
| 32 | nc_open(filename, NC_WRITE, &ncid); |
---|
| 33 | nc_inq_varid(ncid, "val", &valid); |
---|
| 34 | nc_put_var_double(ncid, valid, val); |
---|
| 35 | nc_close(ncid); |
---|
| 36 | } |
---|
| 37 | |
---|
| 38 | |
---|
| 39 | |
---|
| 40 | void compute_distribution(int nGlobalElts, int &start, int &nLocalElts) |
---|
| 41 | { |
---|
| 42 | int mpiSize, mpiRank; |
---|
| 43 | MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); |
---|
| 44 | MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); |
---|
| 45 | start = 0; |
---|
| 46 | nLocalElts = 0; |
---|
| 47 | for (int i = 0; i <= mpiRank; i++) |
---|
| 48 | { |
---|
| 49 | start += nLocalElts; |
---|
| 50 | nLocalElts = nGlobalElts/mpiSize; |
---|
| 51 | if (i < nGlobalElts % mpiSize) nLocalElts++; |
---|
| 52 | } |
---|
| 53 | } |
---|
| 54 | |
---|
| 55 | int main() |
---|
| 56 | { |
---|
| 57 | int interpOrder = 2; |
---|
| 58 | /* low resolution */ |
---|
| 59 | /* |
---|
| 60 | char srcFile[] = "h14.nc"; |
---|
| 61 | char dstFile[] = "r180x90.nc"; |
---|
| 62 | double srcPole[] = {0, 0, 0}; |
---|
| 63 | double dstPole[] = {0, 0, 1}; |
---|
| 64 | int nSrcElt = 13661; |
---|
| 65 | int nDstElt = 16200; |
---|
| 66 | */ |
---|
| 67 | /* high resolution */ |
---|
| 68 | |
---|
| 69 | char srcFile[] = "t740.nc"; |
---|
| 70 | char dstFile[] = "r1440x720.nc"; |
---|
| 71 | double srcPole[] = {0, 0, 0}; |
---|
| 72 | double dstPole[] = {0, 0, 1}; |
---|
| 73 | int nSrcElt = 741034; |
---|
| 74 | int nDstElt = 1036800; |
---|
| 75 | |
---|
| 76 | int nVert = 10; |
---|
| 77 | int nSrcLocal, nDstLocal, startSrc, startDst; |
---|
| 78 | int nWeight; |
---|
| 79 | int mpi_rank, mpi_size ; |
---|
| 80 | |
---|
| 81 | MPI_Init(NULL, NULL); |
---|
| 82 | MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank) ; |
---|
| 83 | MPI_Comm_size(MPI_COMM_WORLD, &mpi_size) ; |
---|
| 84 | |
---|
| 85 | double *srcLon = (double *) malloc(nSrcElt*nVert*sizeof(double)); |
---|
| 86 | double *srcLat = (double *) malloc(nSrcElt*nVert*sizeof(double)); |
---|
| 87 | double *srcVal = (double *) malloc(nSrcElt*sizeof(double)); |
---|
| 88 | double *dstLon = (double *) malloc(nDstElt*nVert*sizeof(double)); |
---|
| 89 | double *dstLat = (double *) malloc(nDstElt*nVert*sizeof(double)); |
---|
| 90 | double *dstVal = (double *) malloc(nDstElt*sizeof(double)); |
---|
| 91 | double *globalVal = (double *) malloc(nDstElt*sizeof(double)); |
---|
| 92 | |
---|
| 93 | ncread(srcLon, srcLat, srcVal, nSrcElt, srcFile); |
---|
| 94 | ncread(dstLon, dstLat, dstVal, nDstElt, dstFile); |
---|
| 95 | |
---|
| 96 | compute_distribution(nSrcElt, startSrc, nSrcLocal); |
---|
| 97 | compute_distribution(nDstElt, startDst, nDstLocal); |
---|
| 98 | |
---|
| 99 | |
---|
| 100 | |
---|
| 101 | Mapper mapper(MPI_COMM_WORLD); |
---|
| 102 | mapper.setVerbosity(PROGRESS) ; |
---|
| 103 | mapper.setSourceMesh(srcLon + startSrc*nVert, srcLat + startSrc*nVert, nVert, nSrcLocal, srcPole ) ; |
---|
| 104 | mapper.setTargetMesh(dstLon + startDst*nVert, dstLat + startDst*nVert, nVert, nDstLocal, dstPole ) ; |
---|
| 105 | |
---|
| 106 | for(int i=0;i<nDstElt;i++) dstVal[i]=0 ; |
---|
| 107 | mapper.setSourceValue(srcVal+ startSrc) ; |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | vector<double> timings = mapper.computeWeights(interpOrder); |
---|
| 111 | |
---|
| 112 | |
---|
| 113 | for(int i=0;i<mapper.nWeights;i++) dstVal[mapper.dstAddress[i]]=dstVal[mapper.dstAddress[i]]+mapper.remapMatrix[i]*srcVal[mapper.sourceWeightId[i]]; |
---|
| 114 | |
---|
| 115 | /* |
---|
| 116 | remap_get_num_weights(srcLon + startSrc*nVert, srcLat + startSrc*nVert, nVert, nSrcLocal, srcPole, |
---|
| 117 | dstLon + startDst*nVert, dstLat + startDst*nVert, nVert, nDstLocal, dstPole, |
---|
| 118 | interpOrder, &nWeight); |
---|
| 119 | |
---|
| 120 | double *weights = (double *) malloc(nWeight*sizeof(double)); |
---|
| 121 | int *srcIdx = (int *) malloc(nWeight*sizeof(int)); |
---|
| 122 | int *srcRank = (int *) malloc(nWeight*sizeof(int)); |
---|
| 123 | int *dstIdx = (int *) malloc(nWeight*sizeof(int)); |
---|
| 124 | |
---|
| 125 | remap_get_weights(weights, srcIdx, dstIdx); |
---|
| 126 | |
---|
| 127 | free(srcLon); free(srcLat); |
---|
| 128 | free(dstLon); free(dstLat); |
---|
| 129 | free(srcIdx); free(dstIdx); |
---|
| 130 | free(srcRank); |
---|
| 131 | free(weights); |
---|
| 132 | |
---|
| 133 | #ifdef DEBUG |
---|
| 134 | memory_report(); |
---|
| 135 | #endif |
---|
| 136 | */ |
---|
| 137 | int* displs=new int[mpi_size] ; |
---|
| 138 | int* recvCount=new int[mpi_size] ; |
---|
| 139 | |
---|
| 140 | MPI_Gather(&startDst,1,MPI_INT,displs,1,MPI_INT,0,MPI_COMM_WORLD) ; |
---|
| 141 | MPI_Gather(&nDstLocal,1,MPI_INT,recvCount,1,MPI_INT,0,MPI_COMM_WORLD) ; |
---|
| 142 | MPI_Gatherv(dstVal,nDstLocal,MPI_DOUBLE,globalVal,recvCount,displs,MPI_DOUBLE,0,MPI_COMM_WORLD) ; |
---|
| 143 | if (mpi_rank==0) ncwriteValue(globalVal, dstFile); |
---|
| 144 | |
---|
| 145 | MPI_Finalize(); |
---|
| 146 | return 0; |
---|
| 147 | } |
---|