sufctanismod.c

su 的源代码库

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

/* SUFCTANISMOD: $Revision: 1.3 $ ; $Date: 2006/10/31 22:18:47 $	*/

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

/************************** self documentation *************************/
char *sdoc[] = {
"									",
" SUFCTANISMOD - Flux-Corrected Transport correction applied to the 2D",
"	  elastic wave equation for finite difference modeling in 	",
"	  anisotropic media						",
"									",
" sufctanismod > outfile [optional parameters]			",
"		outfile is the wavefield snapshot x-component		",
"		x-component of wavefield snapshot is in snapshotx.data	",
"		y-component of wavefield snapshot is in snapshoty.data	",
"		z-component of wavefield snapshot is in snapshotz.data	",
" 									",
" Optional Output Files:						",
" reflxfile=	reflection seismogram file name for x-component		",
"		no output produced if no name specified	 		",
" reflyfile=	reflection seismogram file name for y-component		",
"		no output produced if no name specified	 		",
" reflzfile=	reflection seismogram file name for z-component		",
"		no output produced if no name specified	 		",
" vspxfile=	VSP seismogram file name for x-component		",
"		no output produced if no name specified	 		",
" vspyfile=	VSP seismogram file name for y-component		",
"		no output produced if no name specified	 		",
" vspzfile=	VSP seismogram file name for z-component		",
"		no output produced if no name specified	 		",
"									",
" suhead=1      To get SU-header output seismograms (else suhead=0)	",
"									",
" New parameter:							",
" receiverdepth=0  depth of horizontal receivers (in gridpoints)      ",
"									",
"	               							",     
" Optional Parameters:						",
" dofct=1 	1 do the FCT correction					",
"		0 do not do the FCT correction 				",
" FCT Related parameters:						",
" eta0=0.03	diffusion coefficient					",
"		typical values ranging from 0.008 to 0.06		",
"		about 0.03 for the second-order method 			",
"		about 0.012 for the fourth-order method 		",
" eta=0.04	anti-diffusion coefficient 				",
"		typical values ranging from 0.008 to 0.06		",
"		about 0.04 for the second-order method  		",
"		about 0.015 for the fourth-order method 		",
" fctxbeg=0 	x coordinate to begin applying the FCT correction	",
" fctzbeg=0 	z coordinate to begin applying the FCT correction	",
" fctxend=nx 	x coordinate to stop applying the FCT correction	",
" fctzend=nz 	z coordinate to stop applying the FCT correction	",
" 									",
" deta0dx=0.0	gradient of eta0 in x-direction  d(eta0)/dx		",
" deta0dz=0.0	gradient of eta0 in z-direction  d(eta0)/dz		",
" detadx=0.0	gradient of eta in x-direction 	 d(eta)/dx		",
" detadz=0.0	gradient of eta in z-direction 	 d(eta)/dz		",
"									",
" General Parameters:							",
" order=2	2 second-order finite-difference 			",
"		4 fourth-order finite-difference 			",
"									",
" nt=200        number of time steps 			 		",
" dt=0.004	time step  						",
"									",
" nx=100 	number of grid points in x-direction 			",
" nz=100 	number of grid points in z-direction 			",
"									",
" dx=0.02	spatial step in x-direction 				",
" dz=0.02	spatial step in z-direction 				",
"									",
" sx=nx/2	source x-coordinate (in gridpoints)			",
" sy=nz/2	source z-coordinate (in gridpoints)			",
"									",
" fpeak=20	peak frequency of the wavelet 				",
"									",
" wavelet=1	1 AKB wavelet						",
" 		2 Ricker wavelet 					",
"		3 impulse 						",
"		4 unity 						",
"									",
" isurf=2	1 absorbing surface condition 				",
"		2 free surface condition 				",
"		3 zero surface condition 				",
"									",
" source=1	1 point source 						",
" 		2 sources are located on a given refelector 	        ", 
"			(two horizontal and one dipping reflectors) 	",
" 		3 sources are located on a given dipping refelector     ", 
"									",
" sfile= 	the name of input source file, if no name specified then",
"		use default source location. (source=1 or 2) 		",
"									",
" Density and Elastic Parameters:					",
" dfile= 	the name of input density file,                         ",
"               if no name specified then                             ",
"		assume a linear density profile with ...		",
" rho00=2.0	density at (0, 0) 					",
" drhodx=0.0	density gradient in x-direction  d(rho)/dx		",
" drhodz=0.0	density gradient in z-direction  d(rho)/dz		",
"									",
" afile= 	name of input elastic param.  (c11) aa file, if no name ",
"		specified then, assume a linear profile with ...	",
" aa00=2.0	elastic parameter at (0, 0) 				",
" daadx=0.0	parameter gradient in x-direction  d(aa)/dx		",
" daadz=0.0	parameter gradient in z-direction  d(aa)/dz		",
"									",
" cfile= 	name of input elastic param. (c33)  cc file, if no name ",
"		specified then, assume a linear profile with ...	",
" cc00=2.0	elastic parameter at (0, 0) 				",
" dccdx=0.0	parameter gradient in x-direction  d(cc)/dx		",
" dccdz=0.0	parameter gradient in z-direction  d(cc)/dz		",
"									",
" ffile= 	name of input elastic param.  (c13) ff file, if no name ",
"		specified then, assume a linear profile with ...	",
" ff00=2.0	elastic parameter at (0, 0) 				",
" dffdx=0.0	parameter gradient in x-direction  d(ff)/dx		",
" dffdz=0.0	parameter gradient in z-direction  d(ff)/dz		",
"									",
" lfile= 	name of input elastic param.  (c44) ll file, if no name ",
"		specified then, assume a linear profile with ...	",
" ll00=2.0	elastic parameter at (0, 0) 				",
" dlldx=0.0	parameter gradient in x-direction  d(ll)/dx		",
" dlldz=0.0	parameter gradient in z-direction  d(ll)/dz		",
"									",
" nfile= 	name of input elastic param. (c66)  nn file, if no name ",
"		specified then, assume a linear profile with ...	",
" nn00=2.0	elastic parameter at (0, 0) 				",
" dnndx=0.0	parameter gradient in x-direction  d(nn)/dx		",
" dnndz=0.0	parameter gradient in z-direction  d(nn)/dz		",
"									",
" Optimizations:							",
" The moving boundary option permits the user to restrict the computations",
" of the wavefield to be confined to a specific range of spatial coordinates.",
" The boundary of this restricted area moves with the wavefield		",
" movebc=0	0 do not use moving boundary optimization		",
"		1 use moving boundaries					",
"  If movebc=1 then specify:						",
" mbx1=0 	initial left side of moving boundary			",
" mbz1=0 	initial top of moving boundary 				",
" mbx2=nx 	initial right side of moving boundary 			",
" mbz2=nz 	initial bottom of moving boundary 			",
"									",
NULL};

