source: XIOS/dev/dev_cmip6_omp/extern/remap/test-main.cpp @ 1606

Last change on this file since 1606 was 694, checked in by rlacroix, 9 years ago

Fix compilation issues caused by the new "remap" library.

Use our MPI header so that C++ bindings are always disabled.

File size: 4.4 KB
Line 
1#include "mpi.hpp"
2#include <netcdf.h>
3#include "errhandle.hpp"
4#include "libmapper.hpp"
5#include "mapper.hpp"
6
7#include "node.hpp"
8
9using namespace sphereRemap ;
10
11int ncread(double *lon, double *lat, double* val, int nElt, const char* filename)
12{
13        int ncid, blonid, blatid, valid, neltid;
14        size_t nEltCheck;
15        exit_on_failure(nc_open(filename, NC_NOWRITE, &ncid), std::string("Cannot open netCDF file ") + filename);
16        exit_on_failure(nc_inq_dimid(ncid, "elt", &neltid), std::string("No dimension elt in file ") + filename);
17        nc_inq_dimlen(ncid, neltid, &nEltCheck);
18        exit_on_failure(nElt != nEltCheck, std::string("Array sizes do not match!"));
19
20        nc_inq_varid(ncid, "bounds_lon", &blonid);
21        nc_inq_varid(ncid, "bounds_lat", &blatid);
22        nc_inq_varid(ncid, "val", &valid);
23        nc_get_var_double(ncid, blonid, lon);
24        nc_get_var_double(ncid, blatid, lat);
25        nc_get_var_double(ncid, valid, val);
26        nc_close(ncid);
27}
28
29int ncwriteValue(double *val, const char* filename)
30{
31  int ncid, valid ;
32  nc_open(filename, NC_WRITE, &ncid);
33  nc_inq_varid(ncid, "val", &valid);
34  nc_put_var_double(ncid, valid, val);
35  nc_close(ncid);
36}
37 
38
39
40void compute_distribution(int nGlobalElts, int &start, int &nLocalElts)
41{
42        int mpiSize, mpiRank;
43        MPI_Comm_size(MPI_COMM_WORLD, &mpiSize);
44        MPI_Comm_rank(MPI_COMM_WORLD, &mpiRank);
45        start = 0;
46        nLocalElts = 0;
47        for (int i = 0; i <= mpiRank; i++)
48        {
49                start += nLocalElts;
50                nLocalElts = nGlobalElts/mpiSize;
51                if (i < nGlobalElts % mpiSize) nLocalElts++;
52        }
53}
54
55int main()
56{
57        int interpOrder = 2;
58/* low resolution */
59/*
60        char srcFile[] = "h14.nc";
61        char dstFile[] = "r180x90.nc";
62        double srcPole[] = {0, 0, 0};
63        double dstPole[] = {0, 0, 1};
64        int nSrcElt = 13661;
65        int nDstElt = 16200;
66*/
67/* high resolution */ 
68
69        char srcFile[] = "t740.nc";
70        char dstFile[] = "r1440x720.nc";
71        double srcPole[] = {0, 0, 0};
72        double dstPole[] = {0, 0, 1};
73        int nSrcElt = 741034;
74        int nDstElt = 1036800;
75
76        int nVert = 10;
77        int nSrcLocal, nDstLocal, startSrc, startDst;
78        int nWeight;
79  int mpi_rank, mpi_size ;
80
81        MPI_Init(NULL, NULL);
82  MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank) ;
83  MPI_Comm_size(MPI_COMM_WORLD, &mpi_size) ;
84
85        double *srcLon = (double *) malloc(nSrcElt*nVert*sizeof(double));       
86        double *srcLat = (double *) malloc(nSrcElt*nVert*sizeof(double));       
87        double *srcVal = (double *) malloc(nSrcElt*sizeof(double));     
88        double *dstLon = (double *) malloc(nDstElt*nVert*sizeof(double));       
89        double *dstLat = (double *) malloc(nDstElt*nVert*sizeof(double));       
90        double *dstVal = (double *) malloc(nDstElt*sizeof(double));
91        double *globalVal = (double *) malloc(nDstElt*sizeof(double));
92 
93        ncread(srcLon, srcLat, srcVal, nSrcElt, srcFile);
94        ncread(dstLon, dstLat, dstVal, nDstElt, dstFile);
95
96        compute_distribution(nSrcElt, startSrc, nSrcLocal);
97        compute_distribution(nDstElt, startDst, nDstLocal);
98
99
100
101  Mapper mapper(MPI_COMM_WORLD);
102  mapper.setVerbosity(PROGRESS) ;
103  mapper.setSourceMesh(srcLon + startSrc*nVert, srcLat + startSrc*nVert, nVert, nSrcLocal, srcPole ) ;
104  mapper.setTargetMesh(dstLon + startDst*nVert, dstLat + startDst*nVert, nVert, nDstLocal, dstPole ) ;
105
106  for(int i=0;i<nDstElt;i++) dstVal[i]=0 ;
107  mapper.setSourceValue(srcVal+ startSrc) ;
108 
109 
110  vector<double> timings = mapper.computeWeights(interpOrder);
111
112
113  for(int i=0;i<mapper.nWeights;i++)  dstVal[mapper.dstAddress[i]]=dstVal[mapper.dstAddress[i]]+mapper.remapMatrix[i]*srcVal[mapper.sourceWeightId[i]];
114
115/*
116        remap_get_num_weights(srcLon + startSrc*nVert, srcLat + startSrc*nVert, nVert, nSrcLocal, srcPole,
117                              dstLon + startDst*nVert, dstLat + startDst*nVert, nVert, nDstLocal, dstPole,
118                              interpOrder, &nWeight);
119
120        double *weights = (double *) malloc(nWeight*sizeof(double));   
121        int    *srcIdx  = (int *)    malloc(nWeight*sizeof(int));       
122        int    *srcRank = (int *)    malloc(nWeight*sizeof(int));       
123        int    *dstIdx  = (int *)    malloc(nWeight*sizeof(int));       
124
125        remap_get_weights(weights, srcIdx, dstIdx);
126
127        free(srcLon); free(srcLat);
128        free(dstLon); free(dstLat);
129        free(srcIdx); free(dstIdx);
130        free(srcRank);
131        free(weights);
132
133#ifdef DEBUG
134        memory_report();
135#endif
136*/
137  int* displs=new int[mpi_size] ;
138  int* recvCount=new int[mpi_size] ;
139 
140  MPI_Gather(&startDst,1,MPI_INT,displs,1,MPI_INT,0,MPI_COMM_WORLD) ;
141  MPI_Gather(&nDstLocal,1,MPI_INT,recvCount,1,MPI_INT,0,MPI_COMM_WORLD) ;
142  MPI_Gatherv(dstVal,nDstLocal,MPI_DOUBLE,globalVal,recvCount,displs,MPI_DOUBLE,0,MPI_COMM_WORLD) ;
143  if (mpi_rank==0)  ncwriteValue(globalVal, dstFile);
144 
145        MPI_Finalize();
146        return 0;
147}
Note: See TracBrowser for help on using the repository browser.