fdmod.c

地震波正演和显示模块

}

/* 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

******************************************************************************/
{
	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;
		}
	}
}


/*

   pml_absorb uses the perfectly matched layer absorbing boundary condition.
   The PML formulation is specialized to the acoustic case.



   The size of the new arrays and the location of their (0,0)
   element in the coordinate space of the main p arrays are:


   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        "               "
 
*/

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)
{
        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]);
      }
   }


   /* Calculate w for left edge of bottom pad */

   for (iz=0, ix=nx+pml_thickness-1; iz<pml_thickness+2; ++iz) {
      w_b [ 0][iz] = caz_b [ 0][iz] * w_b [ 0][iz] +
                     cbz_b [ 0][iz] * (ux_b [ix][iz] + uz_b [ix][iz]
                                      -ux_b [ 0][iz] - uz_b [ 0][iz]);
   }


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

   for (ix=1; ix<nx+pml_thickness; ++ix) {
      for (iz=0; iz<pml_thickness+2; ++iz) {
         w_b [ix][iz] = caz_b [ix][iz] * w_b [ix][iz] +
                        cbz_b [ix][iz] * (ux_b [ix-1][iz] + uz_b [ix-1][iz]
                                         -ux_b [ix  ][iz] - uz_b [ix  ][iz]);
      }
   }


   /* Calculate v along top and bottom edge of right pad */

   for (ix=0, jx=nx-1, jz=pml_thickness; ix<pml_thickness+2; ++ix, ++jx) {
      if (jx == nx+pml_thickness) jx = 0;

      v_r [ix][   0] = cax_r [ix][   0] * v_r [ix][   0] +
                       cbx_r [ix][   0] * (ux_r [ix][ 0] + uz_r [ix][ 0]
                                          -ux_b [jx][jz] - uz_b [jx][jz]);

      v_r [ix][nz-1] = cax_r [ix][nz-1] * v_r [ix][nz-1] +
                       cbx_r [ix][nz-1] * (ux_b [jx][   0] + uz_b [jx][   0]
                                          -ux_r [ix][nz-2] - uz_r [ix][nz-2]);
   }


   /* Calculate v in rest of right pad */

   for (ix=0; ix<pml_thickness+2; ++ix) {
      for (iz=1; iz<nz-1; ++iz) {
         v_r [ix][iz] = cax_r [ix][iz] * v_r [ix][iz] +
                        cbx_r [ix][iz] * (ux_r [ix][iz  ] + uz_r [ix][iz  ]
                                         -ux_r [ix][iz-1] - uz_r [ix][iz-1]);
      }
   }


   /* Calculate w along left and right sides of right pad */

   for (iz=0, jx=pml_thickness+2; iz<nz; ++iz) {
      w_r [ 0][iz] = caz_r [ 0][iz] * w_r [ 0][iz] +
                     cbz_r [ 0][iz] * (((abs[3]!=0) ? p [nx-2][iz] : 0.0)
                                      -ux_r [   0][iz] - uz_r [   0][iz]);

      w_r [jx][iz] = caz_r [jx][iz] * w_r [jx][iz] +
                     cbz_r [jx][iz] * (ux_r [jx-1][iz] + uz_r [jx-1][iz]
                                      -((abs[1]!=0) ? p [1][iz] : 0.0));
   }


   /* Calculate w in main part of right pad */

   for (ix=1; ix<pml_thickness+2; ++ix) {