- Timestamp:
- 01/31/19 12:12:52 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
XIOS/dev/dev_trunk_omp/src/transformation/domain_algorithm_interpolate.cpp
r1630 r1646 20 20 #include "interpolate_domain.hpp" 21 21 #include "grid.hpp" 22 #ifdef _usingEP 22 23 using namespace ep_lib; 24 #endif 23 25 24 26 namespace xios { … … 32 34 std::map<int, int>& elementPositionInGridDst2AxisPosition, 33 35 std::map<int, int>& elementPositionInGridDst2DomainPosition) 36 TRY 34 37 { 35 38 std::vector<CDomain*> domainListDestP = gridDst->getDomains(); … … 42 45 return (new CDomainAlgorithmInterpolate(domainListDestP[domainDstIndex], domainListSrcP[domainSrcIndex], interpolateDomain)); 43 46 } 47 CATCH 44 48 45 49 bool CDomainAlgorithmInterpolate::registerTrans() 50 TRY 46 51 { 47 52 CGridTransformationFactory<CDomain>::registerTransformation(TRANS_INTERPOLATE_DOMAIN, create); 48 53 } 54 CATCH 49 55 50 56 CDomainAlgorithmInterpolate::CDomainAlgorithmInterpolate(CDomain* domainDestination, CDomain* domainSource, CInterpolateDomain* interpDomain) 51 57 : CDomainAlgorithmTransformation(domainDestination, domainSource), interpDomain_(interpDomain), writeToFile_(false), readFromFile_(false) 58 TRY 52 59 { 53 60 CContext* context = CContext::getCurrent(); … … 94 101 95 102 } 103 CATCH 96 104 97 105 /*! … … 99 107 */ 100 108 void CDomainAlgorithmInterpolate::computeRemap() 109 TRY 101 110 { 102 111 using namespace sphereRemap; … … 305 314 CArray<double,2> boundsLonSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 306 315 CArray<double,2> boundsLatSrcUnmasked(nVertexSrc,nSrcLocalUnmasked); 316 CArray<double,1> areaSrcUnmasked(nSrcLocalUnmasked); 317 307 318 long int * globalSrcUnmasked = new long int [nSrcLocalUnmasked]; 308 319 309 320 nSrcLocalUnmasked=0 ; 321 bool hasSrcArea=domainSrc_->hasArea && !domainSrc_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ; 322 double srcAreaFactor ; 323 if (hasSrcArea) srcAreaFactor=1./(domainSrc_->radius*domainSrc_->radius) ; 324 310 325 for (int idx=0 ; idx < nSrcLocal; idx++) 311 326 { … … 317 332 boundsLatSrcUnmasked(n,nSrcLocalUnmasked) = boundsLatSrc(n,idx) ; 318 333 } 334 if (hasSrcArea) areaSrcUnmasked(nSrcLocalUnmasked) = domainSrc_->areavalue(idx)*srcAreaFactor ; 319 335 globalSrcUnmasked[nSrcLocalUnmasked]=globalSrc[idx] ; 320 336 ++nSrcLocalUnmasked ; 321 337 } 322 338 } 323 339 324 340 325 341 int nDstLocalUnmasked = 0 ; … … 328 344 CArray<double,2> boundsLonDestUnmasked(nVertexDest,nDstLocalUnmasked); 329 345 CArray<double,2> boundsLatDestUnmasked(nVertexDest,nDstLocalUnmasked); 346 CArray<double,1> areaDstUnmasked(nDstLocalUnmasked); 347 330 348 long int * globalDstUnmasked = new long int [nDstLocalUnmasked]; 331 349 332 350 nDstLocalUnmasked=0 ; 351 bool hasDstArea=domainDest_->hasArea && !domainDest_->radius.isEmpty() && !interpDomain_->use_area.isEmpty() && interpDomain_->use_area==true ; 352 double dstAreaFactor ; 353 if (hasDstArea) dstAreaFactor=1./(domainDest_->radius*domainDest_->radius) ; 333 354 for (int idx=0 ; idx < nDstLocal; idx++) 334 355 { … … 340 361 boundsLatDestUnmasked(n,nDstLocalUnmasked) = boundsLatDest(n,idx) ; 341 362 } 363 if (hasDstArea) areaDstUnmasked(nDstLocalUnmasked) = domainDest_->areavalue(idx)*dstAreaFactor ; 342 364 globalDstUnmasked[nDstLocalUnmasked]=globalDst[idx] ; 343 365 ++nDstLocalUnmasked ; … … 345 367 } 346 368 347 mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 348 mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 369 double* ptAreaSrcUnmasked = NULL ; 370 if (hasSrcArea) ptAreaSrcUnmasked=areaSrcUnmasked.dataFirst() ; 371 372 double* ptAreaDstUnmasked = NULL ; 373 if (hasDstArea) ptAreaDstUnmasked=areaDstUnmasked.dataFirst() ; 374 375 mapper.setSourceMesh(boundsLonSrcUnmasked.dataFirst(), boundsLatSrcUnmasked.dataFirst(), ptAreaSrcUnmasked, nVertexSrc, nSrcLocalUnmasked, &srcPole[0], globalSrcUnmasked); 376 mapper.setTargetMesh(boundsLonDestUnmasked.dataFirst(), boundsLatDestUnmasked.dataFirst(), ptAreaDstUnmasked, nVertexDest, nDstLocalUnmasked, &dstPole[0], globalDstUnmasked); 349 377 350 378 std::vector<double> timings = mapper.computeWeights(orderInterp,renormalize,quantity); … … 400 428 401 429 } 430 CATCH 402 431 403 432 void CDomainAlgorithmInterpolate::processPole(std::map<int,std::vector<std::pair<int,double> > >& interMapValuePole, 404 433 int nbGlobalPointOnPole) 434 TRY 405 435 { 406 436 CContext* context = CContext::getCurrent(); … … 468 498 469 499 } 500 CATCH 470 501 471 502 /*! … … 473 504 */ 474 505 void CDomainAlgorithmInterpolate::computeIndexSourceMapping_(const std::vector<CArray<double,1>* >& dataAuxInputs) 506 TRY 475 507 { 476 508 if (readFromFile_) … … 481 513 } 482 514 } 515 CATCH 483 516 484 517 void CDomainAlgorithmInterpolate::writeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue) 518 TRY 485 519 { 486 520 writeInterpolationInfo(fileToReadWrite_, interpMapValue); 487 521 } 522 CATCH 488 523 489 524 void CDomainAlgorithmInterpolate::readRemapInfo() 525 TRY 490 526 { 491 527 std::map<int,std::vector<std::pair<int,double> > > interpMapValue; … … 494 530 exchangeRemapInfo(interpMapValue); 495 531 } 532 CATCH 496 533 497 534 void CDomainAlgorithmInterpolate::convertRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue) 535 TRY 498 536 { 499 537 CContext* context = CContext::getCurrent(); … … 520 558 } 521 559 } 560 CATCH 522 561 523 562 /*! … … 525 564 */ 526 565 void CDomainAlgorithmInterpolate::exchangeRemapInfo(std::map<int,std::vector<std::pair<int,double> > >& interpMapValue) 566 TRY 527 567 { 528 568 CContext* context = CContext::getCurrent(); … … 706 746 delete [] recvBuff; 707 747 } 748 CATCH 708 749 709 750 /*! Redefined some functions of CONetCDF4 to make use of them */ … … 716 757 int CDomainAlgorithmInterpolate::WriteNetCdf::addDimensionWrite(const StdString& name, 717 758 const StdSize size) 759 TRY 718 760 { 719 761 return CONetCDF4::addDimension(name, size); 720 762 } 763 CATCH 721 764 722 765 int CDomainAlgorithmInterpolate::WriteNetCdf::addVariableWrite(const StdString& name, nc_type type, 723 766 const std::vector<StdString>& dim) 767 TRY 724 768 { 725 769 return CONetCDF4::addVariable(name, type, dim); 726 770 } 771 CATCH 727 772 728 773 void CDomainAlgorithmInterpolate::WriteNetCdf::endDefinition() 774 TRY 729 775 { 730 776 CONetCDF4::definition_end(); 731 777 } 778 CATCH 732 779 733 780 void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<int,1>& data, const StdString& name, … … 735 782 const std::vector<StdSize>* start, 736 783 const std::vector<StdSize>* count) 784 TRY 737 785 { 738 786 CONetCDF4::writeData<int,1>(data, name, collective, record, start, count); 739 787 } 788 CATCH 740 789 741 790 void CDomainAlgorithmInterpolate::WriteNetCdf::writeDataIndex(const CArray<double,1>& data, const StdString& name, … … 743 792 const std::vector<StdSize>* start, 744 793 const std::vector<StdSize>* count) 794 TRY 745 795 { 746 796 CONetCDF4::writeData<double,1>(data, name, collective, record, start, count); 747 797 } 798 CATCH 748 799 749 800 /* … … 754 805 void CDomainAlgorithmInterpolate::writeInterpolationInfo(std::string& filename, 755 806 std::map<int,std::vector<std::pair<int,double> > >& interpMapValue) 807 TRY 756 808 { 757 809 CContext* context = CContext::getCurrent(); … … 806 858 std::vector<StdSize> count(1, localNbWeight); 807 859 860 int my_rank; 861 MPI_Comm_rank(client->intraComm, &my_rank); 862 863 #ifdef _usingEP 808 864 int my_rank_loc = client->intraComm->ep_comm_ptr->size_rank_info[1].first; 809 int my_rank = client->intraComm->ep_comm_ptr->size_rank_info[0].first; 810 865 #elif _usingMPI 866 int my_rank_loc = 0; 867 #endif 811 868 812 869 813 870 WriteNetCdf *netCdfWriter; 814 871 872 #ifdef _usingEP 815 873 MPI_Barrier_local(client->intraComm); 816 874 #endif 875 817 876 if(my_rank_loc==0) 818 877 { 819 info(100)<<"rank "<< my_rank <<" create weight info file"<< std::endl; 878 #pragma omp critical (_output) 879 { 880 info(100)<<"rank "<< my_rank <<" create weight info file"<< std::endl; 881 } 820 882 821 883 WriteNetCdf my_writer(filename, client->intraComm); 822 info(100)<<"rank "<< my_rank <<" file created"<< std::endl; 884 #pragma omp critical (_output) 885 { 886 info(100)<<"rank "<< my_rank <<" file created"<< std::endl; 887 } 823 888 netCdfWriter = &my_writer; 824 889 … … 827 892 netCdfWriter->addDimensionWrite("n_dst", n_dst); 828 893 netCdfWriter->addDimensionWrite("n_weight", globalNbWeight); 829 info(100)<<"rank "<< my_rank <<" addDimensionWrite : n_src, n_dst, n_weight"<< std::endl; 894 #pragma omp critical (_output) 895 { 896 info(100)<<"rank "<< my_rank <<" addDimensionWrite : n_src, n_dst, n_weight"<< std::endl; 897 } 830 898 831 899 std::vector<StdString> dims(1,"n_weight"); … … 836 904 netCdfWriter->addVariableWrite("weight", NC_DOUBLE, dims); 837 905 838 info(100)<<"rank "<< my_rank <<" addVariableWrite : src_idx, dst_idx, weight"<< std::endl; 906 #pragma omp critical (_output) 907 { 908 info(100)<<"rank "<< my_rank <<" addVariableWrite : src_idx, dst_idx, weight"<< std::endl; 909 } 839 910 840 911 // End of definition 841 912 netCdfWriter->endDefinition(); 842 info(100)<<"rank "<< my_rank <<" endDefinition"<< std::endl; 913 #pragma omp critical (_output) 914 { 915 info(100)<<"rank "<< my_rank <<" endDefinition"<< std::endl; 916 } 843 917 844 918 netCdfWriter->closeFile(); 845 info(100)<<"rank "<< my_rank <<" file closed"<< std::endl; 846 } 847 919 #pragma omp critical (_output) 920 { 921 info(100)<<"rank "<< my_rank <<" file closed"<< std::endl; 922 } 923 } 924 925 #ifdef _usingEP 848 926 MPI_Barrier_local(client->intraComm); 927 #endif 849 928 850 929 #pragma omp critical (write_weight_data) 851 930 { 852 931 // open file 853 info(100)<<"rank "<< my_rank <<" writing in weight info file"<< std::endl; 932 #pragma omp critical (_output) 933 { 934 info(100)<<"rank "<< my_rank <<" writing in weight info file"<< std::endl; 935 } 854 936 855 937 WriteNetCdf my_writer(filename, true, client->intraComm); 856 info(100)<<"rank "<< my_rank <<" file opened"<< std::endl; 938 #pragma omp critical (_output) 939 { 940 info(100)<<"rank "<< my_rank <<" file opened"<< std::endl; 941 } 857 942 netCdfWriter = &my_writer; 858 943 … … 864 949 netCdfWriter->writeDataIndex(weights, "weight", false, 0, &start, &count); 865 950 866 info(100)<<"rank "<< my_rank <<" WriteDataIndex : src_idx, dst_idx, weight"<< std::endl; 951 #pragma omp critical (_output) 952 { 953 info(100)<<"rank "<< my_rank <<" WriteDataIndex : src_idx, dst_idx, weight"<< std::endl; 954 } 867 955 } 868 956 869 957 netCdfWriter->closeFile(); 870 info(100)<<"rank "<< my_rank <<" file closed"<< std::endl; 958 #pragma omp critical (_output) 959 { 960 info(100)<<"rank "<< my_rank <<" file closed"<< std::endl; 961 } 871 962 872 963 } 873 964 965 #ifdef _usingEP 874 966 MPI_Barrier_local(client->intraComm); 875 876 877 } 967 #endif 968 969 970 } 971 CATCH 878 972 879 973 /*! … … 885 979 void CDomainAlgorithmInterpolate::readInterpolationInfo(std::string& filename, 886 980 std::map<int,std::vector<std::pair<int,double> > >& interpMapValue) 981 TRY 887 982 { 888 983 int ncid ; … … 946 1041 } 947 1042 } 1043 CATCH 948 1044 949 1045 void CDomainAlgorithmInterpolate::apply(const std::vector<std::pair<int,double> >& localIndex, … … 952 1048 std::vector<bool>& flagInitial, 953 1049 bool ignoreMissingValue, bool firstPass ) 1050 TRY 954 1051 { 955 1052 int nbLocalIndex = localIndex.size(); … … 993 1090 } 994 1091 } 1092 CATCH 995 1093 996 1094 void CDomainAlgorithmInterpolate::updateData(CArray<double,1>& dataOut) 1095 TRY 997 1096 { 998 1097 if (detectMissingValue) … … 1001 1100 size_t nbIndex=dataOut.numElements() ; 1002 1101 1003 for (int idx = 0; idx < nbIndex; ++idx) 1004 { 1005 if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan 1102 if (allMissing.numElements()>0) 1103 { 1104 for (int idx = 0; idx < nbIndex; ++idx) 1105 { 1106 if (allMissing(idx)) dataOut(idx) = defaultValue; // If all data source are nan then data destination must be nan 1107 } 1006 1108 } 1007 1109 … … 1013 1115 } 1014 1116 } 1015 1016 } 1117 CATCH 1118 1119 }
Note: See TracChangeset
for help on using the changeset viewer.