![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 25494
- Timestamp:
- 08/27/20 11:45:27 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp
r25480 r25494 10 10 IssmDouble* pdouble, FILE* fid,char* field,int step); 11 11 int OutputStatistics(Parameters* parameters,Results* results); 12 int ComputeHistogram(Parameters* parameters,Results* results); 12 13 13 14 int main(int argc,char **argv){ /*{{{*/ … … 576 577 maxmeans=xNew<IssmDouble*>(nfields); 577 578 minmeans=xNew<IssmDouble*>(nfields); 578 meanxtype=xNew< IssmDouble*>(nfields);579 meanxsize=xNew< IssmDouble*>(nfields);579 meanxtype=xNew<int>(nfields); 580 meanxsize=xNew<int>(nfields); 580 581 581 582 maxxs=xNew<IssmDouble*>(nfields*nsteps); … … 597 598 _printf0_(" opening file:\n"); 598 599 FILE* fid=fopen(file,"rb"); 600 if(fid==NULL)_error_("cound not open file: " << file << "\n"); 599 601 600 602 /*figure out size of file, and read the whole thing:*/ … … 616 618 for (int f=0;f<nfields;f++){ 617 619 char* field=fields[f]; 620 fseek(fid,0,SEEK_SET); 618 621 for (int j=0;j<nsteps;j++){ 619 622 int counter=f*nsteps+j; … … 653 656 } 654 657 } 655 658 _printf0_(" average in time:\n"); 659 656 660 /*Deal with average in time: */ 661 fseek(fid,0,SEEK_SET); 657 662 for (int f=0;f<nfields;f++){ 658 663 char* field=fields[f]; 659 664 meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]); 660 if(i==(lower_row+1)){ 661 if(meanxtype[f]==1){ 662 maxmeans[f]=xNewZeroInit<IssmDouble>(1); 663 minmeans[f]=xNewZeroInit<IssmDouble>(1); 664 meanxsize[f]=1; 665 } 666 else{ 667 maxmeans[f]=xNewZeroInit<IssmDouble>(doublematsize); 668 minmeans[f]=xNewZeroInit<IssmDouble>(doublematsize); 669 meanxsize[f]=doublematsize; 670 } 671 672 } 665 673 666 if(meanxtype[f]==1){ 674 667 IssmDouble timemean=0; 668 fseek(fid,0,SEEK_SET); 675 669 for (int j=0;j<nsteps;j++){ 676 670 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]); 677 671 timemean+=scalar/nsteps; 678 672 } 673 679 674 /*Figure out max and min of time means: */ 680 *maxmeans[f]=max(*maxmeans[f],timemean); 681 *minmeans[f]=min(*minmeans[f],timemean); 675 if(i==(lower_row+1)){ 676 maxmeans[f]=xNewZeroInit<IssmDouble>(1); 677 minmeans[f]=xNewZeroInit<IssmDouble>(1); 678 *maxmeans[f]=timemean; 679 *minmeans[f]=timemean; 680 } 681 else{ 682 *maxmeans[f]=max(*maxmeans[f],timemean); 683 *minmeans[f]=min(*minmeans[f],timemean); 684 } 682 685 } 683 686 else{ 684 IssmDouble* maxx=maxmeans[f]; 685 IssmDouble* minx=minmeans[f]; 686 for(int k=0;k<doublematsize;k++){ 687 maxx[k]=max(maxx[k],doublemat[k]); 688 minx[k]=min(minx[k],doublemat[k]); 689 } 690 maxmeans[f]=maxx; 691 minmeans[f]=minx; 692 } 687 fseek(fid,0,SEEK_SET); 688 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize); 689 for (int j=0;j<nsteps;j++){ 690 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]); 691 for (int k=0;i<doublematsize;k++){ 692 timemean[k]+=doublemat[k]/nsteps; 693 } 694 } 695 696 if(i==(lower_row+1)){ 697 maxmeans[f]=timemean; 698 minmeans[f]=timemean; 699 } 700 else{ 701 IssmDouble* maxx=maxmeans[f]; 702 IssmDouble* minx=minmeans[f]; 703 704 for(int k=0;k<doublematsize;k++){ 705 maxx[k]=max(maxx[k],timemean[k]); 706 minx[k]=min(minx[k],timemean[k]); 707 } 708 maxmeans[f]=maxx; 709 minmeans[f]=minx; 710 } 711 } 693 712 } 694 713 fclose(fid); … … 715 734 IssmDouble sumscalar_alltimes=0; 716 735 717 ISSM_MPI_ Reduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());718 ISSM_MPI_ Reduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());736 ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); 737 ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); 719 738 720 739 /*add to results:*/ … … 727 746 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,j+1,0)); 728 747 } 748 749 /*Store broadcasted value for later computation of histograms:*/ 750 *maxxs[counter]=allmaxscalar; 751 *minxs[counter]=allminscalar; 752 729 753 } 730 754 } /*}}}*/ … … 736 760 737 761 /*we are broadcasting double arrays:*/ 738 IssmDouble* maxx= =maxxs[counter];739 IssmDouble* minx= =minxs[counter];762 IssmDouble* maxx=maxxs[counter]; 763 IssmDouble* minx=minxs[counter]; 740 764 741 765 IssmDouble* allmax=xNew<IssmDouble>(size); 742 766 IssmDouble* allmin=xNew<IssmDouble>(size); 743 767 744 ISSM_MPI_ Reduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());745 ISSM_MPI_ Reduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());768 ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); 769 ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); 746 770 747 771 /*add to results:*/ … … 754 778 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,j+1,0)); 755 779 } 780 /*Store broadcasted value for later computation of histograms:*/ 781 maxxs[counter]=allmax; 782 minxs[counter]=allmin; 756 783 } 757 784 } /*}}}*/ … … 768 795 IssmDouble allminscalar; 769 796 770 ISSM_MPI_ Reduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());771 ISSM_MPI_ Reduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());797 ISSM_MPI_Allreduce(&maxscalar,&allmaxscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); 798 ISSM_MPI_Allreduce(&minscalar,&allminscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); 772 799 773 800 /*add to results at time step 1:*/ … … 780 807 results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,allminscalar,1,0)); 781 808 } 809 /*Store for later use in histogram computation:*/ 810 *maxmeans[f]=allmaxscalar; 811 *minmeans[f]=allminscalar; 812 782 813 } /*}}}*/ 783 814 else{ /*deal with arrays:{{{*/ … … 786 817 787 818 /*we are broadcasting double arrays:*/ 788 IssmDouble* maxx= =maxmeans[f];789 IssmDouble* minx= =minmeans[f];819 IssmDouble* maxx=maxmeans[f]; 820 IssmDouble* minx=minmeans[f]; 790 821 791 822 IssmDouble* allmax=xNew<IssmDouble>(size); 792 823 IssmDouble* allmin=xNew<IssmDouble>(size); 793 824 794 ISSM_MPI_ Reduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm());795 ISSM_MPI_ Reduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,0,IssmComm::GetComm());825 ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); 826 ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); 796 827 797 828 /*add to results at step 1:*/ … … 800 831 801 832 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); 802 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1, j,0));833 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmax,size,1,1,0)); 803 834 sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin"); 804 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,j,0)); 805 } 835 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allmin,size,1,1,0)); 836 } 837 838 /*Store for later use in histogram computation:*/ 839 maxmeans[f]=allmax; 840 minmeans[f]=allmin; 841 806 842 } /*}}}*/ 807 843 } 808 844 809 845 /*Now that we have the min and max, we can start binning. First allocate 846 * histograms, then start filling them:*/ 847 IssmDouble** histogram=xNew<IssmDouble*>(nfields*nsteps); 848 IssmDouble** timehistogram=xNew<IssmDouble*>(nfields); 849 850 _printf0_("Start reading files again, this time binning values in the histogram:\n"); 851 /*Start opening files:*/ 852 for (int i=(lower_row+1);i<=upper_row;i++){ 853 _printf0_("reading file #: " << i << "\n"); 854 char file[1000]; 855 int length; 856 char* buffer=NULL; 857 858 /*string:*/ 859 sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); 860 861 /*open file: */ 862 _printf0_(" opening file:\n"); 863 FILE* fid=fopen(file,"rb"); 864 if(fid==NULL)_error_("cound not open file: " << file << "\n"); 865 866 /*figure out size of file, and read the whole thing:*/ 867 _printf0_(" reading file:\n"); 868 fseek (fid, 0, SEEK_END); 869 length = ftell (fid); 870 fseek (fid, 0, SEEK_SET); 871 buffer = xNew<char>(length); 872 fread (buffer, sizeof(char), length, fid); 873 874 /*close file:*/ 875 fclose (fid); 876 877 /*create a memory stream with this buffer:*/ 878 _printf0_(" processing file:\n"); 879 fid=fmemopen(buffer, length, "rb"); 880 881 /*start reading data from the buffer directly:*/ 882 for (int f=0;f<nfields;f++){ 883 char* field=fields[f]; 884 fseek(fid,0,SEEK_SET); 885 for (int j=0;j<nsteps;j++){ 886 int counter=f*nsteps+j; 887 xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]); 888 if(i==(lower_row+1)){ 889 if(xtype[counter]==1){ 890 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins); 891 IssmDouble ma=*maxxs[counter]; 892 IssmDouble mi=*minxs[counter]; 893 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; 894 localhistogram[index]++; 895 histogram[counter]=localhistogram; 896 } 897 else if (xtype[counter]==3){ 898 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins); 899 IssmDouble* ma=maxxs[counter]; 900 IssmDouble* mi=minxs[counter]; 901 for (int k=0;k<doublematsize;k++){ 902 IssmDouble scalar=doublemat[k]; 903 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1; 904 localhistogram[k*nbins+index]++; 905 } 906 histogram[counter]=localhistogram; 907 } 908 else _error_("cannot carry out statistics on type " << xtype[counter]); 909 } 910 else{ 911 if(xtype[counter]==1){ 912 IssmDouble* localhistogram=histogram[counter]; 913 IssmDouble ma=*maxxs[counter]; 914 IssmDouble mi=*minxs[counter]; 915 int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; 916 localhistogram[index]++; 917 } 918 else if (xtype[counter]==3){ 919 IssmDouble* localhistogram=histogram[counter]; 920 IssmDouble* ma=maxxs[counter]; 921 IssmDouble* mi=minxs[counter]; 922 for (int k=0;k<doublematsize;k++){ 923 IssmDouble scalar=doublemat[k]; 924 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1; 925 localhistogram[k*nbins+index]++; 926 } 927 } 928 else _error_("cannot carry out statistics on type " << xtype[counter]); 929 } 930 } 931 } 932 _printf0_(" average in time:\n"); 933 934 /*Deal with average in time: */ 935 fseek(fid,0,SEEK_SET); 936 for (int f=0;f<nfields;f++){ 937 char* field=fields[f]; 938 meanxtype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]); 939 940 if(meanxtype[f]==1){ 941 IssmDouble timemean=0; 942 fseek(fid,0,SEEK_SET); 943 for (int j=0;j<nsteps;j++){ 944 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]); 945 timemean+=scalar/nsteps; 946 } 947 948 /*Figure out max and min of time means: */ 949 if(i==(lower_row+1)){ 950 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(nbins); 951 IssmDouble ma=*maxmeans[f]; 952 IssmDouble mi=*minmeans[f]; 953 int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; 954 localhistogram[index]++; 955 timehistogram[f]=localhistogram; 956 } 957 else{ 958 IssmDouble* localhistogram=timehistogram[f]; 959 IssmDouble ma=*maxmeans[f]; 960 IssmDouble mi=*minmeans[f]; 961 int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; 962 localhistogram[index]++; 963 } 964 } 965 else{ 966 fseek(fid,0,SEEK_SET); 967 IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize); 968 for (int j=0;j<nsteps;j++){ 969 readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]); 970 for (int k=0;i<doublematsize;k++){ 971 timemean[k]+=doublemat[k]/nsteps; 972 } 973 } 974 975 if(i==(lower_row+1)){ 976 IssmDouble* localhistogram=xNewZeroInit<IssmDouble>(doublematsize*nbins); 977 IssmDouble* ma=maxmeans[f]; 978 IssmDouble* mi=minmeans[f]; 979 980 for (int k=0;k<doublematsize;k++){ 981 IssmDouble scalar=timemean[k]; 982 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1; 983 localhistogram[k*nbins+index]++; 984 } 985 timehistogram[f]=localhistogram; 986 } 987 else{ 988 989 IssmDouble* localhistogram=timehistogram[f]; 990 IssmDouble* ma=maxmeans[f]; 991 IssmDouble* mi=minmeans[f]; 992 993 for (int k=0;k<doublematsize;k++){ 994 IssmDouble scalar=timemean[k]; 995 int index=(scalar-mi[k])/(ma[k]-mi[k])*nbins; if (index==nbins)index=nbins-1; 996 localhistogram[k*nbins+index]++; 997 } 998 } 999 } 1000 } 1001 fclose(fid); 1002 1003 /*delete buffer:*/ 1004 xDelete<char>(buffer); 1005 } 1006 1007 /*We have agregated histograms across the cluster, now gather them across the cluster onto 1008 *cpu0: */ 1009 for (int f=0;f<nfields;f++){ 1010 int counter0=f*nsteps+0; 1011 if (xtype[counter0]==1){ /*deal with scalars {{{*/ 1012 for (int j=0;j<nsteps;j++){ 1013 int counter=f*nsteps+j; 1014 1015 /*we are broadcasting doubles:*/ 1016 IssmDouble* histo=histogram[counter]; //size nbins 1017 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins); 1018 1019 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); 1020 1021 /*add to results:*/ 1022 if(my_rank==0){ 1023 char fieldname[1000]; 1024 1025 sprintf(fieldname,"%s%s",fields[f],"Histogram"); 1026 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0)); 1027 } 1028 } 1029 } /*}}}*/ 1030 else{ /*deal with arrays:{{{*/ 1031 1032 int size=xsize[counter0]; 1033 for (int j=0;j<nsteps;j++){ 1034 int counter=f*nsteps+j; 1035 1036 /*we are broadcasting double arrays:*/ 1037 IssmDouble* histo=histogram[counter]; 1038 IssmDouble* allhisto=xNew<IssmDouble>(size*nbins); 1039 1040 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); 1041 1042 /*add to results:*/ 1043 if(my_rank==0){ 1044 char fieldname[1000]; 1045 1046 sprintf(fieldname,"%s%s",fields[f],"Histogram"); 1047 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0)); 1048 } 1049 } 1050 } /*}}}*/ 1051 } 1052 1053 /*Now do the same for the time mean fields:*/ 1054 for (int f=0;f<nfields;f++){ 1055 if (meanxtype[f]==1){ /*deal with scalars {{{*/ 1056 1057 /*we are broadcasting doubles:*/ 1058 IssmDouble* histo=timehistogram[f]; 1059 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(nbins); 1060 1061 ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); 1062 1063 /*add to results at time step 1:*/ 1064 if(my_rank==0){ 1065 char fieldname[1000]; 1066 1067 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); 1068 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,1,nbins,1,0)); 1069 } 1070 } /*}}}*/ 1071 else{ /*deal with arrays:{{{*/ 1072 1073 int size=meanxsize[f]; 1074 1075 /*we are broadcasting double arrays:*/ 1076 IssmDouble* histo=timehistogram[f]; 1077 IssmDouble* allhisto=xNewZeroInit<IssmDouble>(size*nbins); 1078 1079 ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); 1080 /*add to results at step 1:*/ 1081 if(my_rank==0){ 1082 char fieldname[1000]; 1083 1084 sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); 1085 results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allhisto,size,nbins,1,0)); 1086 } 1087 } /*}}}*/ 1088 } 810 1089 } /*}}}*/ 811 1090 int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/
Note:
See TracChangeset
for help on using the changeset viewer.
