sufctanismod.c

su 的源代码库

*************************************************************************
* Input: 
*	int nx		number of grids in x direction 
*	int nz		number of grids in z direction
*	float **cc	elastic parameter profile
*	float **rho	density reverse
*
*	int mbx1		left boundary position
*	int mby1		top boundary position
*	int mbx2		right boundary postion 
*	int mby2		bottom boundary position
*	
*************************************************************************
* Output: 
*	float *vell	velocity at the left boundary 
*	float *velr	velocity at the right boundary 
*	float *velt	velocity at the top boundary 
*	float *velb	velocity at the bottom boundary 
*
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	boundary_vel(float **cc, float **rho, 
		     float *vell, float *velr, float *velt, float *velb, 
		     int nx, int nz, int mbx1, int mbz1, int mbx2, int mbz2)
{
  int 	i, k;

  /* velocity at side boundary */
  /* note:  here rho = 1/rho  */
  for (k=0; k<nz; ++k) {
    vell[k] = sqrt(cc[mbx1][k]*rho[mbx1][k]);
    velr[k] = sqrt(cc[mbx2-1][k]*rho[mbx2-1][k]);
  }

  /* velocity at top and bottom boundary */
  for (i=0; i<nx; ++i) {
    velt[i] = sqrt(cc[i][mbz1]*rho[i][mbz1]);
    velb[i] = sqrt(cc[i][mbz2-1]*rho[i][mbz2-1]);
  }
}




