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