suea2df.c

su 的源代码库

		ttr[j-1]=ttr[j]+dz*ca/vr[j];
	}

	warn("vmin=%f vmax=%f",vlim[0],vlim[1]);
	warn("model set");
	
	/* calculate stability and dispersion */
	warn("for stablity dt should be < %f: x-dir and %f: z-dir",
		0.6*dx/vlim[1],0.6*dz/vlim[1]);
	warn("for low dispersion fmax should be < %f ",
		vlim[0]/(5*dx));
	
		
	/* get source function */
	ns=NINT(ts/dt)+1;
	source = alloc1float(ns);
	i=get_source(dt,ts,favg,wtype,source);
	warn("source set");
	if (*sofile!='\0') {
		fprintf(soeisfp,"%s %s %d %d %f %f\n",wtype,stype,ns,i,favg,ts);
		for (i=0; i<ns; i++)
			fprintf(soeisfp,"%f %f\n",i*dt,source[i]);
 	}

	pfac=1.;
	for (i=0; i<ns; i++) source[i]*=pfac;

 	warn("%s %s %d %d %f %f",wtype,stype,ns,i,favg,ts);
	
	/* allocate space for variables*/
	/*  CJ: ADJUSTMENT MADE 030502 - txx,tzz and w extended for symmetry */
	nbytes=ngrids*nxpadded*nzpadded*FSIZE;
	u = alloc2float(nxpadded,nzpadded);
	nbytes=nbytes+nxpadded*nzpadded*FSIZE;
	w = alloc2float(nxpadded+1,nzpadded);
	nbytes=nbytes+(nxpadded+1)*nzpadded*FSIZE;
	txx = alloc2float(nxpadded+1,nzpadded);
	nbytes=nbytes+(nxpadded+1)*nzpadded*FSIZE;
	tzz = alloc2float(nxpadded+1,nzpadded);
	nbytes=nbytes+(nxpadded+1)*nzpadded*FSIZE;
	txz = alloc2float(nxpadded,nzpadded);
	nbytes=nbytes+nxpadded*nzpadded*FSIZE;
	ngrids+=5;
	
	warn("%f MBytes allocated for model + variables\n",
			0.000001*nbytes);

	/* scale Q-factor */
	if(qsw==1) {
		for (i=0; i<nxpadded; i++){
			for (j=0; j<nzpadded; j++) {
				q[j][i]=exp(-PI*dt*favg/q[j][i]); 
			}
		}	
	}
	
	/* zero arrays */
	for (i=0; i<nxpadded; i++){
		for (j=0; j<nzpadded; j++) {
			u[j][i]=0;
			w[j][i]=0;
			txx[j][i]=0;
			tzz[j][i]=0;
			txz[j][i]=0;
		}
	}

	/* write character headers */
	strcpy(tfile,hsfile);
	hcfile=strcat(tfile,".chd");
	if((hchdfp=fopen(hcfile,"w+"))==NULL)
		err("cannot open hcfile=%s\n",hcfile);
	write_chd(nxpadded,nzpadded,nt,dx,dz,dt,fx,xmax,fz,zmax,
			sx,sz,favg,ts,wtype,stype,
			bc,qsw,aniso,hsz,hchdfp,1,mfile);
	fclose(hchdfp);

	strcpy(tfile,vrslfile);vrshfile=strcat(tfile,".chd");
	if((vchdfp=fopen(vrshfile,"w+"))==NULL)
		err("cannot open vrshfile=%s\n",vrshfile);
	write_chd(nxpadded,nzpadded,nt,dx,dz,dt,fx,xmax,fz,zmax,
			sx,sz,favg,ts,wtype,stype,
			bc,qsw,aniso,vsx,vchdfp,2,mfile);
	fclose(vchdfp);

	/* begin wave propagation */
	k=0;
	for (t=ft; t<=lt; t=t+dt) {

		/* calculate area for active grid */
		calc_area(vlim[1],dtx,dtz,limits,is,js,k,nxpadded,nzpadded);

		ifx=limits[0];
		ilx=limits[1];
		jfz=limits[2];
		jlz=limits[3];

		if (strcmp(stype,"pw")==0){
			ifx=2;
			ilx=nxpadded-3;
			jfz=2;
			jlz=nzpadded-3;
		}

		/* introduce velocity source */
		if (k<ns) {
			if (strcmp(stype,"v")==0)
				add_v_source(u,w,dt*source[k],is,js);
		}

		/* update velocities */
		if(tsw==0)
			update_vel(ifx,ilx,jfz,jlz,u,w,txx,tzz,txz,
					rho,c55,dtx,dtz);
		if(tsw==1)
			update_vel_tsw(ifx,ilx,jfz,jlz,u,w,txx,tzz,
					txz,rho,c55,dtx,dtz);
	
		/* introduce plane wave velocity source */
		if (strcmp(stype,"pw")==0){
		  		add_pw_source_V(u,w,txx,tzz,txz,c11,c55,
						source,sang,wbc[1],
						nxpadded-wbc[3]-1,
						wbc[0],nzpadded-wbc[2]-1,
						ns,k,dx,dz,dt,vl,vb,vr,
						ttl,ttb,ttr);
		}

		/* apply free surface bc for velocities */
		if(bc[0]==2) fs4v_bc_top(u,w,c11,c55,nxpadded,nzpadded);
	
		/* apply symmetry bc for velocities */
		if(bc[0]==1) sym4v_bc_top(u,w,nxpadded,nzpadded);
		if(bc[1]==1) sym4v_bc_left(u,w,nxpadded,nzpadded);
		if(bc[2]==1) sym4v_bc_bot(u,w,nxpadded,nzpadded);
		if(bc[3]==1) sym4v_bc_right(u,w,nxpadded,nzpadded);
	
		/* introduce stress source */
		if (k<ns) {
	  		if (strcmp(stype,"p")==0)
				add_p_source(txx,tzz,dt*source[k],is,js);
		}
	
		/* update stresses */
		if(!aniso) {  /* isotropic */
			update_stress_iso(ifx,ilx,jfz,jlz,u,w,txx,tzz,txz,
						c11,c55,dtx,dtz);
		} else {
			update_stress_ani(ifx,ilx,jfz,jlz,u,w,txx,tzz,txz,
						c11,c55,c33,c13,c15,c35,
						dtx,dtz);
		}

		/* introduce plane wave stress source */
		if (strcmp(stype,"pw")==0) {
			add_pw_source_S(u,w,txx,tzz,txz,c11,c55,
						source,sang,wbc[1],
						nxpadded-wbc[3]-1,wbc[0],
						nzpadded-wbc[2]-1,ns,k,
						dx,dz,dt,vl,vb,vr,ttl,
						ttb,ttr);
		}

		/* apply free surface bc for stresses */
		if(bc[0]==2) fs4s_bc_top(txx,tzz,txz,nxpadded,nzpadded);

		/* apply symmetry bc for stresses */
		if(bc[0]==1) sym4s_bc_top(txx,tzz,txz,nxpadded,nzpadded);
		if(bc[1]==1) sym4s_bc_left(txx,tzz,txz,nxpadded,nzpadded);
		if(bc[2]==1) sym4s_bc_bot(txx,tzz,txz,nxpadded,nzpadded);
		if(bc[3]==1) sym4s_bc_right(txx,tzz,txz,nxpadded,nzpadded);

		/* apply absorbing bc */
		if(bc[0]>2)
			abs_bc_top(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc0,bc,wbc,ifx,ilx);
		if(bc[1]>2)
			abs_bc_left(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc1,bc,wbc,jfz,jlz);
		if(bc[2]>2)
			abs_bc_bot(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc2,bc,wbc,ifx,ilx);
		if(bc[3]>2)
			abs_bc_right(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc3,bc,wbc,jfz,jlz);

		/* attenuation */
		if(qsw==1) {
			for (j=jfz-2; j<jlz+2; j++){
				for (i=ifx-2; i<ilx+2; i++) {
			 		u[j][i]*=q[j][i];
					w[j][i]*=q[j][i];
			 		txx[j][i]*=q[j][i];
					tzz[j][i]*=q[j][i];
					txz[j][i]*=q[j][i];
				}
			}	
		}
	
		/* "energy" */
		energy=0;
		for (i=ifx; i<ilx; i++){
			for (j=jfz; j<jlz; j++) {
				energy=energy+u[j][i]*u[j][i]+w[j][i]*w[j][i]; 
			}
		}
	
		/* store wave motion on file (seismograms) */
		for (i=0 ; i < nxpadded ; ++i)
			fwrite(&u[hs1][i], 4, 1, hseisfp);
		for (i=0 ; i < nxpadded ; ++i) 
			fwrite(&w[hs1][i], 4, 1, hseisfp);
		for (j=0 ; j < nzpadded ; ++j) 
			fwrite(&u[j][vs2], 4, 1, vseisfp);
		for (j=0 ; j < nzpadded ; ++j) 
			fwrite(&w[j][vs2], 4, 1, vseisfp);
				
		/* make snapshot */
		for (i=0; i<nsnap; i++) {
			if (!mode) {
				if(k==isnap[i] && *snfile!='\0') {
					make_snap(u,w,sx,sz,dx,dz,
							nxpadded,nzpadded,
							nx,nz,
							k,t,fx,fz,tr,
							sneisfp,wbc);
		 		} 
			} else {
				if(k==isnap[i] && *snfile!='\0') {
		  			make_stress_snap(txx,tzz,txz,sx, sz,
								dx, dz,
								nxpadded, 
								nzpadded,
								k,t,fx,fz,tr,
								sneisfp,wbc);
				}
		 	}
	 	}

		if((verbose==1) && (k%50==0))
			fprintf(stderr,"%d %f %f %d %d %d %d \n",
 				k,t,energy*efac,ifx,ilx,jfz,jlz);


		k=k+1;
	}

	warn("propagation completed\n");

	/* free space */
	free2float(u);
	free2float(w);
	free2float(txx);
	free2float(tzz);
	free2float(txz);
	warn("propagation space freed\n");
	
	/* read in wave motion from file and make seismograms */
	rewind(hseisfp);	/*horizontal line*/ 

	ut = alloc2float(nxpadded,nt); 
	wt = alloc2float(nxpadded,nt);

	for (k=0 ; k < nt ; ++k){
		for (i=0 ; i < nxpadded ; ++i) 
			fread(&ut[k][i], 4, 1, hseisfp);
		for (i=0 ; i < nxpadded ; ++i)
			fread(&wt[k][i], 4, 1, hseisfp);
	}
	fclose(hseisfp);
	hseisfp=fopen(hsfile,"w");

	make_seis(ut,wt,sx,sz,dx,dt,nxpadded,nt,hsz,fx,fz,
			trh,hseisfp,wbc[1],wbc[3],1,verbose);
	fclose(hseisfp);
	free2float(ut);
	free2float(wt);
	warn("u time section output\n");
	
	rewind(vseisfp);	/*vertical line*/
	ut = alloc2float(nzpadded,nt); wt = alloc2float(nzpadded,nt);
	for (k=0 ; k < nt ; ++k){
	  for (j=0 ; j < nzpadded ; ++j) fread(&ut[k][j], 4, 1, vseisfp);
	  for (j=0 ; j < nzpadded ; ++j) fread(&wt[k][j], 4, 1, vseisfp);
	}
	fclose(vseisfp); vseisfp=fopen(vrslfile,"w");

	make_seis(ut,wt,sx,sz,dz,dt,nzpadded,nt,vsx,fx,fz,
			trv,vseisfp,wbc[0],wbc[2],2,verbose);

	fclose(vseisfp); free2float(ut); free2float(wt);
	warn("w time section output\n");

	fclose(sneisfp);

	return(CWP_Exit());
}

