Changeset 28262

Timestamp:

05/08/24 13:21:46 (12 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:

• analyses/MasstransportAnalysis.cpp (modified) (1 diff)
• classes/IoModel.cpp (modified) (1 diff)
• classes/IoModel.h (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

@@ -203 +203 @@
                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum);
                         iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.melt_anomaly",BasalforcingsIsmip6MeltAnomalyEnum,0.);
-                        IssmDouble** array3d = NULL; int* Ms = NULL; int* Ns = NULL; int K;
-                        iomodel->FetchData(&array3d,&Ms,&Ns,&K,"md.basalforcings.tf");
-                        if(!array3d) _error_("md.basalforcings.tf not found in binary file");
-                        for(Object* & object : elements->objects){
-                                Element*       element = xDynamicCast<Element*>(object);
-                                if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
-                                for(int kk=0;kk<K;kk++){
-                                        element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array3d[kk],inputs,iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmip6TfEnum,kk);
+
+                        /*Deal with tf...*/
+                        IssmDouble* array2d = NULL; int M,N,K; IssmDouble* temp = NULL;
+                        iomodel->FetchData(&temp,&M,&K,"md.basalforcings.tf_depths"); xDelete<IssmDouble>(temp);
+                        _assert_(M==1); _assert_(K>=1);
+                        for(int kk=0;kk<K;kk++){
+
+                                /*Fetch TF for this depth*/
+                                iomodel->FetchData(&array2d, &M, &N, kk, "md.basalforcings.tf");
+                                if(!array2d) _error_("md.basalforcings.tf not found in binary file");
+                                for(Object* & object : elements->objects){
+                                        Element*  element = xDynamicCast<Element*>(object);
+                                        if(iomodel->domaintype!=Domain2DhorizontalEnum && !element->IsOnBase()) continue;
+                                        element->DatasetInputAdd(BasalforcingsIsmip6TfEnum,array2d,inputs,iomodel,M,N,1,BasalforcingsIsmip6TfEnum,kk);
                                 }
-                        }
-                        xDelete<int>(Ms); xDelete<int>(Ns);
-                        for(int i=0;i<K;i++) xDelete<IssmDouble>(array3d[i]);
-                        xDelete<IssmDouble*>(array3d);
-                        }
+                        xDelete<IssmDouble>(array2d);
+                        }
+                                                                                          }
                         break;
                 case BeckmannGoosseFloatingMeltRateEnum:

issm/trunk-jpl/src/c/classes/IoModel.cpp

@@ -1565 +1565 @@
 }
 /*}}}*/
+void  IoModel::FetchData(IssmDouble** pmatrix,int* pM,int* pN, int layer_number,const char* data_name){/*{{{*/
+        /*Same function as above, but here we want to fetch only one specific "layer" of the dataset to avoid memory issues*/
+
+        /*output: */
+        int          M=0,N=0;
+        IssmPDouble *matrix = NULL;
+
+        /*intermediary: */
+        int  numrecords=0;
+        int  code,i;
+
+        /*recover my_rank:*/
+        int my_rank=IssmComm::GetRank();
+
+        /*Set file pointer to beginning of the data: */
+        fid=this->SetFilePointerToData(&code,NULL,data_name);
+        if(code!=8)_error_("expecting a IssmDouble mat array for \""<<data_name<<"\"");
+
+        if(my_rank==0){
+                if(fread(&numrecords,sizeof(int),1,fid)!=1) _error_("could not read number of records in matrix array ");
+        }
+        ISSM_MPI_Bcast(&numrecords,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+        if(numrecords){
+
+                /*Check consistency of layer number*/
+                _assert_(layer_number>=0);
+                if(layer_number>=numrecords) _error_("layer number of "<<data_name<<" cannot exceed "<< numrecords-1);
+
+                /*Skip until we get to the right layer*/
+                if(my_rank==0){
+                        for(i=0;i<layer_number;i++){
+                                if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
+                                if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns in " << i << "th matrix of matrix array");
+                                fseek(fid,M*N*sizeof(double),SEEK_CUR);
+                        }
+                }
+
+                /*fetch this slice!*/
+                if(my_rank==0){
+                        if(fread(&M,sizeof(int),1,fid)!=1) _error_("could not read number of rows in " << i << "th matrix of matrix array");
+                }
+                ISSM_MPI_Bcast(&M,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+                if(my_rank==0){
+                        if(fread(&N,sizeof(int),1,fid)!=1) _error_("could not read number of columns in " << i << "th matrix of matrix array");
+                }
+                ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+
+                /*Now allocate matrix: */
+                if(M*N){
+                        matrix=xNew<IssmPDouble>(M*N);
+
+                        /*Read matrix on node 0, then broadcast: */
+                        if(my_rank==0) if(fread(matrix,M*N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read matrix ");
+                        ISSM_MPI_Bcast(matrix,M*N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
+
+                        *pmatrix=xNew<IssmDouble>(M*N);
+                        for(int i=0;i<M*N;++i) (*pmatrix)[i]=matrix[i];
+                        xDelete<IssmPDouble>(matrix);
+                }
+                else{
+                        *pmatrix=NULL;
+                }
+        }
+
+        /*Assign output pointers: */
+        if(pM) *pM = M;
+        if(pN) *pN = N;
+}
+/*}}}*/
 void  IoModel::FetchData(Options* options,const char* lastnonoption){/*{{{*/
 

issm/trunk-jpl/src/c/classes/IoModel.h

@@ -144 +144 @@
 #endif
                 void        FetchData(IssmDouble*** pmatrixarray,int** pmdims,int** pndims, int* pnumrecords,const char* data_name);
+                void        FetchData(IssmDouble** pmatrix,int* pM,int* pN, int layer_number,const char* data_name);
                 void        FetchData(Options *options,const char* data_name);
                 void        FetchData(int num,...);