readmod_visco.c 4.75 KB
Newer Older
Tilman Steinweg's avatar
Tilman Steinweg committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
/*------------------------------------------------------------------------
 * Copyright (C) 2011 For the list of authors, see file AUTHORS.
 *
 * This file is part of SOFI2D.
 * 
 * SOFI2D 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.
 * 
 * SOFI2D 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 SOFI2D. See file COPYING and/or 
  * <http://www.gnu.org/licenses/gpl-2.0.html>.
--------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
 * Read elastic model properties (vp,vs,density) from files
 * This file contains function readmod, which has the purpose
 * to read data from model-files for viscoelastic simulation
 *
 * ----------------------------------------------------------------------*/

#include "fd.h"

#include "fd.h"

void readmod_visco(float  **  rho, float **  pi, float **  u,
float **  taus, float **  taup, float *  eta){

	extern float DT, *FL, TAU, TS;
	extern int NX, NY, NXG, NYG,  POS[3], L, MYID;
	extern char  MFILE[STRING_SIZE];	
	extern FILE *FP;

		
	/* local variables */
	float rhov, muv, piv, vp, vs, qp, qs;
	float *pts, ts, tp, sumu, sumpi, ws;
	int i, j, l, ii, jj;
	FILE *fp_vs, *fp_vp, *fp_rho, *fp_qp ,*fp_qs;
	char filename[STRING_SIZE];



	/* vector for maxwellbodies */
	pts=vector(1,L);
	for (l=1;l<=L;l++) {
		pts[l]=1.0/(2.0*PI*FL[l]);
		eta[l]=DT/pts[l];
	}

	ts=TAU;  
	tp=TAU;

	//previously : ws=2.0*PI*FL[1];
	ws=2.0*PI/TS;
	
	sumu=0.0; 
	sumpi=0.0;
	for (l=1;l<=L;l++){
		sumu=sumu+((ws*ws*pts[l]*pts[l]*ts)/(1.0+ws*ws*pts[l]*pts[l]));
		sumpi=sumpi+((ws*ws*pts[l]*pts[l]*tp)/(1.0+ws*ws*pts[l]*pts[l]));
	}

		

	   fprintf(FP,"\n...reading model information from model-files...\n");

	   fprintf(FP,"\t P-wave velocities:\n\t %s.vp\n\n",MFILE);
	   sprintf(filename,"%s.vp",MFILE);
	   fp_vp=fopen(filename,"r");
Tilman Steinweg's avatar
Tilman Steinweg committed
75
	   if ((fp_vp==NULL) && (MYID==0)) declare_error(" Could not open model file for P velocities ! ");
Tilman Steinweg's avatar
Tilman Steinweg committed
76
77
78
79
80


	   fprintf(FP,"\t Shear wave velocities:\n\t %s.vs\n\n",MFILE);
	   sprintf(filename,"%s.vs",MFILE);
	   fp_vs=fopen(filename,"r");
Tilman Steinweg's avatar
Tilman Steinweg committed
81
	   if ((fp_vs==NULL) && (MYID==0)) declare_error(" Could not open model file for shear velocities ! ");
Tilman Steinweg's avatar
Tilman Steinweg committed
82
83
84
85

	   fprintf(FP,"\t Density:\n\t %s.rho\n\n",MFILE);
	   sprintf(filename,"%s.rho",MFILE);
	   fp_rho=fopen(filename,"r");
Tilman Steinweg's avatar
Tilman Steinweg committed
86
	   if ((fp_rho==NULL) && (MYID==0)) declare_error(" Could not open model file for densities ! ");
Tilman Steinweg's avatar
Tilman Steinweg committed
87
88
89
90

	   fprintf(FP,"\t Qp:\n\t %s.qp\n\n",MFILE);
	   sprintf(filename,"%s.qp",MFILE);
	   fp_qp=fopen(filename,"r");
Tilman Steinweg's avatar
Tilman Steinweg committed
91
	   if ((fp_qp==NULL) && (MYID==0)) declare_error(" Could not open model file for Qp-values ! ");
Tilman Steinweg's avatar
Tilman Steinweg committed
92
93
94
95

	   fprintf(FP,"\t Qs:\n\t %s.qs\n\n",MFILE);
	   sprintf(filename,"%s.qs",MFILE);
	   fp_qs=fopen(filename,"r");
Tilman Steinweg's avatar
Tilman Steinweg committed
96
	   if ((fp_qs==NULL) && (MYID==0)) declare_error(" Could not open model file for Qs-values ! ");
Tilman Steinweg's avatar
Tilman Steinweg committed
97
98
99
100
101
102
103
104
105
106
107
108
	   

	/* loop over global grid */
		for (i=1;i<=NXG;i++){
			for (j=1;j<=NYG;j++){
			fread(&vp, sizeof(float), 1, fp_vp);
			fread(&vs, sizeof(float), 1, fp_vs);
			fread(&rhov, sizeof(float), 1, fp_rho);
			fread(&qp, sizeof(float), 1, fp_qp);
			fread(&qs, sizeof(float), 1, fp_qs);
				
		if ((isnan(vp)) && (MYID==0)) {
Tilman Steinweg's avatar
Tilman Steinweg committed
109
                	declare_error(" Found NaN-Values in Vp-Model !");}
Tilman Steinweg's avatar
Tilman Steinweg committed
110
111

                if ((isnan(vs)) && (MYID==0)) {
Tilman Steinweg's avatar
Tilman Steinweg committed
112
                	declare_error(" Found NaN-Values in Vs-Model !");}
Tilman Steinweg's avatar
Tilman Steinweg committed
113
114

                if ((isnan(rhov)) && (MYID==0)) {
Tilman Steinweg's avatar
Tilman Steinweg committed
115
                	declare_error(" Found NaN-Values in Rho-Model !");}
Tilman Steinweg's avatar
Tilman Steinweg committed
116
117

		if ((isnan(qp)) && (MYID==0)) {
Tilman Steinweg's avatar
Tilman Steinweg committed
118
                        declare_error(" Found NaN-Values in Vs-Model !");}
Tilman Steinweg's avatar
Tilman Steinweg committed
119
120

                if ((isnan(qs)) && (MYID==0)) {
Tilman Steinweg's avatar
Tilman Steinweg committed
121
                        declare_error(" Found NaN-Values in Rho-Model !");}
Tilman Steinweg's avatar
Tilman Steinweg committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170


			muv=vs*vs*rhov/(1.0+sumu);
			piv=vp*vp*rhov/(1.0+sumpi);

			/* 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;

					taus[jj][ii]=2.0/qs;
					taup[jj][ii]=2.0/qp;
					u[jj][ii]=muv;
					rho[jj][ii]=rhov;
					pi[jj][ii]=piv;
				}
			}
		}
	




	fclose(fp_vp);
	fclose(fp_vs);
	fclose(fp_rho);
	fclose(fp_qp);
	fclose(fp_qs);
	
	
	/* each PE writes his model to disk */
	   
	   
	sprintf(filename,"%s.sofi2D.rho",MFILE);

	writemod(filename,rho,3);

	MPI_Barrier(MPI_COMM_WORLD);

	if (MYID==0) mergemod(filename,3);

	free_vector(pts,1,L);
}