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