/*
 * Author: Tong Fei,	Center for Wave Phenomena, 
 *		Colorado School of Mines, Dec 1993
 * Some additional features by: Stig-Kyrre Foss, CWP
 *		Colorado School of Mines, Oct 2001
 * New features (Oct 2001): 
 * - setting receiver depth
 * - outputfiles with SU-headers
 * - additional commentary
 */

/* 
Notes:
	This program performs seismic modeling for elastic anisotropic 
	media with vertical axis of symmetry.  
	The finite-difference method with the FCT correction is used.

	Stability condition:	vmax*dt /(sqrt(2)*min(dx,dz)) < 1
	
	Two major stages are used in the algorithm:
	(1) conventional finite-difference wave extrapolation
	(2) followed by an FCT correction 

References:
	The detailed algorithm description is given in the article
	"Elimination of dispersion in finite-difference modeling 
	and migration"	in CWP-137, project review, page 155-174.

	Original reference to the FCT method:
	Boris, J., and Book, D., 1973, Flux-corrected transport. I.
	SHASTA, a fluid transport algorithm that works: 
	Journal of Computational Physics, vol. 11, p. 38-69.

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

/* function prototypes for subroutines used internally */
void 	read_parameter(int nx, int nz, float dx, float dz, 
		       float p00, float dpdx, float dpdz, char *file, float **pp);
void	anis_solver2(int it, float **u, float **v, float **w, 
		     float **e11, float **e33, float **e23, float **e12, float **e13,  
		     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);
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);
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);
void	boundary_vel(float **cc, float **rho, 
		     float *vell, float *velr, float *velt, float *velb, 
		     int nx, int nz, int mbx1, int mbz1, int mbx2, int mbz2);
void	locate_source(int nx, int nz, int sx, int sz, 
		      float **xzsource, int source);
void tforce_ricker(int n, float *tforce, float dt, float fpeak);
void tforce_akb(int n, float *tforce, float dt, float fpeak);
void tforce_spike(int n, float *tforce, float dt, float fpeak);
void tforce_unit(int n, float *tforce, float dt, float fpeak);
void	moving_bc (int it, int nx, int nz, int sx, int sz, 
		   float dx, float dz, int impulse, int movebc, 
		   float *t, float vmax, int *mbx1, int *mbz1, int *mbx2, int *mbz2);
void	moving_fctbc (int mbx1, int mbz1, int mbx2, int mbz2, 
		      int nxcc1, int nzcc1, int nxcc2, int nzcc2, 
		      int *fctxbeg, int *fctzbeg, int *fctxend, int *fctzend);
void	strain2_x(int mbx1, int mbx2, int mbz1, int mbz2, 
		  float dx, float **a, float **da);
