source: XIOS/dev/branch_openmp/extern/remap/src/libmapper.cpp @ 1335

Last change on this file since 1335 was 1335, checked in by yushan, 6 years ago

dev_omp

File size: 5.5 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
17using namespace ep_lib;
18
19using 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) */
23Mapper *mapper;
24#pragma omp threadprivate(mapper)
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*/
30extern "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{
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);
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);
83  mapper->setVerbosity(PROGRESS) ;
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
114extern "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/*
129extern "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
139extern "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
147extern "C" void mpi_init()
148{
149        /*int argc = 0;
150        char **argv = NULL;
151        MPI_Init(&argc, &argv);*/
152        MPI_Init(NULL, NULL);
153}
154
155extern "C" int mpi_rank()
156{
157        int rank;
158        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
159        return rank;
160}
161
162extern "C" int mpi_size()
163{
164        int size;
165        MPI_Comm_size(MPI_COMM_WORLD, &size);
166        return size;
167}
168
169extern "C" void mpi_finalize()
170{
171        MPI_Finalize();
172}
Note: See TracBrowser for help on using the repository browser.