Changeset 28262


Ignore:
Timestamp:
05/08/24 13:21:46 (11 months ago)
Author:
Mathieu Morlighem
Message:

CHG: making TF a lot less memory intensive to load by reading layer by layer instead of the whole 3D matrix

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r28155 r28262  
    203203                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
    204204                        iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.melt_anomaly",BasalforcingsIsmip6MeltAnomalyEnum,0.);
    205                         IssmDouble** array3d = NULL; int* Ms = NULL; int* Ns = NULL; int K;
    206                         iomodel->FetchData(&array3d,&Ms,&Ns,&K,"md.basalforcings.tf");
    207                         if(!array3d) _error_("md.basalforcings.tf not found in binary file");
    208                         for(Object* & object : elements->objects){
    209                                 Element*       element = xDynamicCast<Element*>(object);
    210                                 if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
    211                                 for(int kk=0;kk<K;kk++){
    212                                         element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array3d[kk],inputs,iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmip6TfEnum,kk);
     205
     206                        /*Deal with tf...*/
     207                        IssmDouble* array2d = NULL; int M,N,K; IssmDouble* temp = NULL;
     208                        iomodel->FetchData(&temp,&M,&K,"md.basalforcings.tf_depths"); xDelete<IssmDouble>(temp);
     209                        _assert_(M==1); _assert_(K>=1);
     210                        for(int kk=0;kk<K;kk++){
     211
     212                                /*Fetch TF for this depth*/
     213                                iomodel->FetchData(&array2d, &M, &N, kk, "md.basalforcings.tf");
     214                                if(!array2d) _error_("md.basalforcings.tf not found in binary file");
     215                                for(Object* & object : elements->objects){
     216                                        Element*  element = xDynamicCast<Element*>(object);
     217                                        if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
     218                                        element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array2d,inputs,iomodel,M,N,1,BasalforcingsIsmip6TfEnum,kk);
    213219                                }
    214                         }
    215                         xDelete<int>(Ms); xDelete<int>(Ns);
    216                         for(int i=0;i<K;i++) xDelete<IssmDouble>(array3d[i]);
    217                         xDelete<IssmDouble*>(array3d);
    218                         }
     220                        xDelete<IssmDouble>(array2d);
     221                        }
     222                                                                                          }
    219223                        break;
    220224                case BeckmannGoosseFloatingMeltRateEnum:
  • issm/trunk-jpl/src/c/classes/IoModel.cpp

    r28145 r28262  
    15651565}
    15661566/*}}}*/
     1567void  IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN, int layer_number,const char* data_name){/*{{{*/
     1568        /*Same function as above, but here we want to fetch only one specific "layer" of the dataset to avoid memory issues*/
     1569
     1570        /*output: */
     1571        int          M=0,N=0;
     1572        IssmPDouble *matrix = NULL;
     1573
     1574        /*intermediary: */
     1575        int  numrecords=0;
     1576        int  code,i;
     1577
     1578        /*recover my_rank:*/
     1579        int my_rank=IssmComm::GetRank();
     1580
     1581        /*Set file pointer to beginning of the data: */
     1582        fid=this->SetFilePointerToData(&code,NULL,data_name);
     1583        if(code!=8)_error_("expecting a IssmDouble mat array for \""<<data_name<<"\"");
     1584
     1585        if(my_rank==0){
     1586                if(fread(&numrecords,sizeof(int),1,fid)!=1) _error_("could not read number of records in matrix array ");
     1587        }
     1588        ISSM_MPI_Bcast(&numrecords,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     1589
     1590        if(numrecords){
     1591
     1592                /*Check consistency of layer number*/
     1593                _assert_(layer_number>=0);
     1594                if(layer_number>=numrecords) _error_("layer number of "<<data_name<<" cannot exceed "<< numrecords-1);
     1595
     1596                /*Skip until we get to the right layer*/
     1597                if(my_rank==0){
     1598                        for(i=0;i<layer_number;i++){
     1599                                if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
     1600                                if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns in " << i << "th matrix of matrix array");
     1601                                fseek(fid,M*N*sizeof(double),SEEK_CUR);
     1602                        }
     1603                }
     1604
     1605                /*fetch this slice!*/
     1606                if(my_rank==0){
     1607                        if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
     1608                }
     1609                ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     1610
     1611                if(my_rank==0){
     1612                        if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns in " << i << "th matrix of matrix array");
     1613                }
     1614                ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     1615
     1616                /*Now allocate matrix: */
     1617                if(M*N){
     1618                        matrix=xNew<IssmPDouble>(M*N);
     1619
     1620                        /*Read matrix on node 0, then broadcast: */
     1621                        if(my_rank==0) if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
     1622                        ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
     1623
     1624                        *pmatrix=xNew<IssmDouble>(M*N);
     1625                        for(int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
     1626                        xDelete<IssmPDouble>(matrix);
     1627                }
     1628                else{
     1629                        *pmatrix=NULL;
     1630                }
     1631        }
     1632
     1633        /*Assign output pointers: */
     1634        if(pM) *pM = M;
     1635        if(pN) *pN = N;
     1636}
     1637/*}}}*/
    15671638void  IoModel::FetchData(Options* options,const char* lastnonoption){/*{{{*/
    15681639
  • issm/trunk-jpl/src/c/classes/IoModel.h

    r27709 r28262  
    144144#endif
    145145                void        FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name);
     146                void        FetchData(IssmDouble** pmatrix,int* pM,int* pN, int layer_number,const char* data_name);
    146147                void        FetchData(Options *options,const char* data_name);
    147148                void        FetchData(int num,...);
Note: See TracChangeset for help on using the changeset viewer.