toy_example_ac_start.c 7.84 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
/*-----------------------------------------------------------------------------------------
 * Copyright (C) 2005  <name of author>
 *
 * This file is part of DENISE.
 * 
 * DENISE 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.
 * 
 * DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/

/* $Id: hh_el.c,v 2.3 2007/08/21 13:16:19 tbohlen Exp $*/
/*
 *   Model defined by flnode file. 
 */

#include "fd.h"
25
#include "../src/fd.h"
Tilman Steinweg's avatar
Tilman Steinweg committed
26 27 28 29

void model_acoustic(float **rho, float **pi){

	/* extern variables */
30
	extern int NX, NY, NXG, NYG, POS[3], MYID, FDORDER;
Tilman Steinweg's avatar
Tilman Steinweg committed
31
	extern float DH, DT;
32 33 34
	extern char MFILE[STRING_SIZE];
	extern char INV_MODELFILE[STRING_SIZE];
	extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
35
	extern int SWS_TAPER_FILE;
Tilman Steinweg's avatar
Tilman Steinweg committed
36

37 38 39
	/* local variables */
	float piv, vp, rhov, taperv, **taper;
	int i, j, ii, jj, l, nd;
Tilman Steinweg's avatar
Tilman Steinweg committed
40 41 42 43 44 45
	char modfile[STRING_SIZE];
	
	FILE *flfile;
	int nodes;
	char cline[256];
	
46
	float *fldepth, *flrho, *flvp, *fltaper;
Tilman Steinweg's avatar
Tilman Steinweg committed
47
	
48
	nd = FDORDER/2;
Tilman Steinweg's avatar
Tilman Steinweg committed
49 50
	
	/*read FL nodes from File*/
51 52
	nodes=5;
	
Tilman Steinweg's avatar
Tilman Steinweg committed
53 54
	fldepth=vector(1,nodes);
	flrho=vector(1,nodes);
55 56 57
	flvp=vector(1,nodes);
	fltaper=vector(1,nodes);
	
58
	if(SWS_TAPER_FILE) taper=matrix(-nd+1,NY+nd,-nd+1,NX+nd);
59
	
Tilman Steinweg's avatar
Tilman Steinweg committed
60 61 62
	flfile=fopen("model_true/flnodes.toy_example_ac.start","r");
	if (flfile==NULL) err(" FL-file could not be opened !");
	
63
	/* Read parameters */
Tilman Steinweg's avatar
Tilman Steinweg committed
64 65 66
	for (l=1;l<=nodes;l++){
		fgets(cline,255,flfile);
		if (cline[0]!='#'){
67
			sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &fltaper[l]);
Tilman Steinweg's avatar
Tilman Steinweg committed
68 69 70 71
		}
		else l=l-1;
	}
	
72
	/* Print parameters */
Tilman Steinweg's avatar
Tilman Steinweg committed
73
	if(MYID==0){
74 75 76 77 78 79 80 81
		printf(" ------------------------------------------------------------------ \n\n");
		printf(" Information of FL nodes: \n\n");
		printf(" \t depth \t vp \n\n");
		
		for (l=1;l<=nodes;l++){
			printf(" \t %f \t %f\n\n",fldepth[l],flvp[l]);
		}
		printf(" ------------------------------------------------------------------ \n\n");
Tilman Steinweg's avatar
Tilman Steinweg committed
82 83
	}
	