/********************************************************************
* fct2d1o - 2-D FCT correction for the FIRST-ORDER wave equation:
********************************************************************
* Note: second order accuracy in space
********************************************************************
* Input:
*
*	for the 2nd-order wave equation
*	float **u	solution at previous time level
*	float **u1	solution at next time level
*
*	float dx	spacing in x-direction
*	float dz	spacing in z-direction
*	int	mbx1	sample start in x-direction
*	int	mbx2	sample end in x-direction
*	int	mbz1	sample start in z-direction
*	int	mbz2	sample end in z-direction
*	float eta0	diffusion coefficient
*	float deta0dx	gradient of eta0 in x-direction
*	float deta0dz	gradient of eta0 in z-direction
*	float eta	anti-diffusion coefficient
*	float detadx	gradient of eta in x-direction
*	float detadz	gradient of eta in z-direction
*
********************************************************************
* Output:
*
*	for 2nd-order wave equation
*	float **u1	solution after the FCT correction
*
*	Working Arrays 
*	f0[][]	
*	f1[][]
*	 
********************************************************************
* References:
*	The detailed algorithm description is given in the article
*	"Elimination of dispersion in finite-difference modeling 
*	and migration"	in CWP-137, project review, page 155-174.
*
*	Original reference to FCT method:
*	Boris, J., and Book, D., 1973, Flux-corrected transport. I.
*	SHASTA, a fluid transport algorithm that works: 
*	Journal of Computational Physics, vol. 11, p. 38-69.
*
********************************************************************
* Author: Tong Fei, Colorado School of Mines, 1993.
********************************************************************/
void	fct2d1o(float **u, float **u1, int nx, int nz, 
		float eta0, float eta, float deta0dx, float deta0dz, 
		float detadx, float detadz, float dx, float dz, 
		int mbx1, int mbx2, int mbz1, int mbz2, 
		float **f0, float **f1)
{
  int	ix, iz;
  float	sf, eta00, eta11, min;
  float	deta0dxdx=deta0dx*dx, deta0dzdz=deta0dz*dz, 
    detadxdx=detadx*dx, detadzdz=detadz*dz;
	
  /*  diffusive fluxes from first differences of u[][] and u1[][] */
  for (ix=mbx1; ix<mbx2-1; ix++)	
    for (iz=mbz1; iz<mbz2-1; iz++)
      {
	eta00=eta0*(1.0+deta0dxdx*(ix-1)+deta0dzdz*(iz-1));
	eta11=eta*(1.0+detadxdx*(ix-1)+detadzdz*(iz-1));
	f0[ix][iz]=eta00*(u[ix+1][iz]-u[ix][iz]);
	f1[ix][iz]=eta11*(u1[ix+1][iz]-u1[ix][iz]);
      }
  for (ix=mbx1; ix<mbx2-1; ix++)	
    {
      eta00=eta0*(1.0+deta0dxdx*(ix-1)+deta0dzdz*(mbz2-1));
      eta11=eta*(1.0+detadxdx*(ix-1)+detadzdz*(mbz2-1));
      f0[ix][mbz2-1]=eta00*(u[ix+1][mbz2-1]-u[ix][mbz2-1]);
      f1[ix][mbz2-1]=eta11*(u1[ix+1][mbz2-1]-u1[ix][mbz2-1]);
    }

  /*  diffusive fluxes from first differences of u1[][]  */

  /* diffuse the first transported u[][] using saved first differences */
  for (ix=mbx1+1; ix<mbx2-1; ix++)
    for (iz=mbz1+1; iz<mbz2-1; iz++)
      u1[ix][iz]=u1[ix][iz]+(f0[ix][iz]-f0[ix-1][iz]);

  /*  take first differences of the diffused u[][]  */
  for (ix=mbx1; ix<mbx2-1; ix++)
    for (iz=mbz1; iz<mbz2-1; iz++)
      {
	f0[ix][iz]=u1[ix+1][iz]-u1[ix][iz];
      }
  for (ix=mbx1; ix<mbx2-1; ix++)
    f0[ix][mbz2-1]=u1[ix+1][mbz2-1]-u1[ix][mbz2-1];
			

  /*  limit the antidiffusive fluxes  */

  for (ix=mbx1+1; ix<mbx2-2; ix++)
    for (iz=mbz1; iz<mbz2; iz++)
      {
	/*		sf=SIGN(f1[ix][iz]);
			f1[ix][iz]=sf*MAX(0.0, MIN(MIN(sf*f0[ix-1][iz], 
			sf*f1[ix][iz]), sf*f0[ix+1][iz]));
	*/
	if (f1[ix][iz] > 0) {
	  sf=1.0;
	} else {
	  sf=-1.0;
	}

	if (sf*f0[ix-1][iz] > sf*f1[ix][iz]) {
	  min = sf*f1[ix][iz];
	} else {
	  min = sf*f0[ix-1][iz];
	}

	if (min > sf*f0[ix+1][iz]) {
	  min = sf*f0[ix+1][iz];
	} 
			
	if (min > 0.0) {
	  f1[ix][iz]=sf*min;
	} else {
	  f1[ix][iz]=0.0;
	}
      }
		  
  for (iz=mbz1; iz<mbz2; iz++)
    {
      /*	sf=SIGN(f1[mbx1][iz]);
		f1[mbx1][iz]=sf*MAX(0.0, MIN(sf*f1[mbx1][iz], 
		sf*f0[mbx1+1][iz])); 
		sf=SIGN(f1[mbx2-2][iz]);
		f1[mbx2-2][iz]=sf*MAX(0.0, MIN(sf*f0[mbx2-3][iz], 
		sf*f1[mbx2-2][iz]));  
      */
      if (f1[mbx1][iz] > 0.0) {
	sf=1.0;
      } else {
	sf=-1.0;
      }

      if (sf*f1[mbx1][iz] > sf*f0[mbx1+1][iz]) {
	min = sf*f0[mbx1+1][iz];
      } else {
	min = sf*f1[mbx1][iz];
      }

      if (min > 0.0) {
	f1[mbx1][iz]=sf*min;
      } else {
	f1[mbx1][iz]=0.0;
      }

      if (f1[mbx2-2][iz] > 0.0) {
	sf=1.0;
      } else {
	sf=-1.0;
      }

      if (sf*f0[mbx2-3][iz] > sf*f1[mbx2-2][iz]) {
	min = sf*f1[mbx2-2][iz];
      } else {
	min = sf*f0[mbx2-3][iz];
      } 

      if (min > 0.0) {
	f1[mbx2-2][iz]=sf*min;
      } else {
	f1[mbx2-2][iz]=0.0;
      }
    }
	
			  
  /*  antidiffuse with the limited fluxes	  */
  for (ix=mbx1+1; ix<mbx2-1; ix++)
    for (iz=mbz1+1; iz<mbz2-1; iz++)
      u1[ix][iz]=u1[ix][iz]-(f1[ix][iz]-f1[ix-1][iz]);

  /*  diffusive fluxes from first differences of u[][] and u1[][] */
  for (ix=mbx1; ix<mbx2-1; ix++)	
    for (iz=mbz1; iz<mbz2-1; iz++)
      {
	eta00=eta0*(1.0+deta0dxdx*(ix-1)+deta0dzdz*(iz-1));
	eta11=eta*(1.0+detadxdx*(ix-1)+detadzdz*(iz-1));
	f0[ix][iz]=eta00*(u[ix][iz+1]-u[ix][iz]);
	f1[ix][iz]=eta11*(u1[ix][iz+1]-u1[ix][iz]);
      }
  for (iz=mbz1; iz<mbz2-1; iz++)
    {	
      eta00=eta0*(1.0+deta0dxdx*(mbx2-1)+deta0dzdz*(mbz2-1));
      eta11=eta*(1.0+detadxdx*(mbx2-1)+detadzdz*(mbz2-1));
      f0[mbx2-1][iz]=eta00*(u[mbx2-1][iz+1]-u[mbx2-1][iz]);
      f1[mbx2-1][iz]=eta11*(u1[mbx2-1][iz+1]-u1[mbx2-1][iz]);
    }	

  /*  diffusive fluxes from first differences of u1[][]  */

  /* diffuse the first transported u[][] using saved first differences */
  for (ix=mbx1+1; ix<mbx2-1; ix++)
    for (iz=mbz1+1; iz<mbz2-1; iz++)
      u1[ix][iz]=u1[ix][iz]+(f0[ix][iz]-f0[ix][iz-1]);

  /*  take first differences of the diffused u[][]  */
  for (ix=mbx1; ix<mbx2-1; ix++)
    for (iz=mbz1; iz<mbz2-1; iz++)
      {
	f0[ix][iz]=u1[ix][iz+1]-u1[ix][iz]; 
      }
  for (iz=mbx1; iz<mbz2-1; iz++)
    f0[mbx2-1][iz]=u1[mbx2-1][iz+1]-u1[mbx2-1][iz]; 
			

  /*  limit the antidiffusive fluxes  */

	
  for (ix=mbx1; ix<mbx2; ix++)
    for (iz=mbz1+1; iz<mbz2-2; iz++)
      {
	/*		sf=SIGN(f1[ix][iz]);
			f1[ix][iz]=sf*MAX(0.0, MIN(MIN(sf*f0[ix][iz-1], 
			sf*f1[ix][iz]), sf*f0[ix][iz+1]));
	*/
	if (f1[ix][iz] > 0) {
	  sf=1.0;
	} else {
	  sf=-1.0;
	}

	if (sf*f0[ix][iz-1] > sf*f1[ix][iz]) {
	  min = sf*f1[ix][iz];
	} else {
	  min = sf*f0[ix][iz-1];
	}

	if (min > sf*f0[ix][iz+1]) {
	  min = sf*f0[ix][iz+1];
	} 
			
	if (min > 0.0) {
	  f1[ix][iz]=sf*min;
	} else {
	  f1[ix][iz]=0.0;
	}

      }
		
  for (ix=mbx1; ix<mbx2; ix++)
    {
      /*	sf=SIGN(f1[ix][mbz1]);
		f1[ix][mbz1]=sf*MAX(0.0, MIN(sf*f1[ix][mbz1], 
		sf*f0[ix][mbz1+1])); 
		sf=SIGN(f1[ix][mbz2-2]);
		f1[ix][mbz2-2]=sf*MAX(0.0, MIN(sf*f0[ix][mbz2-3], 
		sf*f1[ix][mbz2-2])); 
      */
      if (f1[ix][mbz1] > 0.0) {
	sf=1.0;
      } else {
	sf=-1.0;
      }

      if (sf*f1[ix][mbz1] > sf*f0[ix][mbz1+1]) {
	min = sf*f0[ix][mbz1+1];
      } else {
	min = sf*f1[ix][mbz1];
      }

      if (min > 0.0) {
	f1[ix][mbz1]=sf*min;
      } else {
	f1[ix][mbz1]=0.0;
      }

      if (f1[ix][mbz2-2] > 0.0) {
	sf=1.0;
      } else {
	sf=-1.0;
      }

      if (sf*f0[ix][mbz2-3] > sf*f1[ix][mbz2-2]) {
	min = sf*f1[ix][mbz2-2];
      } else {
	min = sf*f0[ix][mbz2-3];
      } 

      if (min > 0.0) {
	f1[ix][mbz2-2]=sf*min;
      } else {
	f1[ix][mbz2-2]=0.0;
      }
 
    }
			  
  /*  antidiffuse with the limited fluxes	  */
  for (ix=mbx1+1; ix<mbx2-1; ix++)
    for (iz=mbz1+1; iz<mbz2-1; iz++)
      u1[ix][iz]=u1[ix][iz]-(f1[ix][iz]-f1[ix][iz-1]);

}

