source: XIOS/dev/branch_yushan_merged/extern/remap/src/libmapper.cpp @ 1176

Last change on this file since 1176 was 1176, checked in by yushan, 5 years ago

save modif

File size: 6.0 KB
Line 
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
20using 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) */
31Mapper *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*/
38extern "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
123extern "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/*
139extern "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
149extern "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
158extern "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
169extern "C" int mpi_rank()
170{
171        int rank;
172        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
173        return rank;
174}
175
176extern "C" int mpi_size()
177{
178        int size;
179        MPI_Comm_size(MPI_COMM_WORLD, &size);
180        return size;
181}
182
183extern "C" void mpi_finalize()
184{
185        MPI_Finalize();
186}
Note: See TracBrowser for help on using the repository browser.