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

Last change on this file since 1205 was 1205, checked in by yushan, 4 years ago

branch merged with trunk @1200

File size: 5.6 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 <stdio.h>
17
18using namespace sphereRemap ;
19
20/* mapper is a ponter to a global class instance whoes members are allocated in the first step (finding the sizes of the weight arrays)
21   and deallocated during the second step (computing the weights) */
22Mapper *mapper;
23#pragma omp threadprivate(mapper)
24
25/** xxx_bounds_yyy is of length n_vert_per_cell_xxx*n_cell_xxx
26   This function computes the weights and returns number of weights is returned through the last argument.
27   To get the actuall weights, use `remap_get_weights`.
28*/
29extern "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,
30                     const double* src_pole,
31                     const double* dst_bounds_lon, const double* dst_bounds_lat, int n_vert_per_cell_dst, int n_cell_dst,
32                     const double* dst_pole,
33                     int order, int* n_weights)
34{
35  assert(src_bounds_lon);
36  assert(src_bounds_lat);
37  assert(n_vert_per_cell_src >= 3);
38  assert(n_cell_src >= 4);
39  assert(dst_bounds_lon);
40  assert(dst_bounds_lat);
41  assert(n_vert_per_cell_dst >= 3);
42  assert(n_cell_dst >= 4);
43  assert(1 <= order && order <= 2);
44
45  mapper = new Mapper(MPI_COMM_WORLD);
46  mapper->setVerbosity(PROGRESS) ;
47  mapper->setSourceMesh(src_bounds_lon, src_bounds_lat, n_vert_per_cell_src, n_cell_src, src_pole ) ;
48  mapper->setTargetMesh(dst_bounds_lon, dst_bounds_lat, 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
113extern "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/*
128extern "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
138extern "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
146extern "C" void mpi_init()
147{
148        /*int argc = 0;
149        char **argv = NULL;
150        MPI_Init(&argc, &argv);*/
151        //MPI_Init(NULL, NULL);
152        int provided;
153        MPI_Init_thread(NULL, NULL, 3, &provided);
154        assert(provided >= 3);
155}
156
157extern "C" int mpi_rank()
158{
159        int rank;
160        MPI_Comm_rank(MPI_COMM_WORLD, &rank);
161        return rank;
162}
163
164extern "C" int mpi_size()
165{
166        int size;
167        MPI_Comm_size(MPI_COMM_WORLD, &size);
168        return size;
169}
170
171extern "C" void mpi_finalize()
172{
173        MPI_Finalize();
174}
Note: See TracBrowser for help on using the repository browser.