Improvement to readmod

The readmod functions verify now the dimension of the model file. Therefore IFOS will abort, if the model file is to small or too big.
Moreover, the write_matrix_disk function is used to write the models to disk.

src/readmod.c

@@ -50,7 +50,7 @@ void readmod(float  **  rho, float **  pi, float **  u, float ** taus, float **
         eta[l]=DT/pts[l];
     }
     
-	   fprintf(FP,"\n...reading model information from modell-files...\n");
+	   fprintf(FP,"\n...reading model information from model-files...\n");
     
 	   /* read density and seismic velocities */
 	   /* ----------------------------------- */
@@ -144,15 +144,42 @@ void readmod(float  **  rho, float **  pi, float **  u, float ** taus, float **
         for (j=1;j<=NYG;j++){
             
             if(WAVETYPE==1||WAVETYPE==3){
+                
+                if(feof(fp_vp)){
+                    err("Model file VP is to small. Check dimensions NX*NY of file.");
+                }
                 fread(&vp, sizeof(float), 1, fp_vp);
+                
                 if (sw_Qp){
-                    fread(&qp, sizeof(float), 1, fp_qp);}
+                    
+                    if(feof(fp_qp)){
+                        err("Model file QP is to small. Check dimensions NX*NY of file.");
+                    }
+                    
+                    fread(&qp, sizeof(float), 1, fp_qp);
+                }
+            }
+            
+            if(feof(fp_vs)){
+                err("Model file VS is to small. Check dimensions NX*NY of file.");
             }
             
             fread(&vs, sizeof(float), 1, fp_vs);
+            
+            if(feof(fp_rho)){
+                err("Model file RHO is to small. Check dimensions NX*NY of file.");
+            }
+            
             fread(&rhov, sizeof(float), 1, fp_rho);
+            
             if (sw_Qs){
-                fread(&qs, sizeof(float), 1, fp_qs);}
+                
+                if(feof(fp_vs)){
+                    err("Model file QS is to small. Check dimensions NX*NY of file.");
+                }
+                
+                fread(&qs, sizeof(float), 1, fp_qs);
+            }
             
             /* only the PE which belongs to the current global gridpoint
              is saving model parameters in his local arrays */
@@ -179,51 +206,65 @@ void readmod(float  **  rho, float **  pi, float **  u, float ** taus, float **
         }
     }
     
-    
-    
+
     
     if(WAVETYPE==1||WAVETYPE==3){
+        
+        fread(&vp, sizeof(float), 1, fp_vp);
+        if(!feof(fp_vp)){
+            err("Model file VP is to big. Check dimensions NX*NY of file.");
+        }
         fclose(fp_vp);
-        if (sw_Qp)  fclose(fp_qp);
+        
+        if (sw_Qp) {
+            fread(&qp, sizeof(float), 1, fp_qp);
+            if(!feof(fp_qp)){
+                err("Model file QP is to big. Check dimensions NX*NY of file.");
+            }
+            fclose(fp_qp);
+        }
+    }
+    
+    fread(&vs, sizeof(float), 1, fp_vs);
+    if(!feof(fp_vs)){
+        err("Model file VS is to big. Check dimensions NX*NY of file.");
     }
     fclose(fp_vs);
+    
+    fread(&rho, sizeof(float), 1, fp_rho);
+    if(!feof(fp_rho)){
+        err("Model file RHO is to big. Check dimensions NX*NY of file.");
+    }
     fclose(fp_rho);
-    if (sw_Qs)  fclose(fp_qs);
+    
+    if (sw_Qs){
+        fread(&qs, sizeof(float), 1, fp_qs);
+        if(!feof(fp_qs)){
+            err("Model file QS is to big. Check dimensions NX*NY of file.");
+        }
+        fclose(fp_qs);
+    }
     
     
     /* each PE writes his model to disk */
     if(WAVETYPE==1||WAVETYPE==3){
         if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vp",MFILE);
         if(PARAMETERIZATION==3) sprintf(filename,"%s.out.pi",MFILE);
-        writemod(filename,pi,3);
-        MPI_Barrier(MPI_COMM_WORLD);
-        if (MYID==0) mergemod(filename,3);
-        MPI_Barrier(MPI_COMM_WORLD);
-        if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vp.%i.%i",MFILE,POS[1],POS[2]);
-        if(PARAMETERIZATION==3) sprintf(filename,"%s.out.pi.%i.%i",MFILE,POS[1],POS[2]);
-        remove(filename);
+        write_matrix_disk(pi, filename);
+        
+        sprintf(filename,"%s.out.qp",MFILE);
+        write_matrix_disk(taup, filename);
     }
     
     if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vs",MFILE);
     if(PARAMETERIZATION==3) sprintf(filename,"%s.out.mu",MFILE);