int get_source(float dt, float ts, float favg, char *wtype, float *source)
{
	int i,ns; float t,x,xx,pi2;
	i=0; t=0; pi2=PI*PI;
	ns=NINT(ts/dt)+1;
	while (i<ns) {
	  if(strcmp(wtype,"ri")==0){ /*Ricker*/
		x=favg*(t-ts/2);xx=x*x;
		source[i]=(1-2*pi2*(xx))*exp(-pi2*xx);
	  } 							  
	  if(strcmp(wtype,"dg")==0){ /*Derivative of Gaussian*/
		x=-4*favg*favg*pi2/log(0.1);
		source[i]=-2*pi2*(t-ts/2)*exp(-x*(t-ts/2)*(t-ts/2));
	  } 							  
	  if(strcmp(wtype,"ga")==0){ /*Gaussian favg is 10% level for fmax*/
		x=-favg*favg*pi2/log(0.1);
		source[i]=exp(-x*(t-ts/2)*(t-ts/2));
	  } 							  
	  if(strcmp(wtype,"sp")==0){ /*spike*/
		source[0]=1;
	  } 						
	  if(strcmp(wtype,"sp2")==0){ /*spike*/
		source[0]=1;
		source[1]=1;
	  } 						
 	  /*warn("%d %f %f\n",i,t,source[i]);*/
	  t=t+dt; i=i+1;
	}
	return(i);
}

