sufdmod2_pml.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* SUFDMOD2_PML: $Revision: 1.7 $ ; $Date: 2006/10/12 18:59:20 $	*/


#include "par.h"
#include "su.h"
#include "segy.h"

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
" SUFDMOD2_PML - Finite-Difference MODeling (2nd order) for acoustic wave",
"    equation with PML absorbing boundary conditions.			",
" Caveat: experimental PML absorbing boundary condition version,	",
"may be buggy!								",
" 									",
" sufdmod2_pml <vfile >wfile nx= nz= tmax= xs= zs= [optional parameters]",
" 									",
" Required Parameters:							",
" <vfile		file containing velocity[nx][nz]		",
" >wfile		file containing waves[nx][nz] for time steps	",
" nx=			number of x samples (2nd dimension)		",
" nz=			number of z samples (1st dimension)		",
" xs=			x coordinates of source				",
" zs=			z coordinates of source				",
" tmax=			maximum time					",
" 									",
" Optional Parameters:							",
" nt=1+tmax/dt		number of time samples (dt determined for stability)",
" mt=1			number of time steps (dt) per output time step	",
" 									",
" dx=1.0		x sampling interval				",
" fx=0.0		first x sample					",
" dz=1.0		z sampling interval				",
" fz=0.0		first z sample					",
" 									",
" fmax = vmin/(10.0*h)	maximum frequency in source wavelet		",
" fpeak=0.5*fmax	peak frequency in ricker wavelet		",
" 									",
" dfile=		input file containing density[nx][nz]		",
" vsx=			x coordinate of vertical line of seismograms	",
" hsz=			z coordinate of horizontal line of seismograms	",
" vsfile=		output file for vertical line of seismograms[nz][nt]",
" hsfile=		output file for horizontal line of seismograms[nx][nt]",
" ssfile=		output file for source point seismograms[nt]	",
" verbose=0		=1 for diagnostic messages, =2 for more		",
" 									",
" abs=1,1,1,1		Absorbing boundary conditions on top,left,bottom,right",
" 			sides of the model. 				",
" 		=0,1,1,1 for free surface condition on the top		",
"                                                                       ",
" ...PML parameters....                                                 ",
" pml_max=1000.0        PML absorption parameter                        ",
" pml_thick=0           half-thickness of pml layer (0 = do not use PML)",
" 									",
" Notes:								",
" This program uses the traditional explicit second order differencing	",
" method. 								",
" 									",
" Two different absorbing boundary condition schemes are available. The ",
" first is a traditional absorbing boundary condition scheme created by ",
" Hale, 1990. The second is based on the perfectly matched layer (PML)	",
" method of Berenger, 1995.						",
" 									",
NULL};

/*
 * Authors:  CWP:Dave Hale
 *           CWP:modified for SU by John Stockwell, 1993.
 *           CWP:added frequency specification of wavelet: Craig Artley, 1993
 *           TAMU:added PML absorbing boundary condition: 
 *               Michael Holzrichter, 1998
 *           CWP/WesternGeco:corrected PML code to handle density variations:
 *               Greg Wimpey, 2006
 *
 * References: (Hale's absobing boundary conditions)
 * 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.
 *
 * References: (PML boundary conditions)
 * 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.
 *
 *
 * Trace header fields set: ns, delrt, tracl, tracr, offset, d1, d2,
 *                          sdepth, trid
 */
/**************** end self doc ********************************/

#define	ABS0	1
#define	ABS1	1
#define	ABS2	1
#define	ABS3	1

/* Prototypes for PML absorbing boundary conditions */
static void pml_init (int nx, int nz, float dx, 
		float dz, float dt, float **dvv, float **od, int verbose);
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 related global variables */
float pml_max=0;
int pml_thick=0;
int pml_thickness=0;

