sufctanismod.c

su 的源代码库

		     float **aa, float **cc, float **ff, float **ll, float **nn, 
		     float **rho, float **xzsource, float *fx, float *fy, float *fz, 
		     float dx, float dz, float dt, int  nx, int  nz, 
		     int dofct, int isurf, float eta0, float eta, 
		     float deta0dx, float deta0dz, float detadx, float detadz, 
		     int mbx1, int mbx2, int mbz1, int mbz2)
{
  int 	ix, iz;
  float	dtdx=dt/dx, dtdz=dt/dz, rdxdz=1.0/dx/dz;

  float	**u1, **f0, **f1;

  /*   allocate temporery space	*/
  u1 = alloc2float(nz, nx);
  f0 = alloc2float(nz, nx);
  f1 = alloc2float(nz, nx);


  /*  solve u at next time step  */
  difference(mbx1, mbx2, mbz1, mbz2, it, u1, u, 
	     xzsource, fx, dt, rdxdz); 
  /*  finite-difference to compute the contribution of e11 to u */
  difference_2x(mbx1, mbx2, mbz1, mbz2, -1, dtdx, u1, e11, aa); 
  /*  finite-difference to compute the contribution of e33 to u */
  difference_2x(mbx1, mbx2, mbz1, mbz2, -1, dtdx, u1, e33, ff); 
  /*  finite-difference to compute the contribution of e13 to u */
  difference_2z(mbx1, mbx2, mbz1, mbz2, -1, 2.0*dtdz, u1, e13, ll); 

  if (dofct) 
    {
      fct2d1o(u, u1, nx, nz, eta0, eta, deta0dx, deta0dz, 
	      detadx, detadz, dx, dz, mbx1, mbx2, mbz1, mbz2, 
	      f0, f1);
    }

  /*  apply the boundary condition  */
  absorb1(it, nx, nz, dx, dz, dt, u, u1, cc, rho, 
	  isurf, mbx1, mbx2, mbz1, mbz2, 1, 1);
  absorb1(it, nx, nz, dx, dz, dt, u, u1, ll, rho, 
	  isurf, mbx1, mbx2, mbz1, mbz2, 1, 1);

  /*  update the solution  */
  for (ix=mbx1; ix < mbx2; ix++)
    for (iz=mbz1; iz < mbz2; iz++)
      {
	u[ix][iz]=u1[ix][iz];
      }

  /*  solve w at next time step  */
  difference(mbx1-1, mbx2, mbz1, mbz2-1, it, u1, w, 
	     xzsource, fz, dt, rdxdz); 
  /*  finite-difference to compute the contribution of e13 to w */
  difference_2x(mbx1-1, mbx2, mbz1, mbz2-1, 0, 2.0*dtdx, u1, e13, ll); 
  /*  finite-difference to compute the contribution of e11 to w */
  difference_2z(mbx1-1, mbx2, mbz1, mbz2-1, 0, dtdx, u1, e11, ff); 
  /*  finite-difference to compute the contribution of e33 to w */
  difference_2z(mbx1-1, mbx2, mbz1, mbz2-1, 0, dtdz, u1, e33, cc); 

  if (dofct) 
    {
      fct2d1o(w, u1, nx, nz, eta0, eta, deta0dx, deta0dz, 
	      detadx, detadz, dx, dz, mbx1, mbx2-1, mbz1, mbz2-1, 
	      f0, f1);
    }

  /*  apply the boundary condition  */
  absorb1(it, nx, nz, dx, dz, dt, w, u1, cc, rho, 
	  isurf, mbx1, mbx2, mbz1, mbz2-1, 0, 1);

  /*  update the solution  */
  for (ix=mbx1; ix < mbx2-1; ix++)
    for (iz=mbz1; iz < mbz2-1; iz++)
      {
	w[ix][iz]=u1[ix][iz];
      }

  /*  solve v at next time step  */
  difference(mbx1-1, mbx2, mbz1, mbz2-1, it, u1, v, 
	     xzsource, fy, dt, rdxdz); 
  /*  finite-difference to compute the contribution of e12 to v */
  difference_2x(mbx1-1, mbx2, mbz1, mbz2-1, 0, 2.0*dtdx, u1, e12, nn); 
  /*  finite-difference to compute the contribution of e23 to v */
  difference_2z(mbx1-1, mbx2, mbz1, mbz2-1, 0, 2.0*dtdz, u1, e23, ll); 

  if (dofct) 
    {
      fct2d1o(v, u1, nx, nz, eta0, eta, deta0dx, deta0dz, 
	      detadx, detadz, dx, dz, mbx1, mbx2-1, mbz1, mbz2-1, 
	      f0, f1);
    }

  /*  apply the boundary condition  */
  absorb1(it, nx, nz, dx, dz, dt, v, u1, ll, rho, 
	  isurf, mbx1, mbx2-1, mbz1, mbz2-1, 0, 1);

  /*  update the solution  */
  for (ix=mbx1; ix < mbx2-1; ix++)
    for (iz=mbz1; iz < mbz2-1; iz++)
      {
	v[ix][iz]=u1[ix][iz];
      }

  /*  compute strain tensor  */
  /*  solve e11 at next time step  */
  difference(mbx1-1, mbx2, mbz1-1, mbz2+1, it, u1, e11, 
	     xzsource, fz, 0.0, 0.0); 
  difference_2x(mbx1-1, mbx2, mbz1-1, mbz2+1, 0, dtdx, u1, u, rho); 
  /*  no BC applied to e11, size (nx-1)*nz  */
  /*  update the solution e11 */
  for (ix=mbx1; ix < mbx2-1; ix++)
    for (iz=mbz1; iz < mbz2; iz++)
      {
	e11[ix][iz]=u1[ix][iz];
      }

  /*  solve e33 at next time step  */
  difference(mbx1-1, mbx2, mbz1, mbz2, it, u1, e33, 
	     xzsource, fz, 0.0, 0.0); 
  difference_2z(mbx1-1, mbx2, mbz1, mbz2, -1, dtdx, u1, w, rho); 
  /*  apply BC to top and bottom  for e33  */
  /*	absorb1(it, nx, nz, dx, dz, dt, e33, u1, cc, rho, 
	isurf, mbx1, mbx2-1, mbz1, mbz2, 0, 1);
  */	/*  update the solution e33 */
  for (ix=mbx1; ix < mbx2-1; ix++)
    for (iz=mbz1+1; iz < mbz2-1; iz++)
      {
	e33[ix][iz]=u1[ix][iz];
      }

  /*  solve e23 at next time step  */
  difference(mbx1-1, mbx2, mbz1, mbz2, it, u1, e23, 
	     xzsource, fz, 0.0, 0.0); 
  difference_2z(mbx1-1, mbx2, mbz1, mbz2, -1, 0.5*dtdx, u1, v, rho); 
  /*  apply BC to top and bottom for e23  */
  /*	absorb1(it, nx, nz, dx, dz, dt, e23, u1, ll, rho, 
	isurf, mbx1, mbx2-1, mbz1, mbz2, 0, 1);
  */	/*  update the solution e23 */
  for (ix=mbx1; ix < mbx2-1; ix++)
    for (iz=mbz1+1; iz < mbz2-1; iz++)
      {
	e23[ix][iz]=u1[ix][iz];
      }

  /*  solve e12 at next time step  */
  difference(mbx1, mbx2, mbz1-1, mbz2, it, u1, e12, 
	     xzsource, fz, 0.0, 0.0); 
  difference_2x(mbx1, mbx2, mbz1-1, mbz2, -1, 0.5*dtdx, u1, v, rho); 
  /*  apply BC to sides for e12  */
  absorb1(it, nx, nz, dx, dz, dt, e12, u1, ll, rho, 
	  isurf, mbx1, mbx2, mbz1, mbz2-1, 1, 0);
  /*  update the solution e12 */
  for (ix=mbx1; ix < mbx2; ix++)
    for (iz=mbz1; iz < mbz2-1; iz++)
      {
	e12[ix][iz]=u1[ix][iz];
      }

  /*  solve e13 at next time step  */
  difference(mbx1, mbx2, mbz1-1, mbz2, it, u1, e13, 
	     xzsource, fz, 0.0, 0.0); 
  difference_2x(mbx1, mbx2, mbz1-1, mbz2, -1, 0.5*dtdx, u1, w, rho); 
  difference_2z(mbx1, mbx2, mbz1-1, mbz2, 0, 0.5*dtdx, u1, u, rho); 
  /*  apply BC to sides for e13  */
  absorb1(it, nx, nz, dx, dz, dt, e13, u1, cc, rho, 
	  isurf, mbx1, mbx2, mbz1, mbz2-1, 1, 0);
  /*  update the solution e13 */
  for (ix=mbx1; ix < mbx2; ix++)
    for (iz=mbz1; iz < mbz2-1; iz++)
      {
	e13[ix][iz]=u1[ix][iz];
      }

  /*  free the temporery space  */
  free2float(u1);
  free2float(f0);	
  free2float(f1);	
}	