void update_vel(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz)
{
	int i,j;
	float dtxx=0,dtzz=0,dtxz=0,dtzx=0;
	 	for (j=jfz; j<=jlz; ++j) {
	 	  for (i=ifx; i<=ilx; ++i) {
	 		dtxx = F1*(txx[j][i+1]-txx[j][i])
		 +F2*(txx[j][i+2]-txx[j][i-1]);
		dtxz = F1*(txz[j+1][i]-txz[j][i])
	 		+F2*(txz[j+2][i]-txz[j-1][i]);
			u[j][i] = u[j][i] + ( dtx*dtxx + dtz*dtxz ) * rho[j][i];
	 	  }
	 	  for (i=ifx; i<=ilx+1; ++i) {
			dtzz = F1*(tzz[j][i]-tzz[j-1][i])
	 		+F2*(tzz[j+1][i]-tzz[j-2][i]);
	 		dtzx = F1*(txz[j][i]-txz[j][i-1])
	 		+F2*(txz[j][i+1]-txz[j][i-2]);
	 		w[j][i] = w[j][i] + ( dtx*dtzx + dtz*dtzz ) * rho[j][i];
	 	  }
	 	}
}
void update_vel_tsw(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz)
{
	int i,j;
	float dtxx=0,dtzz=0,dtxz=0,dtzx=0;
	 	for (j=jfz; j<=jlz; ++j) {
	 	  for (i=ifx; i<=ilx; ++i) {
	 		dtxx = F1*(txx[j][i+1]-txx[j][i])
		 +F2*(txx[j][i+2]-txx[j][i-1]);
/*  		Only use 4th order shear stress derivatives if material is non-fluid, otherwise problems occur at grazing angles near fluid-solid boundary */
	 		if(c55[j][i]>0) { 
		 dtxz = F1*(txz[j+1][i]-txz[j][i])
	 		+F2*(txz[j+2][i]-txz[j-1][i]);
		}
		else{
		 dtxz = (txz[j+1][i]-txz[j][i]);
		}
	 		u[j][i] = u[j][i] + ( dtx*dtxx + dtz*dtxz ) * rho[j][i];
	 	  }
	 	  for (i=ifx; i<=ilx+1; ++i) {
	 		dtzz = F1*(tzz[j][i]-tzz[j-1][i])
	 		+F2*(tzz[j+1][i]-tzz[j-2][i]);
/*  		Only use 4th order shear stress derivatives if material is non-fluid, otherwise problems occur at grazing angles near fluid-solid boundary */
	 		if(c55[j][i]>0) { 
	 		dtzx = F1*(txz[j][i]-txz[j][i-1])
	 		+F2*(txz[j][i+1]-txz[j][i-2]);
		}
		else{
	 		dtzx = (txz[j][i]-txz[j][i-1]);
		}
	 		w[j][i] = w[j][i] + ( dtx*dtzx + dtz*dtzz ) * rho[j][i];
	 	  }
	 	}
}


