Changeset 28262
- Timestamp:
- 05/08/24 13:21:46 (11 months ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
r28155 r28262 203 203 iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.basin_id",BasalforcingsIsmip6BasinIdEnum); 204 204 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); 213 219 } 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 } 219 223 break; 220 224 case BeckmannGoosseFloatingMeltRateEnum: -
issm/trunk-jpl/src/c/classes/IoModel.cpp
r28145 r28262 1565 1565 } 1566 1566 /*}}}*/ 1567 void 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 /*}}}*/ 1567 1638 void IoModel::FetchData(Options* options,const char* lastnonoption){/*{{{*/ 1568 1639 -
issm/trunk-jpl/src/c/classes/IoModel.h
r27709 r28262 144 144 #endif 145 145 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); 146 147 void FetchData(Options *options,const char* data_name); 147 148 void FetchData(int num,...);
Note:
See TracChangeset
for help on using the changeset viewer.