/************************************************************************
* absorb1 -- function to apply boundary condition on one boundary layer,
*	 absorbing boundary condition applied on both sides
*		and bottom.
*************************************************************************
* Notes: 
*		"isurf" can be used to select surface boundary condition	
*
* 		Absorbing boundary conditions on the sides are determined
*	 by the method of Clayton and Engquist, 1980.
*************************************************************************
*************************************************************************
* Input: 
*	int it		time index 
*	int nx		number of grids in x direction 
*	int nz		number of grids in z direction
*	float dx	spatial step in x direction 
*	float dz	spatial step in z direction 
*	float dt	time step	
*	float **u	wavefield at previous time step
*	float **u1	wavefield at current time step
*	float **cc	elastic parameter profile
*	float **rho	density reverse
*	float isurf	1 absorbing BC on surface
*			2 zero gradient BC on surface
*			3 zero value BC on surface
*	int mbx1		left boundary position
*	int mby1		top boundary position
*	int mbx2		right boundary postion 
*	int mby2		bottom boundary position
*	
*************************************************************************
* Output: float **u1	wavefield after applying the boundary condition 
*
*************************************************************************
* Reference: Clayton, R., and Engquist, B., 1980, Absorbing side boundary
*		condition for wave equation migration: Geophysics, vol. 45,
*		p. 895-904.
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	absorb1(int it, int  nx,  int  nz, float dx, 
		float dz, float dt, 
		float **u, float **u1, float **cc, float **rho, 
		int isurf, int mbx1, int mbx2, int mbz1, int mbz2, 
		int side, int tb)
{
  int	i, j;
  float	k, dtdz=dt/dz, dtdx=dt/dx;
  float 	*vell, *velr, *velt, *velb;

  vell = alloc1float(nz);
  velr = alloc1float(nz);
  velt = alloc1float(nx);
  velb = alloc1float(nx);

  /*  compute velocity at the boundary  */
  boundary_vel(cc, rho, vell, velr, velt, velb, 
	       nx, nz, mbx1, mbz1, mbx2, mbz2);

  
  if (tb) {
    for (i=mbx1; i<mbx2; i++)
      {
	k=-dtdz*velt[i];
	if ((k>1)||(k<-1))
	  warn("WARNING:Big k-values.");
	if (isurf == 1)  /* absorbing surface boundary condition */	  {
	    u1[i][mbz1]=u[i][mbz1]-k*(u[i][mbz1+1]-u[i][mbz1]);
	  }
	else if (isurf == 2)  /* zero gradient BC */
	  {
	    u1[i][mbz1]=u1[i][mbz1+1];
	  }	
	else if (isurf == 3)  /*  zero value BC  */
	  {
	    u1[i][mbz1]=0.0;
	  }	

	/* absorbing surface boundary condition */
	k=dtdz*velb[i];
	u1[i][mbz2-1]=u[i][mbz2-1]-k*(u[i][mbz2-1]-u[i][mbz2-2]);
      }
  }
	
  if (side) {
    /* absorbing surface boundary condition */
    for (j=mbz1; j<mbz2; j++)
      {
	k=dtdx*velr[j];
	u1[mbx2-1][j]=u[mbx2-1][j]-k*(u[mbx2-1][j]-u[mbx2-2][j]);	
	k=-dtdx*vell[j];
	u1[mbx1][j]=u[mbx1][j]-k*(u[mbx1+1][j]-u[mbx1][j]);
      }
  }

  free1float(vell);
  free1float(velr);
  free1float(velt);
  free1float(velb);
	
}


