Changeset 25361


Ignore:
Timestamp:
08/08/20 23:25:47 (5 years ago)
Author:
Eric.Larour
Message:

CHG: statistics core now working for average/stddev and for sample retrievals.

Location:
issm/branches/trunk-larour-SLPS2020/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.cpp

    r25317 r25361  
    5454
    5555/*Object constructors and destructor*/
     56FemModel::FemModel(void){ /*{{{*/
     57        /*do nothing:*/
     58} /*}}}*/
    5659FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm,bool trace){/*{{{*/
    5760

  • issm/branches/trunk-larour-SLPS2020/src/c/classes/FemModel.h

    r25066 r25361  
    6767
    6868                /*constructors, destructors: */
     69                FemModel(void);
    6970                FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false);
    7071                FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X);

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

    r25350 r25361  
    33 */
    44
     5/*includes and prototypes:*/
    56#include "./issm.h"
    6 
    7 int main(int argc,char **argv){
     7int ComputeMeanVariance(Parameters* parameters,Results* results);
     8int ComputeSampleSeries(Parameters* parameters,Results* results);
     9int readdata(IssmDouble** pdoublemat, int* pdoublematsize,
     10                     IssmDouble* pdouble, FILE* fid,char* field,int step);
     11int OutputStatistics(Parameters* parameters,Results* results);
     12
     13int main(int argc,char **argv){ /*{{{*/
    814       
    915        char* method=NULL;
    1016
     17        int nfields;
     18        char*  string=NULL;
     19        char** fields=NULL;
     20        char*  field=NULL;
     21        int    offset;
     22        char*  stepstring=NULL;
     23        int    step1,step2;
     24        char*  pattern=NULL;
     25        int    nsteps;
     26        int*   steps=NULL;
     27        int    nindices;
     28        int*   indices=NULL;
     29       
     30        /*datasets*/
     31        Parameters  *parameters   = NULL;
     32        Results  *results   = NULL;
     33
    1134        /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
    12         ISSM_MPI_Comm comm_init=EnvironmentInit(argc,argv);
     35        ISSM_MPI_Comm comm=EnvironmentInit(argc,argv);
    1336
    1437        /*First things first, store the communicator, and set it as a global variable: */
    15         IssmComm::SetComm(comm_init);
     38        IssmComm::SetComm(comm);
     39
     40        /*Initialize and retrieve parameters:{{{*/
     41        parameters   = new Parameters();
     42        results   = new Results();
     43
     44        /*Method and Solution:*/
     45        method=argv[1];
     46        parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum));
     47        parameters->AddObject(new StringParam(AnalysisTypeEnum,method));
     48
     49        /*Directory:*/
     50        parameters->AddObject(new StringParam(DirectoryNameEnum,argv[2]));
     51
     52        /*Model name:*/
     53        parameters->AddObject(new StringParam(InputFileNameEnum,argv[3]));
     54
     55        /*Number of samples:*/
     56        parameters->AddObject(new IntParam(QmuNsampleEnum,atoi(argv[4])));
     57
     58        /*Retrieve fields:*/
     59        nfields=atoi(argv[5]);
     60        for(int i=0;i<nfields;i++){
     61                fields=xNew<char*>(nfields);
     62                string= argv[6+i];
     63                field=xNew<char>(strlen(string)+1);
     64                xMemCpy<char>(field,string,(strlen(string)+1));
     65                fields[i]=field;
     66        }
     67        parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields));
     68        offset=6+nfields;
     69
     70        /*Retrieve time steps:*/
     71        stepstring=argv[offset];
     72        pattern=strstr(stepstring,":");
     73        if (pattern==NULL){
     74                step1=atoi(stepstring);
     75                step2=step1;
     76        }
     77        else{
     78                step2=atoi(pattern+1);
     79                stepstring[pattern-stepstring]='\0';
     80                step1=atoi(stepstring);
     81        }
     82        nsteps=step2-step1+1;
     83        steps=xNew<int>(nsteps);
     84        for (int i=step1;i<=step2;i++)steps[i-step1]=i;
     85        parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps));
     86
     87        offset++;
     88        /*}}}*/
     89       
     90        /*Key off method:*/
     91        if(strcmp(method,"MeanVariance")==0){
     92
     93                ComputeMeanVariance(parameters,results);
     94        }
     95        else if(strcmp(method,"SampleSeries")==0){
     96
     97                /*Retrieve the vertex indices where we'll retrieve our sample time series:*/
     98                nindices=atoi(argv[offset]); offset++;
     99                indices=xNew<int>(nindices);
     100                for (int i=0;i<nindices;i++){
     101                        indices[i]=atoi(argv[offset+i]);
     102                }
     103                parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices));
     104               
     105                ComputeSampleSeries(parameters,results);
     106        }
     107
     108
     109
     110        else _error_("unknown method: " << method << "\n");
     111
     112        /*output results:*/
     113        ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n");
     114        OutputStatistics(parameters,results);
     115
     116        /*Finalize ISSM:*/
     117        ISSM_MPI_Finalize();
     118
     119        /*Return unix success: */
     120        return 0;
     121
     122} /*}}}*/
     123int ComputeMeanVariance(Parameters* parameters,Results* results){ /*{{{*/
     124
     125        int nsamples;
     126        char* directory=NULL;
     127        char* model=NULL;
     128        char** fields=NULL;
     129        int* steps=NULL;
     130        int nsteps;
     131        int nfields;
     132        int range,lower_row,upper_row;
     133       
     134        /*intermediary:*/
     135        IssmDouble* doublemat=NULL;
     136        int         doublematsize;
     137        IssmDouble scalar;
     138
     139        /*computation of average and variance itself:*/
     140        IssmDouble*  x = NULL;
     141        IssmDouble*  x2 = NULL;
     142        IssmDouble** xs = NULL;
     143        IssmDouble** xs2 = NULL;
     144        int*         xtype=NULL;
     145        int*         xsize=NULL;
     146
     147        /*Retrieve parameters:*/
     148        parameters->FindParam(&nsamples,QmuNsampleEnum);
     149        parameters->FindParam(&directory,DirectoryNameEnum);
     150        parameters->FindParam(&model,InputFileNameEnum);
     151        parameters->FindParam(&fields,&nfields,FieldsEnum);
     152        parameters->FindParam(&steps,&nsteps,StepsEnum);
    16153
    17154        /*Get rank:*/
    18155        int my_rank=IssmComm::GetRank();
    19156
    20         /*Initialize parameters:*/
    21         Parameters  *parameters   = new Parameters();
    22 
    23         /*Method:*/
    24         method=argv[1];
    25         parameters->AddObject(new StringParam(SolutionTypeEnum,argv[1]));
    26 
    27         /*Directory:*/
    28         parameters->AddObject(new StringParam(DirectoryEnum,argv[2]));
    29 
    30         /*Model name:*/
    31         parameters->AddObject(new StringParam(InputFileNameEnum,argv[3]));
    32 
    33         /*Key off method:*/
    34         if(strcmp(model,"MeanVariance")==0)
    35         if
    36         switch(method){
    37                 case
    38 
    39         /*Finalize ISSM:*/
    40         ISSM_MPI_Finalize();
    41 
    42         /*Return unix success: */
    43         return 0;
     157        /*Open files and read them complelety, in a distributed way:*/
     158        range=DetermineLocalSize(nsamples,IssmComm::GetComm());
     159        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
     160
     161        /*Initialize arrays:*/
     162        xs=xNew<IssmDouble*>(nfields*nsteps);
     163        xs2=xNew<IssmDouble*>(nfields*nsteps);
     164        xtype=xNew<int>(nfields*nsteps);
     165        xsize=xNew<int>(nfields*nsteps);
     166
     167        /*Start opening files:*/
     168        for (int i=(lower_row+1);i<=upper_row;i++){
     169                _printf0_("reading file #: " << i << "\n");
     170                char file[1000];
     171                int  length;
     172                char* buffer=NULL;
     173
     174                /*string:*/
     175                sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
     176
     177                /*open file: */
     178                _printf0_("    opening file:\n");
     179                FILE* fid=fopen(file,"rb");
     180
     181                /*figure out size of file, and read the whole thing:*/
     182                _printf0_("    reading file:\n");
     183                fseek (fid, 0, SEEK_END);
     184                length = ftell (fid);
     185                fseek (fid, 0, SEEK_SET);
     186                buffer = xNew<char>(length);
     187                fread (buffer, sizeof(char), length, fid);
     188
     189                /*close file:*/
     190                fclose (fid);
     191
     192                /*create a memory stream with this buffer:*/
     193                _printf0_("    processing file:\n");
     194                fid=fmemopen(buffer, length, "rb");
     195
     196                /*start reading data from the buffer directly:*/
     197                for (int f=0;f<nfields;f++){
     198                        char* field=fields[f];
     199                        for (int j=0;j<nsteps;j++){
     200                                int counter=f*nsteps+j;
     201                                xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     202                                if(i==(lower_row+1)){
     203                                        if(xtype[counter]==1){
     204                                                xs[counter]=xNew<IssmDouble>(1);
     205                                                xs2[counter]=xNew<IssmDouble>(1);
     206                                                *xs[counter]=scalar;
     207                                                *xs2[counter]=pow(scalar,2.0);
     208                                                xsize[counter]=1;
     209                                        }
     210                                        else if (xtype[counter]==3){
     211                                                IssmDouble* doublemat2=xNew<IssmDouble>(doublematsize);
     212                                                for(int k=0;k<doublematsize;k++)doublemat2[k]=pow(doublemat[k],2.0);
     213                                                xs[counter]=doublemat;
     214                                                xs2[counter]=doublemat2;
     215                                                xsize[counter]=doublematsize;
     216                                        }
     217                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     218                                }
     219                                else{
     220                                        if(xtype[counter]==1){
     221                                                *xs[counter]+=scalar;
     222                                                *xs2[counter]+=pow(scalar,2.0);
     223                                        }
     224                                        else if (xtype[counter]==3){
     225                                                IssmDouble* newdoublemat=xs[counter];
     226                                                IssmDouble* newdoublemat2=xs2[counter];
     227                                                for(int k=0;k<doublematsize;k++){
     228                                                        newdoublemat[k]+=doublemat[k];
     229                                                        newdoublemat2[k]+=pow(doublemat[k],2.0);
     230                                                }
     231                                                xs[counter]=newdoublemat;
     232                                                xs2[counter]=newdoublemat2;
     233                                        }
     234                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     235                                }
     236                        }
     237                }
     238                fclose(fid);
     239
     240                /*delete buffer:*/
     241                xDelete<char>(buffer);
     242        }
     243        ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD);
     244        _printf0_("Done reading files, now computing mean and variance.\n");
     245
     246        /*We have agregated x and x^2 across the cluster, now consolidate: */
     247        for (int f=0;f<nfields;f++){
     248                for (int j=0;j<nsteps;j++){
     249                        int counter=f*nsteps+j;
     250                        if (xtype[counter]==1){
     251                                /*we are broadcasting doubles:*/
     252                                IssmDouble scalar=*xs[counter];
     253                                IssmDouble scalar2=*xs2[counter];
     254                                IssmDouble sumscalar;
     255                                IssmDouble sumscalar2;
     256                                IssmDouble mean,stddev;
     257
     258                                ISSM_MPI_Reduce(&scalar,&sumscalar,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     259                                ISSM_MPI_Reduce(&scalar2,&sumscalar2,1,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     260                                /*Build average and standard deviation. For standard deviation, use the
     261                                 *following formula: sigma^2=E(x^2)-mu^2:*/
     262                                mean=sumscalar/nsamples;
     263                                stddev=sqrt(sumscalar2/nsamples-pow(mean,2.0));
     264
     265                                /*add to results:*/
     266                                if(my_rank==0){
     267                                        char fieldname[1000];
     268                                       
     269                                        sprintf(fieldname,"%s%s",fields[f],"Mean");
     270                                        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,mean,j+1,0));
     271                                        sprintf(fieldname,"%s%s",fields[f],"Stddev");
     272                                        results->AddResult(new GenericExternalResult<IssmDouble>(results->Size()+1,fieldname,stddev,j+1,0));
     273                                }
     274                        }
     275                        else{
     276                                /*we are broadcasting double arrays:*/
     277                                x=xs[counter];
     278                                x2=xs2[counter];
     279                                int size=xsize[counter];
     280                               
     281                                IssmDouble*  sumx=xNew<IssmDouble>(size);
     282                                IssmDouble*  sumx2=xNew<IssmDouble>(size);
     283                                IssmDouble*  mean=xNew<IssmDouble>(size);
     284                                IssmDouble*  stddev=xNew<IssmDouble>(size);
     285
     286                                ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     287                                ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm());
     288
     289                                /*Build average and standard deviation. For standard deviation, use the
     290                                 *following formula: sigma^2=E(x^2)-mu^2:*/
     291                                for (int k=0;k<size;k++){
     292                                        mean[k]=sumx[k]/nsamples;
     293                                        stddev[k]=sqrt(sumx2[k]/nsamples-pow(mean[k],2.0));
     294                                }
     295
     296                                /*add to results:*/
     297                                if(my_rank==0){
     298                                        char fieldname[1000];
     299
     300                                        sprintf(fieldname,"%s%s",fields[f],"Mean");
     301                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,mean,size,1,j+1,0));
     302                                        sprintf(fieldname,"%s%s",fields[f],"Stddev");
     303                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,stddev,size,1,j+1,0));
     304                                }
     305                        }
     306                }
     307        }
     308
     309
     310} /*}}}*/
     311int ComputeSampleSeries(Parameters* parameters,Results* results){ /*{{{*/
     312
     313        int nsamples;
     314        char* directory=NULL;
     315        char* model=NULL;
     316        char** fields=NULL;
     317        int* steps=NULL;
     318        int nsteps;
     319        int nfields;
     320        int range,lower_row,upper_row;
     321        int* indices=NULL;
     322        int  nindices;
     323       
     324        /*intermediary:*/
     325        IssmDouble* doublemat=NULL;
     326        int         doublematsize;
     327        IssmDouble scalar;
     328
     329        /*computation of average and variance itself:*/
     330        IssmDouble*  x = NULL;
     331        IssmDouble*  allx=NULL;
     332        IssmDouble** xs = NULL;
     333        int*         xtype=NULL;
     334        int*         xsize=NULL;
     335
     336        /*Retrieve parameters:*/
     337        parameters->FindParam(&nsamples,QmuNsampleEnum);
     338        parameters->FindParam(&directory,DirectoryNameEnum);
     339        parameters->FindParam(&model,InputFileNameEnum);
     340        parameters->FindParam(&fields,&nfields,FieldsEnum);
     341        parameters->FindParam(&steps,&nsteps,StepsEnum);
     342        parameters->FindParam(&indices,&nindices,IndicesEnum);
     343
     344        /*Get rank:*/
     345        int my_rank=IssmComm::GetRank();
     346
     347        /*Open files and read them complelety, in a distributed way:*/
     348        range=DetermineLocalSize(nsamples,IssmComm::GetComm());
     349        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm());
     350
     351        /*Initialize arrays:*/
     352        xs=xNew<IssmDouble*>(nfields*nsteps);
     353        xtype=xNew<int>(nfields*nsteps);
     354        xsize=xNew<int>(nfields*nsteps);
     355
     356        /*Start opening files:*/
     357        for (int i=(lower_row+1);i<=upper_row;i++){
     358                _printf0_("reading file #: " << i << "\n");
     359                char file[1000];
     360                int  length;
     361                char* buffer=NULL;
     362
     363                /*string:*/
     364                sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i);
     365
     366                /*open file: */
     367                _printf0_("    opening file:\n");
     368                FILE* fid=fopen(file,"rb");
     369
     370                /*figure out size of file, and read the whole thing:*/
     371                _printf0_("    reading file:\n");
     372                fseek (fid, 0, SEEK_END);
     373                length = ftell (fid);
     374                fseek (fid, 0, SEEK_SET);
     375                buffer = xNew<char>(length);
     376                fread (buffer, sizeof(char), length, fid);
     377
     378                /*close file:*/
     379                fclose (fid);
     380
     381                /*create a memory stream with this buffer:*/
     382                _printf0_("    processing file:\n");
     383                fid=fmemopen(buffer, length, "rb");
     384
     385                /*start reading data from the buffer directly:*/
     386                for (int f=0;f<nfields;f++){
     387                        char* field=fields[f];
     388                        for (int j=0;j<nsteps;j++){
     389                                int counter=f*nsteps+j;
     390                                xtype[counter]=readdata(&doublemat, &doublematsize, &scalar, fid,field,steps[j]);
     391                                if(i==(lower_row+1)){
     392                                        if(xtype[counter]==1){
     393                                                x=xNew<IssmDouble>(range);
     394                                                x[0]=scalar;
     395                                                xs[counter]=x;
     396                                                xsize[counter]=range;
     397                                        }
     398                                        else if (xtype[counter]==3){
     399                                                x=xNew<IssmDouble>(nindices*range);
     400                                                for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
     401                                                xs[counter]=x;
     402                                                xsize[counter]=range*nindices;
     403                                        }
     404                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     405                                }
     406                                else{
     407                                        if(xtype[counter]==1){
     408                                                x=xs[counter];
     409                                                x[i-(lower_row+1)]=scalar;
     410                                                xs[counter]=x;
     411                                        }
     412                                        else if (xtype[counter]==3){
     413                                                x=xs[counter];
     414                                                for(int k=0;k<nindices;k++)x[(i-(lower_row+1))*nindices+k]=doublemat[indices[k]-1];
     415                                                xs[counter]=x;
     416                                        }
     417                                        else _error_("cannot carry out statistics on type " << xtype[counter]);
     418                                }
     419                        }
     420                }
     421                fclose(fid);
     422
     423                /*delete buffer:*/
     424                xDelete<char>(buffer);
     425        }
     426        ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD);
     427        _printf0_("Done reading files, now assembling time series.\n");
     428
     429        for (int f=0;f<nfields;f++){
     430                for (int j=0;j<nsteps;j++){
     431                        int counter=f*nsteps+j;
     432                        if (xtype[counter]==1){
     433                                /*we are broadcasting range times doubles:*/
     434                                x=xs[counter];
     435                                allx=xNew<IssmDouble>(nsamples);
     436                                MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, 1, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
     437                                /*add to results:*/
     438                                if(my_rank==0){
     439                                        char fieldname[1000];
     440                                       
     441                                        sprintf(fieldname,"%s%s",fields[f],"Samples");
     442                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,1,j+1,0));
     443                                }
     444                        }
     445                        else{
     446                                /*we are broadcasting double arrays:*/
     447                                x=xs[counter];
     448                                allx=xNew<IssmDouble>(nsamples*nindices);
     449                               
     450                                MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm());
     451                               
     452                                /*add to results:*/
     453                                if(my_rank==0){
     454                                        char fieldname[1000];
     455                                        sprintf(fieldname,"%s%s",fields[f],"Samples");
     456                                        results->AddResult(new GenericExternalResult<IssmPDouble*>(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0));
     457                                }
     458                        }
     459                }
     460        }
     461
     462
     463} /*}}}*/
     464int OutputStatistics(Parameters* parameters,Results* results){ /*{{{*/
     465       
     466        char   outputfilename[1000];
     467        char* directory=NULL;
     468        char* model=NULL;
     469        char* method=NULL;
     470        int   nsamples;
     471        FemModel* femmodel=new FemModel();
     472
     473       
     474        /*Some parameters that will allow us to use the OutputResultsx module:*/
     475        parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false));
     476        parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true));
     477
     478        parameters->FindParam(&directory,DirectoryNameEnum);
     479        parameters->FindParam(&model,InputFileNameEnum);
     480        parameters->FindParam(&method,AnalysisTypeEnum);
     481        parameters->FindParam(&nsamples,QmuNsampleEnum);
     482        sprintf(outputfilename,"%s/%s-%s.stats-%i",directory,model,method,nsamples);
     483        parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
     484
     485        /*Call OutputResults module:*/
     486        femmodel->parameters=parameters;
     487        femmodel->results=results;
     488
     489        OutputResultsx(femmodel);
     490} /*}}}*/
     491int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/
     492
     493        int length;
     494        char fieldname[1000];
     495        int   fieldname_size;
     496        IssmDouble   rtime;
     497        int          rstep;
     498        int M,N;
     499
     500        //fields that we retrive:
     501        IssmDouble  dfield;
     502        char*       sfield    = NULL;
     503        IssmDouble* dmatfield = NULL;
     504        int*        imatfield = NULL;
     505                       
     506        //type of the returned field:
     507        int type;
     508        int found=0;
     509
     510        while(1){
     511
     512                size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid);
     513                if(ret_code != 1) break; //we are done.
     514
     515                fread(fieldname, sizeof(char), fieldname_size, fid);
     516                //_printf0_("fieldname: " << fieldname << "\n");
     517
     518                fread(&rtime, sizeof(IssmDouble), 1, fid);
     519                fread(&rstep, sizeof(int), 1, fid);
     520
     521                //check on field:
     522                if ((step==rstep) && (strcmp(field,fieldname)==0)){
     523
     524                        //ok, go read the result really:
     525                        fread(&type,sizeof(int),1,fid);
     526                        fread(&M,sizeof(int),1,fid);
     527                        if (type==1){
     528                                fread(&dfield,sizeof(IssmDouble),1,fid);
     529                        }
     530                        else if (type==2){
     531                                fread(&M,sizeof(int),1,fid);
     532                                sfield=xNew<char>(M);
     533                                fread(sfield,sizeof(char),M,fid);
     534                        }
     535                        else if (type==3){
     536                                fread(&N,sizeof(int),1,fid);
     537                                dmatfield=xNew<IssmDouble>(M*N);
     538                                fread(dmatfield,sizeof(IssmDouble),M*N,fid);
     539                        }
     540                        else if (type==4){
     541                                fread(&N,sizeof(int),1,fid);
     542                                imatfield=xNew<int>(M*N);
     543                                fread(imatfield,sizeof(int),M*N,fid);
     544                        }
     545                        else _error_("cannot read data of type " << type << "\n");
     546                        found=1;
     547                        break;
     548                }
     549                else{
     550                        //just skim to next results.
     551                        fread(&type,sizeof(int),1,fid);
     552                        fread(&M,sizeof(int),1,fid);
     553                        if (type==1){
     554                                fseek(fid,sizeof(IssmDouble),SEEK_CUR);
     555                        }
     556                        else if(type==2){
     557                                fseek(fid,M*sizeof(char),SEEK_CUR);
     558                        }
     559                        else if(type==3){
     560                                fread(&N,sizeof(int),1,fid);
     561                                fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR);
     562                        }
     563                        else if(type==4){
     564                                fread(&N,sizeof(int),1,fid);
     565                                fseek(fid,M*N*sizeof(int),SEEK_CUR);
     566                        }
     567                        else _error_("cannot read data of type " << type << "\n");
     568                }
     569        }
     570        if(found==0)_error_("cound not find " << field << " at step " << step  << "\n");
     571
     572        /*assign output pointers:*/
     573        *pdoublemat=dmatfield;
     574        *pdoublematsize=M*N;
     575        *pdouble=dfield;
     576
     577        /*return:*/
     578        return type;
    44579
    45580}
     581/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.