void update_stress_iso(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float dtx, float dtz)
{
	int i,j;
	float dux,duz,dwz,dwx;
	 	for (j=jfz; j<=jlz; ++j) {
	 	  for (i=ifx; i<=ilx+1; ++i) {
	 		dwz = F1*(w[j+1][i]-w[j][i])
	 		+F2*(w[j+2][i]-w[j-1][i]);
	 		dux = F1*(u[j][i]-u[j][i-1])
	 		+F2*(u[j][i+1]-u[j][i-2]);
	 		txx[j][i] = txx[j][i] +
	 		c11[j][i]*dtx*dux + (c11[j][i]-2*c55[j][i])*dtz*dwz;
	 		tzz[j][i] = tzz[j][i] +
	 		c11[j][i]*dtz*dwz + (c11[j][i]-2*c55[j][i])*dtx*dux;
	 	  }
	 	  for (i=ifx; i<=ilx; ++i) {
	 		dwx = F1*(w[j][i+1]-w[j][i])
		 +F2*(w[j][i+2]-w[j][i-1]);
	 		duz = F1*(u[j][i]-u[j-1][i])
	 		+F2*(u[j+1][i]-u[j-2][i]);
	 		txz[j][i] = txz[j][i] +
	 		c55[j][i]*(dtz*duz + dtx*dwx);
	 	  }
	 	}
}

void update_stress_ani(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float **c33, float **c13, 
	float **c15, float **c35, float dtx, float dtz)
{
	int i,j;
	float dux,duz,dwz,dwx;

 	for (j=jfz; j<=jlz; ++j) {
		for (i=ifx; i<=ilx; ++i) {

			dwx = F1*(w[j][i+1]-w[j][i])
				+F2*(w[j][i+2]-w[j][i-1]);

			duz = F1*(u[j][i]-u[j-1][i])
	 					+F2*(u[j+1][i]-u[j-2][i]);

	 				dwz = F1*(w[j+1][i]-w[j][i])
	 					+F2*(w[j+2][i]-w[j-1][i]);

					dux = F1*(u[j][i]-u[j][i-1])
			 		+F2*(u[j][i+1]-u[j][i-2]);

			txx[j][i] = txx[j][i] 
					+ c11[j][i]*dtx*dux 
					+ c13[j][i]*dtz*dwz 
					+ c15[j][i]*(dtz*duz + dtx*dwx);

					tzz[j][i] = tzz[j][i] 
					+ c33[j][i]*dtz*dwz 
					+ c13[j][i]*dtx*dux 
					+ c35[j][i]*(dtz*duz + dtx*dwx); 

			txz[j][i] = txz[j][i] 
					+ c15[j][i]*dtx*dux 
					+ c35[j][i]*dtz*dwz 
					+ c55[j][i]*(dtz*duz + dtx*dwx);
	 	  	}
		i=ilx+1;

		dwx = F1*(w[j][i+1]-w[j][i])