void	strain2_z(int mbx1, int mbx2, int mbz1, int mbz2, 
		  float dz, float **a, float **da);
void	strain2_xy(int mbx1, int mbx2, int mbz1, int mbz2, 
		   float dx, float **ay, float **da);
void	strain2_xz(int mbx1, int mbx2, int mbz1, int mbz2, 
		   float dx, float dz, float **ax, float **az, float **da);
void	strain2_yz(int mbx1, int mbx2, int mbz1, int mbz2, 
		   float dz, float **ay, float **da);
void	strain4_x(int mbx1, int mbx2, int mbz1, int mbz2, 
		  float dx, float **a, float **da);
void	strain4_z(int mbx1, int mbx2, int mbz1, int mbz2, 
		  float dz, float **a, float **da);
void	strain4_xy(int mbx1, int mbx2, int mbz1, int mbz2, 
		   float dx, float **ay, float **da);
void	strain4_xz(int mbx1, int mbx2, int mbz1, int mbz2, 
		   float dx, float dz, float **ax, float **az, float **da);
void	strain4_yz(int mbx1, int mbx2, int mbz1, int mbz2, float dz, float **ay, float **da);
void	difference(int mbx1, int mbx2, int mbz1, int mbz2, int it, 
		   float **u1, float **u, 
		   float **xzsource, float *f, float dt, float rdxdz); 
void	difference_2x(int mbx1, int mbx2, int mbz1, int mbz2, int shift, 
		      float dtdz, float **u, float **e, float **c);
void	difference_2z(int mbx1, int mbx2, int mbz1, int mbz2, int shift, 
		      float dtdz, float **u, float **e, float **c);
void	difference_4x(int mbx1, int mbx2, int mbz1, int mbz2, int shift, 
		      float dtdz, float **u, float **e, float **c);
void	difference_4z(int mbx1, int mbx2, int mbz1, int mbz2, int shift, 
		      float dtdz, float **u, float **e, float **c);

void	fct2d1o(float **u, float **u1, int nx, int nz, 
		float eta0, float eta, float deta0dx, float deta0dz, 
		float detadx, float detadz, float dx, float dz, 
		int mbx1, int mbx2, int mbz1, int mbz2, 
		float **f0, float **f1);

void    su_output(int nx, int sx, int sdepth, int ns, float dx, 
		  float dt, FILE *outputfile, float **data);


