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