最近在將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 frequencyOutput: 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=0verb=1 nt=1000 dt=0.001fmax=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

總結

以上是生活随笔為你收集整理的基于Madagascar的二维地震声波波动方程正演模拟的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。