int
main (int  argc,  char  **argv)
{

  int receiverdepth,suhead;
  int	ix, iz, it, nx, nz,  
    nt, sx, sz, vspnx;
  int	isurf, impulse, dofct, mbx1, mbx2, mbz1, mbz2;
  int 	fctxbeg, fctxend, fctzbeg, fctzend,
    nxcc1=0, nxcc2=0, nzcc1=0, nzcc2=0;
  int 	indexux, indexuy, indexuz, wavelet, movebc;
  int 	source, order; 
  float	dx, dz, dt, fpeak, rho00, vmax;
  float	eta, eta0;
  float	daadx, daadz, 
    dccdx, dccdz, dffdx, dffdz, dlldx, dlldz, 
    dnndx, dnndz, drhodx, drhodz; 
  float	deta0dx, deta0dz, detadx, detadz;
  float 	aa00, cc00, ff00, ll00, nn00;
  float	**u,  **v,  **w,   
    **e11, **e33, **e12, **e13, **e23, 
    **xzsource, **aa, **cc, **ff, **ll, **nn, **rho, 
    *fx, *fy, *fz, *t, 
    **refl_x=NULL, **refl_y=NULL, **refl_z=NULL,
    **vsp_x=NULL, **vsp_y=NULL, **vsp_z=NULL;
  FILE	*outfp=stdout;
  FILE	*outfpx, *outfpy, *outfpz;
  FILE	*outreflx, *outrefly, *outreflz,  
    *outvspx, *outvspy, *outvspz; 
  /* input file names */
  char *sfile="";		/* source file name */
  char *dfile="";		/* density file name */
  char *afile="";		/* file name for elastic parameter aa*/
  char *cfile="";		/* file name for elastic parameter cc*/
  char *ffile="";		/* file name for elastic parameter ff*/
  char *lfile="";		/* file name for elastic parameter ll*/
  char *nfile="";		/* file name for elastic parameter nn*/
  /* output file names */
  char *reflxfile="";	/* reflection seismogram file name, x-comp */
  char *reflyfile="";	/* reflection seismogram file name, y-comp */
  char *reflzfile="";	/* reflection seismogram file name, z-comp */
  char *vspxfile="";	/* VSP seismogram file name, x-comp */
  char *vspyfile="";	/* VSP seismogram file name, y-comp */
  char *vspzfile="";	/* VSP seismogram file name, z-comp */
  
  initargs (argc, argv);
  requestdoc(0);
  
  /* get required parameters  */
  if (!getparint("receiverdepth", &receiverdepth)) receiverdepth=1;
  if (!getparint("nt", &nt)) nt=250;
  if (!getparint("nx", &nx)) nx=300;
  if (!getparint("nz", &nz)) nz=250;
  if (!getparint("sx", &sx)) sx=nx/2;
  if (!getparint("sz", &sz)) sz=nz/2;
  if (!getparint("vspnx", &vspnx)) vspnx=nx/2;
  if (!getparint("impulse", &impulse)) impulse=0;
  if (!getparint("source", &source)) source=0;
  if (!getparint("isurf", &isurf)) isurf=1;
  if (!getparint("dofct", &dofct)) dofct=0;
  if (!getparint("fctxbeg", &fctxbeg)) fctxbeg=0;
  if (!getparint("fctzbeg", &fctzbeg)) fctzbeg=0;
  if (!getparint("fctxend", &fctxend)) fctxend=nx;
  if (!getparint("fctzend", &fctzend)) fctzend=nz;
  if (!getparint("indexux", &indexux)) indexux=0;
  if (!getparint("indexuy", &indexuy)) indexuy=0;
  if (!getparint("indexuz", &indexuz)) indexuz=0;
  if (!getparint("wavelet", &wavelet)) wavelet=1;
  if (!getparint("movebc", &movebc)) movebc=0;
  if (!getparint("order", &order)) order=2;
  if (!getparint("suhead",&suhead)) suhead=1;
  if (!getparfloat("dx", &dx)) dx=0.02;
  if (!getparfloat("dz", &dz)) dz=0.02;
  if (!getparfloat("dt", &dt)) dt=0.002;
  if (!getparfloat("fpeak", &fpeak)) fpeak=20.0;
  if (!getparfloat("aa00", &aa00)) aa00=2.0;
  if (!getparfloat("cc00", &cc00)) cc00=2.0;
  if (!getparfloat("ff00", &ff00)) ff00=2.0;
  if (!getparfloat("ll00", &ll00)) ll00=2.0;
  if (!getparfloat("nn00", &nn00)) nn00=2.0;
  if (!getparfloat("daadx", &daadx)) daadx=0.0;
  if (!getparfloat("daadz", &daadz)) daadz=0.0;
  if (!getparfloat("dccdx", &dccdx)) dccdx=0.0;
  if (!getparfloat("dccdz", &dccdz)) dccdz=0.0;
  if (!getparfloat("dffdx", &dffdx)) dffdx=0.0;
  if (!getparfloat("dffdz", &dffdz)) dffdz=0.0;
  if (!getparfloat("dlldx", &dlldx)) dlldx=0.0;
  if (!getparfloat("dlldz", &dlldz)) dlldz=0.0;
  if (!getparfloat("dnndx", &dnndx)) dnndx=0.0;
  if (!getparfloat("dnndz", &dnndz)) dnndz=0.0;
  if (!getparfloat("drhodx", &drhodx)) drhodx=0.0;
  if (!getparfloat("drhodz", &drhodz)) drhodz=0.0;
  if (!getparfloat("rho00", &rho00)) rho00=1.0;
  if (!getparfloat("eta0", &eta0)) eta0=0.03;
  if (!getparfloat("eta", &eta)) eta=0.04;
  if (!getparfloat("deta0dx", &deta0dx)) deta0dx=0.0;
  if (!getparfloat("deta0dz", &deta0dz)) deta0dz=0.0;
  if (!getparfloat("detadx", &detadx)) detadx=0.0;
  if (!getparfloat("detadz", &detadz)) detadz=0.0;
  getparstring("sfile",&sfile);
  getparstring("dfile",&dfile);
  getparstring("afile",&afile);
  getparstring("cfile",&cfile);
  getparstring("ffile",&ffile);
  getparstring("lfile",&lfile);
  getparstring("nfile",&nfile);
  getparstring("nfile",&nfile);
  getparstring("reflxfile",&reflxfile);
  getparstring("reflyfile",&reflyfile);
  getparstring("reflzfile",&reflzfile);
  getparstring("vspxfile",&vspxfile);
  getparstring("vspyfile",&vspyfile);
  getparstring("vspzfile",&vspzfile);
	
	
  /*   allocate space	*/
  u = alloc2float(nz, nx);
  v = alloc2float(nz,nx);
  w = alloc2float(nz,nx);
  e11 = alloc2float(nz,nx);
  e33 = alloc2float(nz,nx);
  e23 = alloc2float(nz,nx);
  e12 = alloc2float(nz,nx);
  e13 = alloc2float(nz,nx);
  xzsource = alloc2float(nz,nx);	
  aa = alloc2float(nz,nx);	
  cc = alloc2float(nz,nx);	
  ff = alloc2float(nz,nx);