最近在将SU写的地震勘探的程序迁移到Madagascar上，初步尝试，写了一个二维声波方程正演程序，很简单，也很基本，只能输出波场快照，没有吸收边界条件，贴出来，供大家参考。代码和脚本如下：
#include <time.h>
#include "rsf.h"
#define FSIZE sizeof(float)
static float ricker (float t, float fpeak);
void ptsrc (float saf,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 fd2d (float multis,float xs,float zs,int nx,float dx,float fx,int nz,
float dz,float fz,int nt,float dt,float fmax,float fpeak,float tdelay,
float **s,float **v,float **pf,float **p,float **pb);
int main (int argc, char *argv[])
{
bool verb; /* verbose flag */
int T_beg,T_end,T_dur; /* program run time */
int it,iz,ix; /* index variables */
int nz,nx;
float dz,dx;
float fz,fx;
float h;
int nt;
float dt;
float tmax;
float mt;
float t;
int ns;
int is;
float dxs,dzs;
float fxs,fzs;
float *xs,*zs;
float c0,c1,c2; /* Laplacian coefficients */
float vmin,vmax;
float **v;
float fpeak,fmax;
float saf;
float tdelay;
sf_axis at,az,ax; /* cube axes */
FILE *fv=NULL; /* velocitu file */
sf_file fo=NULL; /* output file */
float **s;
float **Pf;
float **P;
float **Pb;
float **Ptmpt;
sf_init(argc,argv);
T_beg=clock();
sf_warning("************ Program BEG! ***********.\n");
/* get required parameters */
if(! sf_getint("nt",&nt)) sf_error(" must specify nt!");
if(! sf_getint("ns",&ns)) sf_error(" must specify ns!");
if(! sf_getint("nx",&nx)) sf_error(" must specify nx!");
if(! sf_getint("nz",&nz)) sf_error(" must specify nz!");
if(! sf_getfloat("dt",&dt)) sf_error(" must specify dt!");
if(! sf_getfloat("dx",&dx)) sf_error(" must specify dx!");
if(! sf_getfloat("dz",&dz)) sf_error(" must specify dz!");
if(! sf_getfloat("dxs",&dxs)) sf_error(" must specify dxs!");
if(! sf_getfloat("dzs",&dzs)) sf_error(" must specify dzs!");
/* Input and output file information */
if(! sf_getbool("verb",&verb)) verb=0;
if(! sf_getfloat("fx",&fx)) fx=0;
if(! sf_getfloat("fz",&fz)) fz=0;
if(! sf_getfloat("fxs",&fxs)) fxs=0;
if(! sf_getfloat("fzs",&fzs)) fzs=0;
if(! sf_getfloat("saf",&saf)) saf=1;
if(! sf_getfloat("tdelay",&tdelay)) tdelay=0.0;
sf_warning("verb=%d\n",verb);
sf_warning("nt=%d,nx=%d,nz=%d\n",nt,nx,nz);
sf_warning("dt=%f,dx=%f,dz=%f\n",dt,dx,dz);
/* allocate wavefield arrays */
xs = sf_floatalloc(ns);
zs = sf_floatalloc(ns);
v = sf_floatalloc2(nz,nx);
s = sf_floatalloc2(nz,nx);
Pf = sf_floatalloc2(nz,nx);
P = sf_floatalloc2(nz,nx);
Pb = sf_floatalloc2(nz,nx);
memset((void *) v[0], 0,FSIZE*nz*nx);
memset((void *) s[0], 0,FSIZE*nz*nx);
memset((void *) Pf[0],0,FSIZE*nz*nx);
memset((void *) P[0], 0,FSIZE*nz*nx);
memset((void *) Pb[0],0,FSIZE*nz*nx);
/* read velocity */
fv=stdin;
fread(v[0],sizeof(float),nx*nz,fv);
sf_warning("*****v[%d][%d]=%f s*****.\n",1,1,v[1][1]);
/* determine minimum and maximum velocities */
vmin = vmax = v[0][0];
for (ix=0; ix<nx;ix++)
{
for (iz=0; iz<nz;iz++)
{
vmin = SF_MIN(vmin,v[ix][iz]);
vmax = SF_MAX(vmax,v[ix][iz]);
}
}
sf_warning("vmax=%f;vmin=%f\n",vmax,vmin);
/* determine mininum spatial sampling interval */
h = SF_MIN(SF_ABS(dx),SF_ABS(dz));
/* determine time sampling interval to ensure stability */
if (dt > h/(2.0*vmax))
{
sf_error(" dt must <= %f!",h/(2.0*vmax));
}
/* determine maximum temporal frequency to avoid dispersion */
if(! sf_getfloat("fmax",&fmax)) sf_error(" must specify fmax!");
if (fmax > vmin/(10.0*h)) sf_error(" fmax must <= %f!",vmin/(10.0*h));
/* compute or set peak frequency for ricker wavelet */
if(! sf_getfloat("fpeak",&fpeak)) sf_error(" must specify fpeak!");
if (SF_NINT(fmax/fpeak) != 2) sf_error(" fpeak must = fmax/2!");
/* determine source coordinates */
for (is=0;is<ns;is++)
{
xs[is] = fxs+dxs*is;
zs[is] = fzs+dzs*is;
// sf_warning("xs=%f;zs=%f",xs[is],zs[is]);
}
/* do finite difference modeling */
sf_warning(" ****************** do FM ***********************");
for ( is = 0; is < ns; ++is)
{
sf_warning(" Forward Progress source:%d/%d",is+1,ns);
fd2d (saf,xs[is],zs[is],nx,dx,fx,nz,dz,fz,nt,dt,fmax,fpeak,tdelay,s,
v,Pf,P,Pb);
}
sf_close();
T_end = clock();
T_dur = T_end - T_beg ;
sf_warning("*****Program END! It takes %f s*****.",T_dur/1000.0);
exit (0);
}
static float ricker (float t, float fpeak)
/*****************************************************************************
Compute Ricker wavelet as a function of time
******************************************************************************
Input:
t time at which to evaluate Ricker wavelet
fpeak peak (dominant) frequency of wavelet
******************************************************************************
Notes:
The amplitude of the Ricker wavelet at a frequency of 2.5*fpeak is
approximately 4 percent of that at the dominant frequency fpeak.
The Ricker wavelet effectively begins at time t = -1.0/fpeak. Therefore,
for practical purposes, a causal wavelet may be obtained by a time delay
of 1.0/fpeak.
The Ricker wavelet has the shape of the second derivative of a Gaussian.
******************************************************************************
Author: Dave Hale, Colorado School of Mines, 04/29/90
******************************************************************************/
{
float x,xx;
x = SF_PI*fpeak*t;
xx = x*x;
/* return (-6.0+24.0*xx-8.0*xx*xx)*exp(-xx); */
/* return SF_PI*fpeak*(4.0*xx*x-6.0*x)*exp(-xx); */
return exp(-xx)*(1.0-2.0*xx);
}
void ptsrc (float saf,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)
/*******************************************************************************
update source pressure function for a point source
********************************************************************************
Input:
xs x coordinate of point source
zs z coordinate of point source
nx number of x samples
dx x sampling interval
fx first x sample
nz number of z samples
dz z sampling interval
fz first z sample
dt time step (ignored)
t time at which to compute source function
fmax maximum frequency (ignored)
fpeak peak frequency
Output:
tdelay time delay of beginning of source function
s array[nx][nz] of source pressure at time t+dt
********************************************************************************
Author: Dave Hale, Colorado School of Mines, 03/01/90
*******************************************************************************/
{
int ix,iz,ixs,izs;
float ts,xn,zn,xsn,zsn;
memset((void *)s[0], (int)'\0', nx*nz*FSIZE);
/* compute time-dependent part of source function */
/* fpeak = 0.5*fmax; this is now getparred */
tdelay = 1.0/fpeak;
if (t>2.0*tdelay) return;
ts = ricker(t-tdelay,fpeak);
/* let source contribute within limited distance */
xsn = (xs-fx)/dx;
zsn = (zs-fz)/dz;
ixs = SF_NINT(xsn);
izs = SF_NINT(zsn);
for (ix=SF_MAX(0,ixs-3); ix<=SF_MIN(nx-1,ixs+3); ++ix)
{
for (iz=SF_MAX(0,izs-3); iz<=SF_MIN(nz-1,izs+3); ++iz)
{
xn = ix-xsn;
zn = iz-zsn;
s[ix][iz] = saf*ts*exp(-xn*xn-zn*zn);
}
}
}
void fd2d (float multis,float xs,float zs,int nx,float dx,float fx,
int nz,float dz,float fz,
int nt,float dt, float fmax, float fpeak, float tdelay, float **s,
float **v,float **pf,float **p,float **pb)
/*******************************************************************************
update source pressure function for a point source
********************************************************************************/
{
int ix,iz,it;
int N;
float a0,a1,a2,a3;
float t;
float **ptemp;
char fname[BUFSIZ];
FILE *fp = NULL;
a0 = -2.7222;
a1 = 1.5000;
a2 = -0.1500;
a3 = 1/90;
N=6;
for ( it = 0,t=0; it < nt; ++it,t+=dt)
{
ptsrc (multis,xs,zs,nx,dx,fx,nz,dz,fz,dt,t,fmax,fpeak,tdelay,s);
for ( ix = N/2; ix < nx-N/2; ++ix)
{
for ( iz = N/2; iz < nz-N/2; ++iz)
{
pf[ix][iz]=2*p[ix][iz]-pb[ix][iz] +
v[ix][iz]*v[ix][iz]*(dt*dt)*
( (a0*p[ix][iz]+
a1*(p[ix+1][iz]+p[ix-1][iz])+
a2*(p[ix+2][iz]+p[ix-2][iz])+
a3*(p[ix+3][iz]+p[ix-3][iz])
)/dx/dx/2+
(a0*p[ix][iz]+
a1*(p[ix][iz+1]+p[ix][iz-1])+
a2*(p[ix][iz+2]+p[ix][iz-2])+
a3*(p[ix][iz+3]+p[ix][iz-3])
)/dz/dz/2
)+
s[ix][iz]*v[ix][iz]*v[ix][iz]*(dt*dt);
}
}
sprintf(fname,"fw_%d.bin",it);
fp=fopen(fname,"wb");
fwrite(p[0],sizeof(float),nx*nz,fp);
fclose(fp);
ptemp=pb;
pb=p;
p=pf;
pf=ptemp;
}
}
声明：本程序借鉴了SU的部分代码，本程序仅限于学习，如有他用，后果自负。
程序运行脚本：
#!/bin/sh
#
vel=Data/model/vel.bin
# model information
n1=201 d1=10 f1=0.0 label1="Depth (km)"
n2=401 d2=10 f2=0.0 label2="Distance (km)"
# seismic source information
ns=1 fxs=2000 fzs=1000 dxs=50 dzs=0
verb=1
nt=1000 dt=0.001
fmax=40 fpeak=20 saf=100 tdelay=0.1
./SRC/fdm2d_cpu < $vel \
verb=$verb nt=$nt dt=$dt \
nx=$n2 dx=$d2 fx=$f2 nz=$n1 dz=$d1 fz=$f1 \
fmax=$fmax fpeak=$fpeak saf=$saf tdelay=$tdelay \
ns=$ns fxs=$fxs fzs=$fzs dxs=$dxs dzs=$dzs \
exit 0
![查看Linux系统文件夹下各文件所占空间大小[图]](data:image/gif;base64,R0lGODlhAQABAIAAAP///wAAACwAAAAAAQABAAACAkQBADs=)