sufdmod2_pml.c

su 的源代码库

	int ix,iz;
	float xscale1,zscale1,xscale2,zscale2;
		
	/* determine constants */
	xscale1 = (dt*dt)/(dx*dx);
	zscale1 = (dt*dt)/(dz*dz);
	xscale2 = 0.25*xscale1;
	zscale2 = 0.25*zscale1;
	
	/* do the finite-difference star */
	for (ix=1; ix<nx-1; ++ix) {
		for (iz=1; iz<nz-1; ++iz) {
			pp[ix][iz] = 2.0*p[ix][iz]-pm[ix][iz] +
				dvv[ix][iz]*(
					od[ix][iz]*(
						xscale1*(
							p[ix+1][iz]+
							p[ix-1][iz]-
							2.0*p[ix][iz]
						) +
						zscale1*(
							p[ix][iz+1]+
							p[ix][iz-1]-
							2.0*p[ix][iz]
						)
					) +
					(
						xscale2*(
							(od[ix+1][iz]-
							od[ix-1][iz]) *
							(p[ix+1][iz]-
							p[ix-1][iz])
						) +
						zscale2*(
							(od[ix][iz+1]-
							od[ix][iz-1])*
							(p[ix][iz+1]-
							p[ix][iz-1])
						)
					)
				) +
				s[ix][iz];
		}
	}
}

/* convolve with finite-difference star for variable density and dx==dz */
static void star2 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp)
{
	int ix,iz;
	float scale1,scale2;
	
	if ( dx != dz ) 
		warn("ASSERT FAILED: dx != dz in star2");
	/* determine constants */
	scale1 = (dt*dt)/(dx*dx);
	scale2 = 0.25*scale1;
	
	/* do the finite-difference star */
	for (ix=1; ix<nx-1; ++ix) {
		for (iz=1; iz<nz-1; ++iz) {
			pp[ix][iz] = 2.0*p[ix][iz]-pm[ix][iz] +
				dvv[ix][iz]*(
					od[ix][iz]*(
						scale1*(
							p[ix+1][iz]+
							p[ix-1][iz]+
							p[ix][iz+1]+
							p[ix][iz-1]-
							4.0*p[ix][iz]
						)
					) +
					(
						scale2*(
							(od[ix+1][iz]-
							od[ix-1][iz]) *
							(p[ix+1][iz]-
							p[ix-1][iz]) +
							(od[ix][iz+1]-
							od[ix][iz-1]) *
							(p[ix][iz+1]-
							p[ix][iz-1])
						)
					)
				) +
				s[ix][iz];
		}
	}
}

/* convolve with finite-difference star for density==1.0 and dx!=dz */
static void star3 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp)
{
	int ix,iz;
	float xscale,zscale;
		
	if ( od != ((float **) NULL) ) 
		warn("ASSERT FAILED: od !=  NULL in star3");
	/* determine constants */
	xscale = (dt*dt)/(dx*dx);
	zscale = (dt*dt)/(dz*dz);
	
	/* do the finite-difference star */
	for (ix=1; ix<nx-1; ++ix) {
		for (iz=1; iz<nz-1; ++iz) {
			pp[ix][iz] = 2.0*p[ix][iz]-pm[ix][iz] +
				dvv[ix][iz]*(
					xscale*(
						p[ix+1][iz]+
						p[ix-1][iz]-
						2.0*p[ix][iz]
					) +
					zscale*(
						p[ix][iz+1]+
						p[ix][iz-1]-
						2.0*p[ix][iz]
					)
				) +
				s[ix][iz];
		}
	}
}

/* convolve with finite-difference star for density==1.0 and dx==dz */
static void star4 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp)
{
	int ix,iz;
	float scale;
	
	/* determine constants */
	if ( od != ((float **) NULL) ) 
		warn("ASSERT FAILED: od !=  NULL in star4");
	if ( dz != dx ) 
		warn("ASSERT FAILED: dz !=  dx in star4");
	scale = (dt*dt)/(dx*dz);
	
	/* do the finite-difference star */
	for (ix=1; ix<nx-1; ++ix) {
		for (iz=1; iz<nz-1; ++iz) {
			pp[ix][iz] = 2.0*p[ix][iz]-pm[ix][iz] +
				scale*dvv[ix][iz]*(
					p[ix+1][iz]+
					p[ix-1][iz]+
					p[ix][iz+1]+
					p[ix][iz-1]-
					4.0*p[ix][iz]
				) +
				s[ix][iz];
		}
	}
}