float **cax_b, **cax_r;
float **cbx_b, **cbx_r;
float **caz_b, **caz_r;
float **cbz_b, **cbz_r;
float **dax_b, **dax_r;
float **dbx_b, **dbx_r;
float **daz_b, **daz_r;
float **dbz_b, **dbz_r;

float **ux_b,  **ux_r;
float **uz_b,  **uz_r;
float **v_b,   **v_r;
float **w_b,   **w_r;

float dvv_0, dvv_1, dvv_2, dvv_3;
float sigma, sigma_ex, sigma_ez, sigma_mx, sigma_mz;


/* Prototypes for finite differencing */
void ptsrc (float xs, float zs,
	int nx, float dx, float fx,
	int nz, float dz, float fz,
	float dt, float t, float fmax, float fpeak, float tdelay, float **s);
void exsrc (int ns, float *xs, float *zs,
	int nx, float dx, float fx,
	int nz, float dz, float fz,
	float dt, float t, float fmax, float **s);
void tstep2 (int nx, float dx, int nz, float dz, float dt,
	float **dvv, float **od, float **s,
	float **pm, float **p, float **pp, int *abs);

/* Globals for trace manipulation */
segy cubetr; 	/* data cube traces */
segy srctr;	/* source seismogram traces */
segy horiztr;	/* horizontal line seismogram traces */
segy verttr;	/* vertical line seismogram traces */

