[694] | 1 | #include "mpi.hpp" |
---|
[688] | 2 | #include <vector> |
---|
| 3 | #include <cassert> |
---|
| 4 | #include <cstring> |
---|
| 5 | #ifdef WRITE_TIMING |
---|
| 6 | #include <fstream> |
---|
| 7 | #endif |
---|
| 8 | |
---|
| 9 | #include "node.hpp" |
---|
| 10 | #include "grid.hpp" |
---|
| 11 | #include "circle.hpp" // cptRadius |
---|
| 12 | #include "elt.hpp" |
---|
| 13 | #include "meshutil.hpp" // cptArea |
---|
| 14 | #include "mapper.hpp" |
---|
| 15 | #include "cputime.hpp" // cputime |
---|
[1155] | 16 | #include "gridRemap.hpp" |
---|
[688] | 17 | |
---|
| 18 | using namespace sphereRemap ; |
---|
| 19 | |
---|
[1155] | 20 | extern CRemapGrid srcGrid; |
---|
| 21 | #pragma omp threadprivate(srcGrid) |
---|
| 22 | |
---|
| 23 | extern CRemapGrid tgtGrid; |
---|
| 24 | #pragma omp threadprivate(tgtGrid) |
---|
| 25 | |
---|
| 26 | |
---|
[688] | 27 | /* mapper is a ponter to a global class instance whoes members are allocated in the first step (finding the sizes of the weight arrays) |
---|
| 28 | and deallocated during the second step (computing the weights) */ |
---|
| 29 | Mapper *mapper; |
---|
[1155] | 30 | #pragma omp threadprivate(mapper) |
---|
[688] | 31 | |
---|
| 32 | /** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx |
---|
| 33 | This function computes the weights and returns number of weights is returned through the last argument. |
---|
| 34 | To get the actuall weights, use `remap_get_weights`. |
---|
| 35 | */ |
---|
| 36 | extern "C" void remap_get_num_weights(const double* src_bounds_lon, const double* src_bounds_lat, int n_vert_per_cell_src, int n_cell_src, |
---|
| 37 | const double* src_pole, |
---|
| 38 | const double* dst_bounds_lon, const double* dst_bounds_lat, int n_vert_per_cell_dst, int n_cell_dst, |
---|
| 39 | const double* dst_pole, |
---|
| 40 | int order, int* n_weights) |
---|
| 41 | { |
---|
| 42 | assert(src_bounds_lon); |
---|
| 43 | assert(src_bounds_lat); |
---|
| 44 | assert(n_vert_per_cell_src >= 3); |
---|
| 45 | assert(n_cell_src >= 4); |
---|
| 46 | assert(dst_bounds_lon); |
---|
| 47 | assert(dst_bounds_lat); |
---|
| 48 | assert(n_vert_per_cell_dst >= 3); |
---|
| 49 | assert(n_cell_dst >= 4); |
---|
| 50 | assert(1 <= order && order <= 2); |
---|
| 51 | |
---|
| 52 | mapper = new Mapper(MPI_COMM_WORLD); |
---|
| 53 | mapper->setVerbosity(PROGRESS) ; |
---|
| 54 | mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, n_vert_per_cell_src, n_cell_src, src_pole ) ; |
---|
| 55 | mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, n_vert_per_cell_dst, n_cell_dst, dst_pole ) ; |
---|
| 56 | |
---|
| 57 | /* |
---|
| 58 | srcGrid.pole = Coord(src_pole[0], src_pole[1], src_pole[2]); |
---|
| 59 | tgtGrid.pole = Coord(dst_pole[0], dst_pole[1], dst_pole[2]); |
---|
| 60 | |
---|
| 61 | int mpiRank, mpiSize; |
---|
| 62 | MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank); |
---|
| 63 | MPI_Comm_size(MPI_COMM_WORLD, &mpiSize); |
---|
| 64 | |
---|
| 65 | std::vector<Elt> src_elt; src_elt.reserve(n_cell_src); |
---|
| 66 | std::vector<Node> src_msh; src_msh.reserve(n_cell_src); |
---|
| 67 | for (int i = 0; i < n_cell_src; i++) |
---|
| 68 | { |
---|
| 69 | int offs = i*n_vert_per_cell_src; |
---|
| 70 | Elt elt(src_bounds_lon + offs, src_bounds_lat + offs, n_vert_per_cell_src); |
---|
| 71 | elt.src_id.rank = mpiRank; |
---|
| 72 | elt.src_id.ind = i; |
---|
| 73 | src_elt.push_back(elt); |
---|
| 74 | src_msh.push_back(Node(elt.x, cptRadius(elt), &src_elt.back())); |
---|
| 75 | cptEltGeom(src_elt[i], Coord(src_pole[0], src_pole[1], src_pole[2])); |
---|
| 76 | } |
---|
| 77 | std::vector<Elt> dst_elt; dst_elt.reserve(n_cell_dst); |
---|
| 78 | std::vector<Node> dst_msh; dst_msh.reserve(n_cell_dst); |
---|
| 79 | for (int i = 0; i < n_cell_dst; i++) |
---|
| 80 | { |
---|
| 81 | int offs = i*n_vert_per_cell_dst; |
---|
| 82 | Elt elt(dst_bounds_lon + offs, dst_bounds_lat + offs, n_vert_per_cell_dst); |
---|
| 83 | dst_elt.push_back(elt); |
---|
| 84 | dst_msh.push_back(Node(elt.x, cptRadius(elt), &dst_elt.back())); |
---|
| 85 | cptEltGeom(dst_elt[i], Coord(dst_pole[0], dst_pole[1], dst_pole[2])); |
---|
| 86 | } |
---|
| 87 | double tic = cputime(); |
---|
| 88 | mapper = new Mapper(MPI_COMM_WORLD); |
---|
| 89 | mapper->setVerbosity(PROGRESS) ; |
---|
| 90 | mapper->buildSSTree(src_msh, dst_msh); |
---|
| 91 | double tac = cputime(); |
---|
| 92 | vector<double> timings = mapper->computeWeights(dst_elt, src_elt, order); |
---|
| 93 | */ |
---|
| 94 | |
---|
| 95 | vector<double> timings = mapper->computeWeights(order); |
---|
| 96 | /* |
---|
| 97 | #ifdef WRITE_TIMING |
---|
| 98 | int nGloWeights, gloSrcSize, gloDstSize; |
---|
| 99 | int locSrcSize = src_elt.size(); |
---|
| 100 | int locDstSize = dst_elt.size(); |
---|
| 101 | MPI_Reduce(&(mapper->nWeights), &nGloWeights, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); |
---|
| 102 | MPI_Reduce(&locSrcSize, &gloSrcSize, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); |
---|
| 103 | MPI_Reduce(&locDstSize, &gloDstSize, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD); |
---|
| 104 | |
---|
| 105 | if (mpiRank == 0) |
---|
| 106 | { |
---|
| 107 | ofstream timeout; |
---|
| 108 | timeout.open("timings.txt", fstream::out | fstream::app); |
---|
| 109 | timeout << mpiSize << " " << src_elt.size() << " " << dst_elt.size() << " " << nGloWeights |
---|
| 110 | << " " << tac-tic << " " << timings[0] << " " << timings[1] << " " << timings[2] << endl; |
---|
| 111 | timeout.close(); |
---|
| 112 | } |
---|
| 113 | #endif |
---|
| 114 | */ |
---|
| 115 | |
---|
| 116 | *n_weights = mapper->nWeights; |
---|
| 117 | |
---|
| 118 | } |
---|
| 119 | |
---|
| 120 | extern "C" void remap_get_barycentres_and_areas(const double* bounds_lon, const double* bounds_lat, int n_vert_per_cell, int n_cell, |
---|
| 121 | const double* pole, |
---|
| 122 | double* centre_lon, double* centre_lat, double* areas) |
---|
| 123 | { |
---|
| 124 | for (int i = 0; i < n_cell; i++) |
---|
| 125 | { |
---|
| 126 | int offs = i*n_vert_per_cell; |
---|
| 127 | Elt elt(bounds_lon + offs, bounds_lat + offs, n_vert_per_cell); |
---|
| 128 | cptEltGeom(elt, Coord(pole[0], pole[1], pole[2])); |
---|
| 129 | if (centre_lon && centre_lat) lonlat(elt.x, centre_lon[i], centre_lat[i]); |
---|
| 130 | if (areas) areas[i] = elt.area; |
---|
| 131 | } |
---|
| 132 | } |
---|
| 133 | |
---|
| 134 | /* |
---|
| 135 | extern "C" void remap_get_weights(double* weights, int* src_indices, int* src_ranks, int* dst_indices) |
---|
| 136 | { |
---|
| 137 | memcpy(weights, mapper->remapMatrix, mapper->nWeights*sizeof(double)); |
---|
| 138 | memcpy(src_indices, mapper->srcAddress, mapper->nWeights*sizeof(int)); |
---|
| 139 | memcpy(src_ranks, mapper->srcRank, mapper->nWeights*sizeof(int)); |
---|
| 140 | memcpy(dst_indices, mapper->dstAddress, mapper->nWeights*sizeof(int)); |
---|
| 141 | delete mapper; |
---|
| 142 | } |
---|
| 143 | */ |
---|
| 144 | |
---|
| 145 | extern "C" void remap_get_weights(double* weights, int* src_indices, int* dst_indices) |
---|
| 146 | { |
---|
| 147 | memcpy(weights, mapper->remapMatrix, mapper->nWeights*sizeof(double)); |
---|
| 148 | memcpy(src_indices, mapper->srcAddress, mapper->nWeights*sizeof(int)); |
---|
| 149 | memcpy(dst_indices, mapper->dstAddress, mapper->nWeights*sizeof(int)); |
---|
| 150 | delete mapper; |
---|
| 151 | } |
---|
| 152 | |
---|
| 153 | extern "C" void mpi_init() |
---|
| 154 | { |
---|
| 155 | /*int argc = 0; |
---|
| 156 | char **argv = NULL; |
---|
| 157 | MPI_Init(&argc, &argv);*/ |
---|
[1155] | 158 | //MPI_Init(NULL, NULL); |
---|
| 159 | int provided; |
---|
| 160 | MPI_Init_thread(NULL, NULL, 3, &provided); |
---|
| 161 | assert(provided >= 3); |
---|
[688] | 162 | } |
---|
| 163 | |
---|
| 164 | extern "C" int mpi_rank() |
---|
| 165 | { |
---|
| 166 | int rank; |
---|
| 167 | MPI_Comm_rank(MPI_COMM_WORLD, &rank); |
---|
| 168 | return rank; |
---|
| 169 | } |
---|
| 170 | |
---|
| 171 | extern "C" int mpi_size() |
---|
| 172 | { |
---|
| 173 | int size; |
---|
| 174 | MPI_Comm_size(MPI_COMM_WORLD, &size); |
---|
| 175 | return size; |
---|
| 176 | } |
---|
| 177 | |
---|
| 178 | extern "C" void mpi_finalize() |
---|
| 179 | { |
---|
| 180 | MPI_Finalize(); |
---|
| 181 | } |
---|