-    writemod(filename,u,3);
-    MPI_Barrier(MPI_COMM_WORLD);
-    
-    if (MYID==0) mergemod(filename,3);
-    MPI_Barrier(MPI_COMM_WORLD);
-    if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vs.%i.%i",MFILE,POS[1],POS[2]);
-    if(PARAMETERIZATION==3) sprintf(filename,"%s.out.mu.%i.%i",MFILE,POS[1],POS[2]);
-    remove(filename);
-    
+    write_matrix_disk(u, filename);
     
     sprintf(filename,"%s.out.rho",MFILE);
-    writemod(filename,rho,3);
-    MPI_Barrier(MPI_COMM_WORLD);
+    write_matrix_disk(rho, filename);
     
-    if (MYID==0) mergemod(filename,3);
-    MPI_Barrier(MPI_COMM_WORLD);
-    sprintf(filename,"%s.out.rho.%i.%i",MFILE,POS[1],POS[2]);
-    remove(filename);
+    sprintf(filename,"%s.out.qs",MFILE);
+    write_matrix_disk(taus, filename);
     
     free_vector(pts,1,L);
 }

src/readmod_acoustic.c

@@ -2,98 +2,105 @@
  * Copyright (C) 2016  For the list of authors, see file AUTHORS.
  *
  * This file is part of IFOS.
- * 
+ *
  * IFOS is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, version 2.0 of the License only.
- * 
+ *
  * IFOS is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  * GNU General Public License for more details.
- * 
+ *
  * You should have received a copy of the GNU General Public License
  * along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
------------------------------------------------------------------------------------------*/
+ -----------------------------------------------------------------------------------------*/
 
 /*------------------------------------------------------------------------
- *   Read acoustic model properties (vp,density) from files  
+ *   Read acoustic model properties (vp,density) from files
  *
  *  ----------------------------------------------------------------------*/
 
 
 /* This file contains function readmod, which has the purpose
-   to read data from model-files for viscoelastic simulation */
+ to read data from model-files for viscoelastic simulation */
 
 #include "fd.h"
 
 void readmod_acoustic(float  **  rho, float **  pi){
-
-	extern int NX, NY, NXG, NYG,  POS[3], MYID, PARAMETERIZATION;
-	extern char  MFILE[STRING_SIZE];	
-	extern FILE *FP;
-
-		
-	/* local variables */
-	float rhov, piv, vp;
-	int i, j, ii, jj;
-	FILE *fp_vp, *fp_rho;
-	char filename[STRING_SIZE];
-
-
-	   fprintf(FP,"\n...reading model information from modell-files...\n");
-           
+    
+    extern int NX, NY, NXG, NYG,  POS[3], MYID, PARAMETERIZATION;
+    extern char  MFILE[STRING_SIZE];
+    extern FILE *FP;
+    
+    
+    /* local variables */
+    float rhov, piv, vp;
+    int i, j, ii, jj;
+    FILE *fp_vp, *fp_rho;
+    char filename[STRING_SIZE];
+    
+    
+	   fprintf(FP,"\n...reading model information from model-files...\n");
+    
 	   /* read density and seismic velocities */
 	   /* ----------------------------------- */
 	   if(PARAMETERIZATION==1){
-	   fprintf(FP,"\t Vp:\n\t %s.vp\n\n",MFILE);
-	   sprintf(filename,"%s.vp",MFILE);
-	   fp_vp=fopen(filename,"r");
-	   if (fp_vp==NULL) err(" Could not open model file for Vp ! ");
-
-	   fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
-	   sprintf(filename,"%s.rho",MFILE);
-	   fp_rho=fopen(filename,"r");
-	   if (fp_rho==NULL) err(" Could not open model file for densities ! ");
-           }
+           fprintf(FP,"\t Vp:\n\t %s.vp\n\n",MFILE);
+           sprintf(filename,"%s.vp",MFILE);
+           fp_vp=fopen(filename,"r");
+           if (fp_vp==NULL) err(" Could not open model file for Vp ! ");
+           
+           fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
+           sprintf(filename,"%s.rho",MFILE);
+           fp_rho=fopen(filename,"r");
+           if (fp_rho==NULL) err(" Could not open model file for densities ! ");
+       }
 	   
