Changeset 21668

Timestamp:

04/12/17 15:49:49 (8 years ago)

Author:

Mathieu Morlighem

Message:

NEW: new marshall format for SMB using uint8

Location:

issm/trunk-jpl/src

Files:

• c/classes/IoModel.cpp (modified) (6 diffs)
• m/classes/SMBforcing.m (modified) (1 diff)
• m/solve/WriteData.m (modified) (2 diffs)

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

@@ -683 +683 @@
                                 case 7: 
                                 case 8: 
+                                case 10:
                                         /*We are not interested in this record, too memory intensive. Skip it: */
                                         /*skip: */
@@ -825 +826 @@
                                         case 7: break; //do nothing. not interested in this type of data, which is memory intensive.
                                         case 8: break; //do nothing. not interested in this type of data, which is memory intensive.
+                                        case 10: break; //do nothing. not interested in this type of data, which is memory intensive.
                                         case 9:
                                                           ISSM_MPI_Bcast(&numstrings,1,ISSM_MPI_INT,0,IssmComm::GetComm()); 
@@ -1132 +1134 @@
         /*Set file pointer to beginning of the data: */
         fid=this->SetFilePointerToData(&code,NULL,data_name);
-        if(code!=5 && code!=6 && code!=7)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")");
+        if(code!=5 && code!=6 && code!=7 && code!=10)_error_("expecting a IssmDouble, integer or boolean matrix for \""<<data_name<<"\""<<" (Code is "<<code<<")");
 
         /*Now fetch: */
@@ -1150 +1152 @@
         /*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;
+                if(code==10){
+                        /*Special case for Compressed mat*/
+                        IssmPDouble offset,range;
+                        if(my_rank==0) if(fread(&offset,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read offset");
+                        ISSM_MPI_Bcast(&offset,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
+
+                        if(my_rank==0) if(fread(&range,sizeof(IssmPDouble),1,fid)!=1) _error_("could not read range");
+                        ISSM_MPI_Bcast(&range,1,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
+
+                        *pmatrix=xNew<IssmDouble>(M*N);
+
+                        /*Read matrix*/
+                        uint8_t* rawmatrix=xNew<uint8_t>((M-1)*N);
+                        if(my_rank==0) if(fread(rawmatrix,(M-1)*N*sizeof(char),1,fid)!=1) _error_("could not read matrix ");
+                        ISSM_MPI_Bcast(rawmatrix,(M-1)*N,ISSM_MPI_CHAR,0,IssmComm::GetComm()); 
+
+                        for(int i=0;i<(M-1)*N;++i) (*pmatrix)[i]=offset+range*reCast<IssmDouble>(rawmatrix[i])/255.;
+                        xDelete<uint8_t>(rawmatrix);
+
+                        /*read time now*/
+                        IssmPDouble* timematrix=xNew<IssmDouble>(N);
+                        if(my_rank==0) if(fread(timematrix,N*sizeof(IssmPDouble),1,fid)!=1) _error_("could not read time in compressed matrix");
+                        ISSM_MPI_Bcast(timematrix,N,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm()); 
+
+                        for(int i=0;i<N;++i) (*pmatrix)[(M-1)*N+i]=timematrix[i];
+                        xDelete<IssmPDouble>(timematrix);
+
+                }
+                else{
+                        /*Read matrix on node 0, then broadcast: */
+                        matrix=xNew<IssmPDouble>(M*N);
+                        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;
@@ -1542 +1572 @@
 
         /*Defaulting only supported for double arrays*/
-        if(code!=7) _error_(vector_name<<" is not a double array");
+        if(code!=7 && code!=10) _error_(vector_name<<" is not a double array");
 
         this->FetchData(&doublearray,&M,&N,vector_name);
@@ -2466 +2496 @@
 
                                 /*if record_code points to a vector, get its type (nodal or elementary): */
-                                if(5<=record_code && record_code<=7){
+                                if((5<=record_code && record_code<=7) || record_code==10){
                                         if(fread(&vector_type,sizeof(int),1,fid)!=1) _error_("Could not read vector_type");
                                 }

issm/trunk-jpl/src/m/classes/SMBforcing.m

@@ -69 +69 @@
                         WriteData(fid,prefix,'name','md.smb.model','data',1,'format','Integer');
                         WriteData(fid,prefix,'object',self,'class','smb','fieldname','mass_balance','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
+                        %WriteData(fid,prefix,'object',self,'class','smb','fieldname','mass_balance','format','CompressedMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
                         
                         %process requested outputs

issm/trunk-jpl/src/m/solve/WriteData.m

@@ -168 +168 @@
                 fwrite(fid,data','double'); %get to the "c" convention, hence the transpose
         end
+        % }}}
+elseif strcmpi(format,'CompressedMat'), % {{{
+
+        %Get size
+        s=size(data);
+
+        if s(1)==1 | s(2)==1,
+                %No need to use Compressed format
+                error('Not needed (should call WriteData with DoubleMat)');
+        end
+
+        %first write length of record
+        recordlength=4+4+8+8+1*(s(1)-1)*s(2)+8*s(2)+4+4; %2 integers (32 bits) + the matrix + code + matrix type
+        if recordlength>2^31; error(['field ' name ' cannot be marshalled because it is larger than 2^31 bytes!']); end
+        fwrite(fid,recordlength,'int');
+
+        %write data code and matrix type: 
+        fwrite(fid,FormatToCode(format),'int'); 
+        fwrite(fid,mattype,'int');
+
+        %write matrix size
+        fwrite(fid,s(1),'int'); 
+        fwrite(fid,s(2),'int'); 
+
+        %Write offset and range
+        A = data(1:end-1,:);
+        offset = min(A(:));
+        range = max(A(:)) - offset;
+        fwrite(fid,offset,'double'); 
+        fwrite(fid,range,'double'); 
+
+        %Convert data to uint8 and write it
+        A=uint8((A-offset)/range*255);
+        fwrite(fid,A','uint8'); %get to the "c" convention, hence the transpose
+
+        %Write last row as double (time)
+        fwrite(fid,data(end,:),'double'); 
         % }}}
 elseif strcmpi(format,'MatArray'), % {{{
@@ -258 +295 @@
         elseif strcmpi(format,'StringArray'),
                 code=9;
+        elseif strcmpi(format,'CompressedMat'),
+                code=10;
         else 
                 error('FormatToCode error message: data type not supported yet!');