int
main(int argc, char **argv)
{
	int ix,iz,it,is;	/* counters */
	int nx,nz,nt,mt;	/* x,z,t,tsizes */

	int verbose;		/* is verbose? */
	int nxs;		/* number of source x coordinates */
	int nzs;		/* number of source y coordinates */
	int ns;			/* total number of sources ns=nxs=nxz */

	int vs2;		/* depth in samples of horiz rec line */
	int hs1;		/* horiz sample of vert rec line */

	float fx;		/* first x value */
	float dx;		/* x sample interval */

	float fz;		/* first z value */
	float dz;		/* z sample interval */
	float h;		/* minumum spatial sample interval */

	float hsz;		/* z position of horiz receiver line */
	float vsx;		/* x position of vertical receiver line */

	float dt;		/* time sample interval */
	float fmax;		/* maximum temporal frequency allowable */
	float fpeak;		/* peak frequency of ricker wavelet */
	float tdelay=0.;	/* time delay of source beginning */

	float vmin;		/* minimum wavespeed in vfile */
	float vmax;		/* maximum wavespeed in vfile */

	float dmin=0.0;		/* minimum density in dfile */
	float dmax=0.0;		/* maximum density in dfile */

	float tmax;		/* maximum time to compute */
	float t;		/* time */
	float *xs;		/* array of source x coordinates */
	float *zs;		/* array of source z coordinates */
	int *ixs;		/* array of source x sample locations */
	int *izs;		/* array of source z sample locations */
	float **s;		/* array of source pressure values */
	float **dvv;		/* array of velocity values from vfile */
	float **od;		/* array of density values from dfile */

	/* pressure field arrays */
	float **pm;		/* pressure field at t-dt */
	float **p;		/* pressure field at t */
	float **pp;		/* pressure field at t+dt */
	float **ptemp;		/* temp pressure array */

	/* output data arrays */
	float **ss;		/* source point seismogram array */
	float **hs;		/* seismograms from horiz receiver line */
	float **vs;		/* seismograms from vert receiver line */

	/* file names */
	char *dfile="";		/* density file name */
	char *vsfile="";	/* vert receiver seismogram line file  name */
	char *hsfile="";	/* horiz receiver seismogram line file name */
	char *ssfile="";	/* source point seismogram file name */

	/* input file pointers */
	FILE *velocityfp=stdin; /* pointer to input velocity data */
	FILE *densityfp;	/* pointer to input density data file */

	/* output file pointers */
	FILE *hseisfp=NULL;	/* pointer to output horiz rec line file  */
	FILE *vseisfp=NULL;	/* pointer to output vert rec line file  */	
	FILE *sseisfp=NULL;	/* pointer to source point seis output file */

	/* SEGY fields */
	long tracl=0;		/* trace number within a line */
	long tracr=0;		/* trace number within a reel */

	/* Absorbing boundary conditions related stuff*/
	int abs[4];		/* absorbing boundary cond. flags */
	int nabs;		/* number of values given */

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);
	
	/* get required parameters */
	/* get dimensions of model, maximum duration */
	if (!getparint("nx",&nx)) err("must specify nx!");
	if (!getparint("nz",&nz)) err("must specify nz!");
	if (!getparfloat("tmax",&tmax)) err("must specify tmax!");

	/* get source information, coordinates */
	nxs = countparval("xs");
	nzs = countparval("zs");
	if (nxs!=nzs)
		err("number of xs = %d must equal number of zs = %d",
			nxs,nzs);
	ns = nxs;

	if (ns==0) err("must specify xs and zs!");
	xs = alloc1float(ns);
	zs = alloc1float(ns);
	ixs = alloc1int(ns);
	izs = alloc1int(ns);
	getparfloat("xs",xs);
	getparfloat("zs",zs);
	
	/* Get absorbing boundary information */
	nabs = countparval("abs");
	if (nabs==4) {
		getparint("abs", abs);	
	} else {
		abs[0] = ABS0;
		abs[1] = ABS1;
		abs[2] = ABS2;
		abs[3] = ABS3;

		if (!((nabs==4) || (nabs==0)) ) 
			warn("Number of abs %d, using abs=1,1,1,1",nabs);
	}
	
	/* get optional parameters */
	if (!getparint("nt",&nt)) nt = 0;
	if (!getparint("mt",&mt)) mt = 1;
	if (!getparfloat("dx",&dx)) dx = 1.0;
	if (!getparfloat("fx",&fx)) fx = 0.0;
	if (!getparfloat("dz",&dz)) dz = 1.0;
	if (!getparfloat("fz",&fz)) fz = 0.0;

        if (!getparfloat("pml_max",&pml_max)) pml_max = 1000.0;
        if (!getparint("pml_thick",&pml_thick)) pml_thick = 0;
        pml_thickness = 2 * pml_thick;

	/* source coordinates in samples */
	for (is=0 ; is < ns ; ++is) {
		ixs[is] = NINT( ( xs[is] - fx )/dx );
		izs[is] = NINT( ( zs[is] - fz )/dx );
	}

	/* z-coorinate of horizontal line of detectors */
	if (!getparfloat("hsz",&hsz)) hsz = 0.0;
	hs1 = NINT( (hsz - fz)/dz );

	/* x-coordinate of vertical line of detectors */
	if (!getparfloat("vsx",&vsx)) vsx = 0.0;
	vs2 = NINT((vsx - fx)/dx );
	
	if (!getparint("verbose",&verbose)) verbose = 0;

	/* Input and output file information */
	getparstring("dfile",&dfile);
	getparstring("hsfile",&hsfile);
	getparstring("vsfile",&vsfile);
	getparstring("ssfile",&ssfile);
	
	/* allocate space */
	s = alloc2float(nz,nx);
	dvv = alloc2float(nz,nx);
	od = alloc2float(nz,nx);
	pm = alloc2float(nz,nx);
	p = alloc2float(nz,nx);
	pp = alloc2float(nz,nx);
	
	/* read velocities */
	fread(dvv[0],sizeof(float),nx*nz,velocityfp);
	
	/* determine minimum and maximum velocities */
	vmin = vmax = dvv[0][0];
	for (ix=0; ix<nx; ++ix) {
		for (iz=0; iz<nz; ++iz) {
			vmin = MIN(vmin,dvv[ix][iz]);
			vmax = MAX(vmax,dvv[ix][iz]);
		}
	}
	
	/* determine mininum spatial sampling interval */
	h = MIN(ABS(dx),ABS(dz));
	
	/* determine time sampling interval to ensure stability */
	dt = h/(2.0*vmax);
	
	/* determine maximum temporal frequency to avoid dispersion */
	if (!getparfloat("fmax", &fmax))	fmax = vmin/(10.0*h);

	/* compute or set peak frequency for ricker wavelet */
	if (!getparfloat("fpeak", &fpeak))	fpeak = 0.5*fmax;

	/* determine number of time steps required to reach maximum time */
	if (nt==0) nt = 1+tmax/dt;

	/* if requested, open file and allocate space for seismograms */
	/* ... horizontal line of seismograms */
	if (*hsfile!='\0') {
		if((hseisfp=fopen(hsfile,"w"))==NULL)
			err("cannot open hsfile=%s",hsfile);
		hs = alloc2float(nt,nx);
	} else {
		hs = NULL;
	}

	/* ... vertical line of seismograms */
	if (*vsfile!='\0') {
		if((vseisfp=fopen(vsfile,"w"))==NULL)
			err("cannot open vsfile=%s",vsfile);
		vs = alloc2float(nt,nz);
	} else {
		vs = NULL;
	}

	/* ...  seismograms at the source point */
	if (*ssfile!='\0') {
		if((sseisfp=fopen(ssfile,"w"))==NULL)
			err("cannot open ssfile=%s",ssfile);
		ss = alloc2float(nt,ns);
	} else {
		ss = NULL;
	}
	
	/* if specified, read densities */
	if (*dfile!='\0') {
		if((densityfp=fopen(dfile,"r"))==NULL)
			err("cannot open dfile=%s",dfile);
		if (fread(od[0],sizeof(float),nx*nz,densityfp)!=nx*nz)
			err("error reading dfile=%s",dfile);
		fclose(densityfp);
		dmin = dmax = od[0][0];
		for (ix=0; ix<nx; ++ix) {
			for (iz=0; iz<nz; ++iz) {
				dmin = MIN(dmin,od[ix][iz]);
				dmax = MAX(dmax,od[ix][iz]);
			}
		}
	}
	
	/* if densities not specified or constant, make densities = 1 */
	if (*dfile=='\0' || dmin==dmax ) {
		for (ix=0; ix<nx; ++ix)
			for (iz=0; iz<nz; ++iz)
				od[ix][iz] = 1.0;
		dmin = dmax = 1.0;
	}
	
	/* compute density*velocity^2 and 1/density and zero time slices */	
	for (ix=0; ix<nx; ++ix) {
		for (iz=0; iz<nz; ++iz) {
			dvv[ix][iz] = od[ix][iz]*dvv[ix][iz]*dvv[ix][iz];
			od[ix][iz] = 1.0/od[ix][iz];
			pp[ix][iz] = p[ix][iz] = pm[ix][iz] = 0.0;
		}
	}
	
	/* if densities constant, free space and set NULL pointer */
	if (dmin==dmax) {
		free2float(od);
		od = NULL;
	}
	
	/* if verbose, print parameters */
	if (verbose) {
		warn("nx = %d",nx);
		warn("dx = %g",dx);
		warn("nz = %d",nz);
		warn("dz = %g",dz);
		warn("nt = %d",nt);
		warn("dt = %g",dt);
		warn("tmax = %g",tmax);
		warn("fmax = %g",fmax);
		warn("fpeak = %g",fpeak);
		warn("vmin = %g",vmin);
		warn("vmax = %g",vmax);
		warn("mt = %d",mt);
                warn("pml_max = %g",pml_max);
                warn("pml_half = %d",pml_thick);
                warn("pml_thickness = %d",pml_thickness);
		if (dmin==dmax) {
			warn("constant density");
		} else {
			warn("dfile=%s",dfile);
			warn("dmin = %g",dmin);
			warn("dmax = %g",dmax);