1 | #include "mpi.hpp" |
---|
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); |
---|
43 | double* src_area=NULL ; |
---|
44 | double* dst_area=NULL ; |
---|
45 | mapper = new Mapper(MPI_COMM_WORLD); |
---|
46 | mapper->setVerbosity(PROGRESS) ; |
---|
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 ) ; |
---|
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 | { |
---|
140 | memcpy(weights, mapper->remapMatrix, mapper->nWeights*sizeof(double)); |
---|
141 | memcpy(src_indices, mapper->srcAddress, mapper->nWeights*sizeof(int)); |
---|
142 | memcpy(dst_indices, mapper->dstAddress, mapper->nWeights*sizeof(int)); |
---|
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 | } |
---|