-
-	/* loop over global grid */
-		for (i=1;i<=NXG;i++){
-			for (j=1;j<=NYG;j++){
-			fread(&vp, sizeof(float), 1, fp_vp);
-			fread(&rhov, sizeof(float), 1, fp_rho);
-				
-			/* only the PE which belongs to the current global gridpoint 
-			is saving model parameters in his local arrays */
-				if ((POS[1]==((i-1)/NX)) && 
-				    (POS[2]==((j-1)/NY))){
-					ii=i-POS[1]*NX;
-					jj=j-POS[2]*NY;
-                                
-                                pi[jj][ii]=vp;
-                                rho[jj][ii]=rhov;
-				
-				}
-			}
-		}
-	
-
-	fclose(fp_vp);
-	fclose(fp_rho);
-	
-	
-	/* each PE writes his model to disk */
-	sprintf(filename,"%s.fdveps.pi",MFILE);
-	writemod(filename,pi,3);
-	MPI_Barrier(MPI_COMM_WORLD);
-
-	if (MYID==0) mergemod(filename,3);
-	
-	sprintf(filename,"%s.fdveps.rho",MFILE);
-	writemod(filename,rho,3);
-	MPI_Barrier(MPI_COMM_WORLD);
-	                        
-	if (MYID==0) mergemod(filename,3);
-
+    
+    /* loop over global grid */
+    for (i=1;i<=NXG;i++){
+        for (j=1;j<=NYG;j++){
+            
+            if(feof(fp_vp) && feof(fp_rho)){
+                err("Model file VP or RHO is to small. Check dimensions NX*NY of file.");
+            }
+            
+            fread(&vp, sizeof(float), 1, fp_vp);
+            fread(&rhov, sizeof(float), 1, fp_rho);
+            
+            /* only the PE which belongs to the current global gridpoint
+             is saving model parameters in his local arrays */
+            if ((POS[1]==((i-1)/NX)) &&
+                (POS[2]==((j-1)/NY))){
+                ii=i-POS[1]*NX;
+                jj=j-POS[2]*NY;
+                
+                pi[jj][ii]=vp;
+                rho[jj][ii]=rhov;
+                
+            }
+        }
+    }
+    
+    
+    fread(&vp, sizeof(float), 1, fp_vp);
+    if(!feof(fp_vp)){
+        err("Model file VP is to big. Check dimensions NX*NY of file.");
+    }
+    fclose(fp_vp);
+    
+    fread(&rho, sizeof(float), 1, fp_rho);
+    if(!feof(fp_rho)){
+        err("Model file RHO is to big. Check dimensions NX*NY of file.");
+    }
+    fclose(fp_rho);
+    
+    /* Write model to file */
+    sprintf(filename,"%s.out.vp",MFILE);
+    write_matrix_disk(pi, filename);
+    
+    sprintf(filename,"%s.out.rho",MFILE);
+    write_matrix_disk(rho, filename);
+    
 }

src/readmod_elastic.c

