Changeset 25496

Timestamp:

08/27/20 18:46:48 (5 years ago)

Author:

Eric.Larour

Message:

CHG: fixed issues with time mean histograms.

File:

• issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp (modified) (8 diffs)

issm/branches/trunk-larour-SLPS2020/src/c/main/issm_post.cpp

@@ -156 +156 @@
         int*         xsize=NULL;
 
+        IssmDouble** meanx=NULL;
+        IssmDouble** meanx2=NULL;
+        int*         meantype=NULL;
+        int*         meansize=NULL;
+
         /*Retrieve parameters:*/
         parameters->FindParam(&nsamples,QmuNsampleEnum);
@@ -175 +180 @@
         xtype=xNew<int>(nfields*nsteps);
         xsize=xNew<int>(nfields*nsteps);
+
+        meantype=xNew<int>(nfields);
+        meansize=xNew<int>(nfields);
+        meanx=xNew<IssmDouble*>(nfields);
+        meanx2=xNew<IssmDouble*>(nfields);
 
         /*Start opening files:*/
@@ -247 +257 @@
                         }
                 }
+                
+                /*Deal with time mean: */
+                for (int f=0;f<nfields;f++){
+                        char* field=fields[f];
+                        fseek(fid,0,SEEK_SET);
+                        meantype[f]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[0]);
+                        if(i==(lower_row+1)){
+                                if(meantype[f]==1){
+                                        meanx[f]=xNewZeroInit<IssmDouble>(1);
+                                        meanx2[f]=xNewZeroInit<IssmDouble>(1);
+                                        meansize[f]=1;
+                                }
+                                else{
+                                        meanx[f]=xNewZeroInit<IssmDouble>(doublematsize);
+                                        meanx2[f]=xNewZeroInit<IssmDouble>(doublematsize);
+                                        meansize[f]=doublematsize;
+                                }
+                        }
+                        fseek(fid,0,SEEK_SET);
+                        if(meantype[f]==1){
+                                IssmDouble sc=0;
+                                IssmDouble sc2=0;
+                                for(int j=0;j<nsteps;j++){
+                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+                                        sc+=scalar/nsteps;
+                                }
+                                sc2+=pow(sc,2.0);
+                                *meanx[f]+=sc;
+                                *meanx2[f]+=sc2;
+                        }
+                        else{
+                                IssmDouble* sc=meanx[f];
+                                IssmDouble* sc2=meanx2[f];
+                                IssmDouble* timemean=xNewZeroInit<IssmDouble>(doublematsize);
+                                IssmDouble* timemean2=xNewZeroInit<IssmDouble>(doublematsize);
+
+                                for(int j=0;j<nsteps;j++){
+                                        readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
+                                        for (int k=0;k<doublematsize;k++){
+                                                timemean[k]+=doublemat[k]/nsteps;
+                                        }
+                                }
+                                for (int k=0;k<doublematsize;k++){
+                                        timemean2[k]=pow(timemean[k],2.0);
+                                }
+                                for (int k=0;k<doublematsize;k++){
+                                        sc[k]+=timemean[k];
+                                        sc2[k]+=timemean2[k];
+                                }
+
+                        }
+
+                }
                 fclose(fid);
 
@@ -258 +321 @@
          *cpu0 and then compute statistics:*/
         for (int f=0;f<nfields;f++){
-
-                IssmDouble sumscalar_alltimes=0;
-                IssmDouble sumscalar2_alltimes=0;
-                IssmDouble* sumarray_alltimes=NULL;
-                IssmDouble* sumarray2_alltimes=NULL;
-                        
                 int counter0=f*nsteps+0;
                 if (xtype[counter0]==1){ /*deal with scalars {{{*/
@@ -293 +350 @@
                                 }
 
-                                /*We are also looking for the mean and stddev  over all time steps, 
-                                 *so accumulate sums again:*/
-                                sumscalar_alltimes+=sumscalar;
-                                sumscalar2_alltimes+=sumscalar2;
-                        }
-                        /*Build averae and standard deviation for all time steps:*/
-                        mean=sumscalar_alltimes/nsamples/nsteps;
-                        stddev=sqrt(sumscalar2_alltimes/nsamples/nsteps-pow(mean,2.0));
-                        /*add to results at step 1:*/
-                        if(my_rank==0){
-                                char fieldname[1000];
-
-                                sprintf(fieldname,"%s%s",fields[f],"TimeMean");
-                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
-                                sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
-                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
                         }
                 } /*}}}*/
@@ -314 +355 @@
 
                         int size=xsize[counter0];
-                        sumarray_alltimes=xNewZeroInit<IssmDouble>(size);
-                        sumarray2_alltimes=xNewZeroInit<IssmDouble>(size);
 
                         IssmDouble*  mean=xNew<IssmDouble>(size);
@@ -349 +388 @@
                                         results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
                                 }
-
-                                /*We want mean and stddev across time too, so accumulate through time: */
-                                for (int k=0;k<size;k++){
-                                        sumarray_alltimes[k]+=sumx[k];
-                                        sumarray2_alltimes[k]+=sumx2[k];
-                                }
-                        }
-
-                        /*build time mean and stddev:*/
+                        }
+                } /*}}}*/
+        }
+        /*Do the same but for the time mean:*/
+        for (int f=0;f<nfields;f++){
+                if (meantype[f]==1){ /*deal with scalars {{{*/
+                        IssmDouble mean,stddev;
+
+                        /*we are broadcasting doubles:*/
+                        IssmDouble scalar=*meanx[f];
+                        IssmDouble scalar2=*meanx2[f];
+                        IssmDouble sumscalar;
+                        IssmDouble sumscalar2;
+
+                        ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+                        ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+                        /*Build average and standard deviation. For standard deviation, use the 
+                         *following formula: sigma^2=E(x^2)-mu^2:*/
+                        mean=sumscalar/nsamples;
+                        stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
+
+                        /*add to results:*/
+                        if(my_rank==0){
+                                char fieldname[1000];
+
+                                sprintf(fieldname,"%s%s",fields[f],"TimeMean");
+                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,1,0));
+                                sprintf(fieldname,"%s%s",fields[f],"TimeStddev");
+                                results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,1,0));
+                        }
+                } /*}}}*/
+                else{ /*deal with arrays:{{{*/
+
+                        int size=meansize[f];
+                        IssmDouble*  mean=xNew<IssmDouble>(size);
+                        IssmDouble*  stddev=xNew<IssmDouble>(size);
+
+                        /*we are broadcasting double arrays:*/
+                        x=meanx[f];
+                        x2=meanx2[f];
+
+                        IssmDouble*  sumx=xNew<IssmDouble>(size);
+                        IssmDouble*  sumx2=xNew<IssmDouble>(size);
+                                
+                        ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+                        ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
+
+                        /*Build average and standard deviation. For standard deviation, use the 
+                         *following formula: sigma^2=E(x^2)-mu^2:*/
                         for (int k=0;k<size;k++){
-                                mean[k]=sumarray_alltimes[k]/nsamples/nsteps;
-                                stddev[k]=sqrt(sumarray2_alltimes[k]/nsamples/nsteps-pow(mean[k],2.0));
-                        }
-
-                        /*add to results at step 1:*/
+                                mean[k]=sumx[k]/nsamples;
+                                stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
+                        }
+
+                        /*add to results:*/
                         if(my_rank==0){
                                 char fieldname[1000];
@@ -373 +452 @@
                         }
                 } /*}}}*/
-
-        }
+        }
+
 
 } /*}}}*/