/************************************************************************
* absorb2 -- function to apply boundary condition on two boundary layers,
*	 absorbing boundary condition applied on both sides
*		and bottom.
*************************************************************************
* Notes: 
*		"isurf" can be used to select surface boundary condition	
*
* 		Absorbing boundary conditions on the sides are determined
*	 by the method of Clayton and Engquist, 1980.
*************************************************************************
*************************************************************************
* Input: 
*	int it		time index 
*	int nx		number of grids in x direction 
*	int nz		number of grids in z direction
*	float dx	spatial step in x direction 
*	float dz	spatial step in z direction 
*	float dt	time step	
*	float **u	wavefield at previous time step
*	float **u1	wavefield at current time step
*	float **cc	elastic parameter profile
*	float **rho	density reverse
*	float isurf	1 absorbing BC on surface
*			2 zero gradient BC on surface
*			3 zero value BC on surface
*	int mbx1		left boundary position
*	int mby1		top boundary position
*	int mbx2		right boundary postion 
*	int mby2		bottom boundary position
*	
*************************************************************************
* Output: float **u1	wavefield after applying the boundary condition 
*
*************************************************************************
* Reference: Clayton, R., and Engquist, B., 1980, Absorbing side boundary
*		condition for wave equation migration: Geophysics, vol. 45,
*		p. 895-904.
*************************************************************************
* Author: Tong Fei, 1993, Colorado School of Mines.
*************************************************************************/
void	absorb2(int it, int  nx,  int  nz, float dx, 
		float dz, float dt, 
		float **u, float **u1, float **cc, float **rho, 
		int isurf, int mbx1, int mbx2, int mbz1, int mbz2, 
		int side, int tb)
{
  int	i, j;
  float	k, dtdz=dt/dz, dtdx=dt/dx;
  float 	*vell, *velr, *velt, *velb;

  vell = alloc1float(nz);
  velr = alloc1float(nz);
  velt = alloc1float(nx);
  velb = alloc1float(nx);

  /*  compute velocity at the boundary  */
  boundary_vel(cc, rho, vell, velr, velt, velb, 
	       nx, nz, mbx1, mbz1, mbx2, mbz2);

  if (tb) {
    for (i=mbx1+2; i<mbx2-2; i++)
      {
	k=-dtdz*velt[i];
	if (isurf == 1)  /* absorbing surface boundary condition */
	  {
	    u1[i][mbz1]=u[i][mbz1]-k*(u[i][mbz1+1]-u[i][mbz1]);
	    u1[i][mbz1+1]=u[i][mbz1+1]-k*(u[i][mbz1+2]-u[i][mbz1+1]);
	  }
	else if (isurf == 2)  /* zero gradient surface BC */
	  {
	    u1[i][mbz1]=u1[i][mbz1+1];
	    u1[i][mbz1+1]=u[i][mbz1+1]-k*(u[i][mbz1+2]-u[i][mbz1+1]);
	  }	
	else if (isurf == 3)  /* zero value surface BC */
	  {
	    u1[i][mbz1]=0.0;
	    u1[i][mbz1+1]=u[i][mbz1+1]-k*(u[i][mbz1+2]-u[i][mbz1+1]);
	  }	
	/*	absorbing bottom boundary condition	*/	
	k=dtdz*velb[i];
	u1[i][mbz2-1]=u[i][mbz2-1]-k*(u[i][mbz2-1]-u[i][mbz2-2]);
	u1[i][mbz2-2]=u[i][mbz2-2]-k*(u[i][mbz2-2]-u[i][mbz2-3]);
      }
  }
	
  /*	absorbing side boundary condition	*/
  if (side) {	
    for (j=mbz1; j<mbz2; j++)
      {
	k=dtdx*velr[j];
	u1[mbx2-1][j]=u[mbx2-1][j]-k*(u[mbx2-1][j]-u[mbx2-2][j]);	
	u1[mbx2-2][j]=u[mbx2-2][j]-k*(u[mbx2-2][j]-u[mbx2-3][j]);	
	k=-dtdx*vell[j];
	u1[mbx1][j]=u[mbx1][j]-k*(u[mbx1+1][j]-u[mbx1][j]);
	u1[mbx1+1][j]=u[mbx1+1][j]-k*(u[mbx1+2][j]-u[mbx1+1][j]);
      }
  }	

  free1float(vell);
  free1float(velr);
  free1float(velt);
  free1float(velb);
	
}

/************************************************************************
* boundary_vel -- function to compute the velocity at the boundary
*************************************************************************