static void absorb (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **pm, float **p, float **pp,
	int *abs)
/*****************************************************************************
absorb - absorbing boundary conditions 
*****************************************************************************
Input:
nx	number of samples in x direction
dx	spatial sampling interval in x direction
nz	number of samples in z direction
dz	spatial sampling interval in z direction
dt	time sampling interval
dvv	array of velocity values from vfile
od	array of density values from dfile
pm	pressure field at time t-1
p	pressure field at time t
pp	pressure field at t+dt
abs	flag indicating to absorb or not to absorb

*****************************************************************************
Notes:
This method is an improvement on the method of Clayton and Engquist, 1977
and 1980. The method is described in Hale, 1990.

*****************************************************************************
References:
Clayton, R. W., and Engquist, B., 1977, Absorbing boundary conditions
for acoustic and elastic wave equations, Bull. Seism. Soc. Am., 6, 1529-1540. 

Clayton, R. W., and Engquist, B., 1980, Absorbing boundary conditions
for wave equation migration, Geophysics, 45, 895-904.

Hale, D.,  1990, Adaptive absorbing boundaries for finite-difference
modeling of the wave equation migration, unpublished report from the
Center for Wave Phenomena, Colorado School of Mines.

Richtmyer, R. D., and Morton, K. W., 1967, Difference methods for
initial-value problems, John Wiley & Sons, Inc, New York.

Thomee, V., 1962, A stable difference scheme for the mixed boundary problem
for a hyperbolic, first-order system in two dimensions, J. Soc. Indust.
Appl. Math., 10, 229-245.

Toldi, J. L., and Hale, D., 1982, Data-dependent absorbing side boundaries,
Stanford Exploration Project Report SEP-30, 111-121.

*****************************************************************************
Author: CWP: Dave Hale 1990 
******************************************************************************/
{
	int ix,iz;
	float ov,ovs,cosa,beta,gamma,dpdx,dpdz,dpdt,dpdxs,dpdzs,dpdts;

	/* solve for upper boundary */
	iz = 1;
	for (ix=0; ix<nx; ++ix) {

		if (abs[0]!=0) {

			if (od!=NULL)
				ovs = 1.0/(od[ix][iz]*dvv[ix][iz]);
			else
				ovs = 1.0/dvv[ix][iz];
			ov = sqrt(ovs);
			if (ix==0)
				dpdx = (p[1][iz]-p[0][iz])/dx;
			else if (ix==nx-1)
				dpdx = (p[nx-1][iz]-p[nx-2][iz])/dx;
			else
				dpdx = (p[ix+1][iz]-p[ix-1][iz])/(2.0*dx);
			dpdt = (pp[ix][iz]-pm[ix][iz])/(2.0*dt);
			dpdxs = dpdx*dpdx;
			dpdts = dpdt*dpdt;
			if (ovs*dpdts>dpdxs)
				cosa = sqrt(1.0-dpdxs/(ovs*dpdts));
			else 
				cosa = 0.0;

			beta = ov*dz/dt*cosa;
			gamma = (1.0-beta)/(1.0+beta);

			pp[ix][iz-1] = gamma*(pp[ix][iz]-p[ix][iz-1])+p[ix][iz];
		} else {
			pp[ix][iz-1] = 0.0;
		}
	}


	/* extrapolate along left boundary */
	ix = 1;
	for (iz=0; iz<nz; ++iz) {
		if (abs[1]!=0) {
			if (od!=NULL)
				ovs = 1.0/(od[ix][iz]*dvv[ix][iz]);
			else
				ovs = 1.0/dvv[ix][iz];
			ov = sqrt(ovs);
			if (iz==0)
				dpdz = (p[ix][1]-p[ix][0])/dz;
			else if (iz==nz-1)
				dpdz = (p[ix][nz-1]-p[ix][nz-2])/dz;
			else
				dpdz = (p[ix][iz+1]-p[ix][iz-1])/(2.0*dz);
			dpdt = (pp[ix][iz]-pm[ix][iz])/(2.0*dt);
			dpdzs = dpdz*dpdz;
			dpdts = dpdt*dpdt;
			if (ovs*dpdts>dpdzs)
				cosa = sqrt(1.0-dpdzs/(ovs*dpdts));
			else
				cosa = 0.0;

			beta = ov*dx/dt*cosa;
			gamma = (1.0-beta)/(1.0+beta);
			pp[ix-1][iz] = gamma*(pp[ix][iz]-p[ix-1][iz])+p[ix][iz];
		} else {
			pp[ix-1][iz] = 0.0;
		}
	}

	/* extrapolate along lower boundary */
	iz = nz-2;
	for (ix=0; ix<nx; ++ix) {
		if (abs[2]!=0) {
			if (od!=NULL)
				ovs = 1.0/(od[ix][iz]*dvv[ix][iz]);
			else
				ovs = 1.0/dvv[ix][iz];
			ov = sqrt(ovs);
			if (ix==0)
				dpdx = (p[1][iz]-p[0][iz])/dx;
			else if (ix==nx-1)
				dpdx = (p[nx-1][iz]-p[nx-2][iz])/dx;
			else
				dpdx = (p[ix+1][iz]-p[ix-1][iz])/(2.0*dx);
			dpdt = (pp[ix][iz]-pm[ix][iz])/(2.0*dt);
			dpdxs = dpdx*dpdx;
			dpdts = dpdt*dpdt;
			if (ovs*dpdts>dpdxs)
				cosa = sqrt(1.0-dpdxs/(ovs*dpdts));
			else 
				cosa = 0.0;

			beta = ov*dz/dt*cosa;
			gamma = (1.0-beta)/(1.0+beta);

			pp[ix][iz+1] = gamma*(pp[ix][iz]-p[ix][iz+1])+p[ix][iz];
		} else {
			pp[ix][iz+1] = 0.0;
		}
	}

	/* extrapolate along right boundary */
	ix = nx-2;
	for (iz=0; iz<nz; ++iz) {
		if (abs[3]!=0) {
			if (od!=NULL)
				ovs = 1.0/(od[ix][iz]*dvv[ix][iz]);
			else
				ovs = 1.0/dvv[ix][iz];
			ov = sqrt(ovs);
			if (iz==0)
				dpdz = (p[ix][1]-p[ix][0])/dz;
			else if (iz==nz-1)
				dpdz = (p[ix][nz-1]-p[ix][nz-2])/dz;
			else
				dpdz = (p[ix][iz+1]-p[ix][iz-1])/(2.0*dz);
			dpdt = (pp[ix][iz]-pm[ix][iz])/(2.0*dt);
			dpdzs = dpdz*dpdz;
			dpdts = dpdt*dpdt;
			if (ovs*dpdts>dpdzs)
				cosa = sqrt(1.0-dpdzs/(ovs*dpdts));
			else
				cosa = 0.0;

			beta = ov*dx/dt*cosa;
			gamma = (1.0-beta)/(1.0+beta);
			pp[ix+1][iz] =gamma*(pp[ix][iz]-p[ix+1][iz])+p[ix][iz];
		} else {
			pp[ix+1][iz] = 0.0;
		}
	}
}