@@ -43,7 +43,7 @@ void readmod_elastic(float  **  rho, float **  pi, float **  u){
     
     
     
-	   fprintf(FP,"\n...reading model information from modell-files...\n");
+	   fprintf(FP,"\n...reading model information from model-files...\n");
     
 	   /* read density and seismic velocities */
 	   /* ----------------------------------- */
@@ -96,9 +96,19 @@ void readmod_elastic(float  **  rho, float **  pi, float **  u){
         for (j=1;j<=NYG;j++){
             
             if(WAVETYPE==1||WAVETYPE==3){
+                
+                if(feof(fp_vp)){
+                    err("Model file VP is to small. Check dimensions NX*NY of file.");
+                }
+                
                 fread(&vp, sizeof(float), 1, fp_vp);
                 
             }
+            
+            if(feof(fp_vs) && feof(fp_rho)){
+                err("Model file VS or RHO is to small. Check dimensions NX*NY of file.");
+            }
+            
             fread(&vs, sizeof(float), 1, fp_vs);
             fread(&rhov, sizeof(float), 1, fp_rho);
             
@@ -122,8 +132,19 @@ void readmod_elastic(float  **  rho, float **  pi, float **  u){
         }
     }
     if(WAVETYPE==1||WAVETYPE==3){
+        
+        fread(&vp, sizeof(float), 1, fp_vp);
+        if(!feof(fp_vp)){
+            err("Model file VP is to big. Check dimensions NX*NY of file.");
+        }
         fclose(fp_vp);
     }
+    
+    fread(&vs, sizeof(float), 1, fp_vs);
+    fread(&rho, sizeof(float), 1, fp_rho);
+    if(!feof(fp_vs) && !feof(fp_rho)){
+        err("Model file VS or RHO is to big. Check dimensions NX*NY of file.");
+    }
     fclose(fp_vs);
     fclose(fp_rho);
     
@@ -133,36 +154,15 @@ void readmod_elastic(float  **  rho, float **  pi, float **  u){
     if(WAVETYPE==1||WAVETYPE==3){
         if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vp",MFILE);
         if(PARAMETERIZATION==3) sprintf(filename,"%s.out.pi",MFILE);
-        writemod(filename,pi,3);
-        MPI_Barrier(MPI_COMM_WORLD);
-        if (MYID==0) mergemod(filename,3);
-        MPI_Barrier(MPI_COMM_WORLD);
-        if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vp.%i.%i",MFILE,POS[1],POS[2]);
-        if(PARAMETERIZATION==3) sprintf(filename,"%s.out.pi.%i.%i",MFILE,POS[1],POS[2]);
-        remove(filename);
+        write_matrix_disk(pi, filename);
     }
     
     if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vs",MFILE);
     if(PARAMETERIZATION==3) sprintf(filename,"%s.out.mu",MFILE);
-    writemod(filename,u,3);
-    MPI_Barrier(MPI_COMM_WORLD);
-    
-    if (MYID==0) mergemod(filename,3);
-    MPI_Barrier(MPI_COMM_WORLD);
-    if(PARAMETERIZATION==1) sprintf(filename,"%s.out.vs.%i.%i",MFILE,POS[1],POS[2]);
-    if(PARAMETERIZATION==3) sprintf(filename,"%s.out.mu.%i.%i",MFILE,POS[1],POS[2]);
-    remove(filename);
-    
+    write_matrix_disk(u, filename);
     
     sprintf(filename,"%s.out.rho",MFILE);
-    writemod(filename,rho,3);
-    MPI_Barrier(MPI_COMM_WORLD);
-    
-    if (MYID==0) mergemod(filename,3);
-    MPI_Barrier(MPI_COMM_WORLD);
-    sprintf(filename,"%s.out.rho.%i.%i",MFILE,POS[1],POS[2]);
-    remove(filename);
-    
+    write_matrix_disk(rho, filename);
     
 }
 

src/readmod_viscac.c

@@ -50,7 +50,7 @@ void readmod_viscac(float  **  rho, float **  pi, float ** taup, float * eta){
 		eta[l]=DT/pts[l];
 	}
 	
-	fprintf(FP,"\n...reading model information from modell-files...\n");
+	fprintf(FP,"\n...reading model information from model-files...\n");
 	
 	/* read density and seismic velocities */
 	/* ----------------------------------- */
@@ -81,9 +81,20 @@ void readmod_viscac(float  **  rho, float **  pi, float ** taup, float * eta){
 	/* loop over global grid */
 	for (i=1;i<=NXG;i++){
 	for (j=1;j<=NYG;j++){
-		fread(&vp, sizeof(float), 1, fp_vp);
+        
+        if(feof(fp_vp) && feof(fp_rho)){
+            err("Model file VP or RHO is to small. Check dimensions NX*NY of file.");
+        }
+		
+        fread(&vp, sizeof(float), 1, fp_vp);
 		fread(&rhov, sizeof(float), 1, fp_rho);
-		if (sw_Qp){
+		
+        if (sw_Qp){
+            
+            if(feof(fp_qp)){
+                err("Model file QP is to small. Check dimensions NX*NY of file.");
+            }
+            
 			fread(&qp, sizeof(float), 1, fp_qp);
 		}
 		
@@ -104,24 +115,36 @@ void readmod_viscac(float  **  rho, float **  pi, float ** taup, float * eta){
 	}
 	}
 	
-	
-	fclose(fp_vp);
-	fclose(fp_rho);
-	if (sw_Qp)  fclose(fp_qp);
-	
-	
-	/* each PE writes his model to disk */
-	sprintf(filename,"%s.fdveps.pi",MFILE);
-	writemod(filename,pi,3);
-	MPI_Barrier(MPI_COMM_WORLD);
+    fread(&vp, sizeof(float), 1, fp_vp);
+    if(!feof(fp_vp)){
+        err("Model file VP is to big. Check dimensions NX*NY of file.");
+    }
+    fclose(fp_vp);
+    
+    fread(&rho, sizeof(float), 1, fp_rho);
+    if(!feof(fp_rho)){
+        err("Model file RHO is to big. Check dimensions NX*NY of file.");
+    }
+    fclose(fp_rho);
 
-	if (MYID==0) mergemod(filename,3);
+    if (sw_Qp){
+        fread(&qp, sizeof(float), 1, fp_qp);
+        if(!feof(fp_qp)){
+            err("Model file QP is to big. Check dimensions NX*NY of file.");
+        }
+        fclose(fp_qp);
+    }
+    
 	
-	sprintf(filename,"%s.fdveps.rho",MFILE);
-	writemod(filename,rho,3);
-	MPI_Barrier(MPI_COMM_WORLD);
-	                        
-	if (MYID==0) mergemod(filename,3);
-
+	/* write model to disk */
+	sprintf(filename,"%s.out.vp",MFILE);
+    write_matrix_disk(pi,filename);
+    
+    sprintf(filename,"%s.out.rho",MFILE);
+    write_matrix_disk(rho,filename);
+    
+    sprintf(filename,"%s.out.qp",MFILE);
+    write_matrix_disk(taup,filename);
+    
 	free_vector(pts,1,L);
 }