fdmod.c

地震波正演和显示模块

#include <stdarg.h>
#include <stdio.h>
#include "../include/par.h"
#include "../lib/getpars.c"
#include "../lib/atopkge.c"
#include "../lib/filestat.c"

/*********************** self documentation ***********************/
char *sdoc[] = {
" 									",
" FDMOD - Finite-Difference MODeling (2nd order) for acoustic wave equation",
" 									",
" fdmod <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				",
" dxs=			x coordinates  interval of source				",
" zs=			z coordinates of source				",
" dzs=			z coordinates  interval of source				",
" ns=			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. 								",
" 									",
NULL};

/**************** end self doc ********************************/

#define	ABS0	1
#define	ABS1	1
#define	ABS2	1
#define	ABS3	1
#define O2 0.5000000
#define O6 0.1666667
#define EXIT_SUCCESS 0

/* Prototypes for PML absorbing boundary conditions */
static void pml_init (int nx, int nz, float dx, 
		float dz, float dt, float **dvv);
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);
static void pml_free();

/* PML related global variables */
float pml_max;
int pml_thick;
int pml_thickness;

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);

void xindex (int nx, float ax[], float x, int *index);
void intcub (int ideriv, int nin, float xin[], float ydin[][4],int nout, float xout[], float yout[]);
void warn(char *fmt, ...);
void strchop(char *s, char *t);
void cmonot (int n, float x[], float y[], float yd[][4]);
void requestdoc(int flag);
void pagedoc(void);
void *alloc1 (size_t n1, size_t size);
void **alloc2 (size_t n1, size_t n2, size_t size);
float *alloc1float(size_t n1);
float **alloc2float(size_t n1, size_t n2);
void free1 (void *p);
void free2 (void **p);
void free1float(float *p);
int *alloc1int(size_t n1);
void free1int(int *p);
void free2float(float **p);

int xargc ;char **xargv;

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;		/* minimum density in dfile */
	float dmax;		/* maximum density in dfile */

	float tmax;		/* maximum time to compute */
	float t;		/* time */
	float *xs;		/* array of source x coordinates */
	float dxs;		/* array of source x coordinates interval */
	float *zs;		/* array of source z coordinates */
	float dzs;		/* array of source z coordinates interval */
	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)) warn("must specify nx!");
	if (!getparint("nz",&nz)) warn("must specify nz!");
	if (!getparfloat("tmax",&tmax)) warn("must specify tmax!");

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

	if (ns==0) warn("must specify xs and zs!");
	xs = alloc1float(ns);
	zs = alloc1float(ns);
	ixs = alloc1int(ns);
	izs = alloc1int(ns);
	getparfloat("xs",xs);
	getparfloat("dxs",&dxs);
	getparfloat("zs",zs);
	getparfloat("dzs",&dzs);
        for(is=1;is<ns;is++){
        xs[is]=xs[0]+is*dxs;
        zs[is]=zs[0]+is*dzs;
        }
 
	
	/* 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 )/dz );
	}

	/* 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) {
		if(dvv[ix][iz]>0.000001)	
                        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;
               warn("dt = %g",dt);
		warn("tmax = %g",tmax);
	/* if requested, open file and allocate space for seismograms */
	/* ... horizontal line of seismograms */
	if (*hsfile!='\0') {
		if((hseisfp=fopen(hsfile,"w"))==NULL)
			warn("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)
			warn("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)
			warn("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)
			warn("cannot open dfile=%s",dfile);
		if (fread(od[0],sizeof(float),nx*nz,densityfp)!=nx*nz)
			warn("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);
		if (dmin==dmax) {