static void pml_absorb (int nx, float dx, int nz, float dz, float dt,
        float **dvv, float **od, float **pm, float **p, float **pp,
        int *abs)
/**************************************************************************
  pml_absorb - uses the perfectly matched layer absorbing boundary condition.
***************************************************************************
Notes:
   The PML formulation is specialized to the acoustic case.


   Array        Size        Location of [0][0] with respect to p [0][0]
   -----   ---------------  -------------------------------------------
   ux_b    (nx+pml, pml+2)  (0, nz-1)
   uz_b         "               "
   dax_b        "               "
   dbx_b        "               "
   daz_b        "               "
   dbz_b        "               "
 
   ux_r    (pml+2, nz)      (nx-1, 0)
   uz_r         "               "
   dax_r        "               "
   dbx_r        "               "
   daz_r        "               "
   dbz_r        "               "

   v_b     (nx+pml, pml+3)  (0, nz-1.5)
   cax_b        "               "
   cbx_b        "               "

   w_b     (nx+pml, pml+2)  (-0.5, nz-1)
   caz_b        "               "
   cbz_b        "               "

   v_r     (pml+2, nz)      (nx-1, -0.5)
   cax_r        "               "
   cbx_r        "               "

   w_r     (pml+3, nz)      (nx-1.5, 0)
   caz_r        "               "
   cbz_r        "               "
***************************************************************************
References:
   Jean-Pierre Berenger, ``A Perfectly Matched Layer for the Absorption of
   Electromagnetic Waves,''  Journal of Computational Physics, vol. 114,
   pp. 185-200.

   Hastings, Schneider, and Broschat, ``Application of the perfectly
   matched layer (PML) absorbing boundary condition to elastic wave
   propogation,''  Journal of the Accoustical Society of America,
   November, 1996.

   Allen Taflove, ``Electromagnetic Modeling:  Finite Difference Time
   Domain Methods'', Baltimore, Maryland: Johns Hopkins University Press,
   1995, chap. 7, pp. 181-195.

   The PML ABC is implemented by extending the modeled region on
   the bottom and right sides and treating the modeled region as
   periodic.

   In the extended region, the differential equations of PML are
   modeled.  The extension is accomplished by using additional
   arrays which record the state in the extended regions.  The
   result is a nasty patchwork of arrays.  (It is possible to use
   the PML differential equations to model the absorbing and
   non-absorbing regions.  This greatly simplifies things at the
   expense of memory.)

   The size of the new arrays and the location of their (0,0)
   element in the coordinate space of the main p arrays are:
*************************************************************************
   Author:  TAMU: Michael Holzrichter, 1998
*************************************************************************/
{
        int ix, iz, jx, jz;

   /* Calculate v for bottom pad above and below main domain */

   for (ix=0, jz=pml_thickness+2; ix<nx; ++ix) {
      v_b[ix][ 0] = cax_b [ix][ 0] * v_b [ix][ 0] +
                      cbx_b [ix][ 0] * (ux_b [ix][   0] + uz_b [ix][   0]
                                       -((abs[2]!=0) ? p [ix][nz-2] : 0.0));

      v_b[ix][jz] = cax_b [ix][jz] * v_b [ix][jz] +
                      cbx_b [ix][jz] * (((abs[0]!=0) ? p [ix][   1] : 0.0)
                                       -ux_b [ix][jz-1] - uz_b [ix][jz-1]);
   }	

   /* Calculate v for bottom pad above and below right pad */

   for (ix=nx, jx=1, jz=pml_thickness+2; ix<nx+pml_thickness; ++ix, ++jx) {
      v_b  [ix][ 0] = cax_b [ix][ 0] * v_b [ix][ 0] +
                      cbx_b [ix][ 0] * (ux_b [ix][   0] + uz_b [ix][   0]
                                       -ux_r [jx][nz-2] - uz_r [jx][nz-2]);

      v_b  [ix][jz] = cax_b [ix][jz] * v_b [ix][jz] +
                      cbx_b [ix][jz] * (ux_r [jx][   1] + uz_r [jx][   1]
                                       -ux_b [ix][jz-1] - uz_b [ix][jz-1]);
   }


   /* Calculate v for main part of bottom pad */

   for (ix=0; ix<nx+pml_thickness; ++ix) {
      for (iz=1; iz<pml_thickness+2; ++iz) {
         v_b [ix][iz] = cax_b [ix][iz] * v_b [ix][iz] +
                        cbx_b [ix][iz] * (ux_b [ix][iz  ] + uz_b [ix][iz  ]
                                         -ux_b [ix][iz-1] - uz_b [ix][iz-1]);
      }