/************************************************************************
* tforce_ricker -- Compute the	time response of a source function as
*	a Ricker wavelet with peak frequency "fpeak" Hz.	
*************************************************************************
*************************************************************************
* Input: 
*	int nt		number of time step
*	float dt 	time step
*	float fpeak	peak frequency of the Ricker wavelet
*	
*************************************************************************
* Output: float *tforce		source array
*
*************************************************************************
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void tforce_ricker(int nt, float *tforce, float dt, float fpeak)
{
  int	it;
  float	t1, t0; 

  t0=1.0/fpeak;

  for (it=0; it<nt; it++) {
    t1=it*dt;
    tforce[it] = exp(-PI*PI*fpeak*fpeak*(t1-t0)*(t1-t0))*
      (1.0-2.*PI*PI*fpeak*fpeak*(t1-t0)*(t1-t0));
  }
}

/************************************************************************
* tforce_akb -- Compute the	time response of a source function as
*	a wavelet based on a wavelet used by Alford, Kelly, and Boore.
*************************************************************************
*************************************************************************
* Input: 
*	int nt		number of time step
*	float dt 	time step
*	float fpeak	peak frequency of the wavelet
*	
*************************************************************************
* Output: float *tforce		source array
*
*************************************************************************
* Reference:
*	Alford, R., Kelly, K., and Boore, D., 1974,
*	Accuracy of finite-difference modeling of the acoustic wave