84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
	/* loop over global grid - vp */
	for (i=1;i<=NXG;i++){
		for (l=1;l<nodes;l++){
			if (fldepth[l]==fldepth[l+1]){
				if ((i==1) && (MYID==0)){
					printf("depth: %f m: double node\n",fldepth[l]);
				}
			}else{
				for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
					
					vp=0.0;
					
					vp=(DH*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
					vp=vp*1000.0;
					
					piv=vp;
					
					/* 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))){
Tilman Steinweg's avatar
Tilman Steinweg committed
104 105
						ii=i-POS[1]*NX;
						jj=j-POS[2]*NY;
106
						
Tilman Steinweg's avatar
Tilman Steinweg committed
107 108 109 110
						pi[jj][ii]=piv;
					}
				}
			}
111 112
		}
		for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
Tilman Steinweg's avatar
Tilman Steinweg committed
113
			
114 115 116 117
			vp=0.0;
			vp=flvp[nodes]*1000.0;
			
			piv=vp;
Tilman Steinweg's avatar
Tilman Steinweg committed
118

119 120 121 122 123
			/* 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;
Tilman Steinweg's avatar
Tilman Steinweg committed
124

125
				pi[jj][ii]=piv;
Tilman Steinweg's avatar
Tilman Steinweg committed
126 127
			}
		}
128
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
129
		
130
	/* loop over global grid - taper */
131
	if(SWS_TAPER_FILE){
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
		for (i=1;i<=NXG;i++){
			for (l=1;l<nodes;l++){
				if (fldepth[l]==fldepth[l+1]){
					if ((i==1) && (MYID==0)){
					printf("depth: %f m: double node\n",fldepth[l]);}}
				else{
					for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
						
						taperv=0.0;
						
						taperv=(DH*(j-1)-fldepth[l])*(fltaper[l+1]-fltaper[l])/(fldepth[l+1]-fldepth[l])+fltaper[l];
						
						/* 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;
						
						taper[jj][ii]=taperv;
						}
					}
				}
			}
			
			for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
158
				
159 160 161 162 163 164 165 166 167 168 169 170 171 172
				taperv=0.0;
				taperv=fltaper[nodes];
				
				/* 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;

						taper[jj][ii]=taperv;
						}
			}
		}
173
	}
174
		
Tilman Steinweg's avatar
Tilman Steinweg committed
175 176 177
	free_vector(fldepth,1,nodes);
	free_vector(flrho,1,nodes);
	free_vector(flvp,1,nodes);
178
	free_vector(fltaper,1,nodes);
Tilman Steinweg's avatar
Tilman Steinweg committed
179 180 181
	
	
	/**************************************************/
182
	/* creation of density models from true model */
Tilman Steinweg's avatar
Tilman Steinweg committed
183 184 185 186 187 188 189 190 191 192
	/**************************************************/
	
	/*read FL nodes from File*/
	nodes=6;
	fldepth=vector(1,nodes);
	flrho=vector(1,nodes);
	flvp=vector(1,nodes);	
	flfile=fopen("model_true/flnodes.toy_example_ac","r");
	if (flfile==NULL) err(" FL-file could not be opened !");
	
193
	/* Read parameters */
Tilman Steinweg's avatar
Tilman Steinweg committed
194 195 196 197 198 199 200 201
	for (l=1;l<=nodes;l++){
		fgets(cline,255,flfile);
		if (cline[0]!='#'){
			sscanf(cline,"%f%f%f",&fldepth[l], &flrho[l], &flvp[l]);
		}
		else l=l-1;
	}
	
202
	/* Print parameters */	
Tilman Steinweg's avatar
Tilman Steinweg committed
203
	if(MYID==0){
204 205 206 207 208 209 210 211
		printf(" ------------------------------------------------------------------ \n\n");
		printf(" Information of FL nodes: \n\n");
		printf(" \t depth \t vp \t rho \n\n");
		
		for (l=1;l<=nodes;l++){
			printf(" \t %f \t %f \t %f \n\n",fldepth[l],flvp[l],flrho[l]);
		}
		printf(" ------------------------------------------------------------------ \n\n");
Tilman Steinweg's avatar
Tilman Steinweg committed
212 213
	}
	
214
	/* loop over global grid - density */
215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232
	for (i=1;i<=NXG;i++){
		for (l=1;l<nodes;l++){
			if (fldepth[l]==fldepth[l+1]){
				if ((i==1) && (MYID==0)){
					printf("depth: %f m: double node\n",fldepth[l]);
				}
			}else{
				for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
		
					rhov=0.0;
				
					rhov=(DH*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
					rhov=rhov*1000.0;
					
					
					/* 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))){
Tilman Steinweg's avatar
Tilman Steinweg committed
233 234 235 236 237 238 239
						ii=i-POS[1]*NX;
						jj=j-POS[2]*NY;

						rho[jj][ii]=rhov;
					}
				}
			}
240 241 242
		}
		
		for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
Tilman Steinweg's avatar
Tilman Steinweg committed
243
			
244 245 246 247 248 249 250 251 252
			rhov=0.0;
			rhov=flrho[nodes]*1000.0;
			
			
			/* 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;
Tilman Steinweg's avatar
Tilman Steinweg committed
253

254
				rho[jj][ii]=rhov;
Tilman Steinweg's avatar
Tilman Steinweg committed
255 256
			}
		}
257
	}
Tilman Steinweg's avatar
Tilman Steinweg committed
258 259 260 261
	
	free_vector(fldepth,1,nodes);
	free_vector(flrho,1,nodes);
	free_vector(flvp,1,nodes);	
262
	
Tilman Steinweg's avatar
Tilman Steinweg committed
263 264 265 266 267 268 269 270
	
        sprintf(modfile,"%s_rho_it0.bin",INV_MODELFILE);
	writemod(modfile,rho,3);
	MPI_Barrier(MPI_COMM_WORLD);
	if (MYID==0) mergemod(modfile,3);
	MPI_Barrier(MPI_COMM_WORLD); 
	sprintf(modfile,"%s_rho_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
	remove(modfile);
271
	
Tilman Steinweg's avatar
Tilman Steinweg committed
272 273 274 275 276 277 278
	sprintf(modfile,"%s_vp_it0.bin",INV_MODELFILE);
	writemod(modfile,pi,3);
	MPI_Barrier(MPI_COMM_WORLD);
	if (MYID==0) mergemod(modfile,3);
	MPI_Barrier(MPI_COMM_WORLD); 
	sprintf(modfile,"%s_vp_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
	remove(modfile);
279
	
280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296
	if(SWS_TAPER_FILE){
		sprintf(modfile,"%s",TAPER_FILE_NAME);
		writemod(modfile,taper,3);
		MPI_Barrier(MPI_COMM_WORLD);
		if (MYID==0) mergemod(modfile,3);
		MPI_Barrier(MPI_COMM_WORLD); 
		sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME,POS[1],POS[2]);
		remove(modfile);
		
		sprintf(modfile,"%s",TAPER_FILE_NAME_RHO);
		writemod(modfile,taper,3);
		MPI_Barrier(MPI_COMM_WORLD);
		if (MYID==0) mergemod(modfile,3);
		MPI_Barrier(MPI_COMM_WORLD); 
		sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME_RHO,POS[1],POS[2]);
		remove(modfile);
	}
297
	
298
	if(SWS_TAPER_FILE) free_matrix(taper,-nd+1,NY+nd,-nd+1,NX+nd);
299
	
Tilman Steinweg's avatar
Tilman Steinweg committed
300
}