基于Madagascar的二维地震声波波动方程正演模拟

最近在将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的二维地震声波波动方程正演模拟相关推荐

二阶声波正演c语言程序,声波波动方程正演模拟程序总结 - 图文
此处,为了成图完整,我用的是t2,而不是 t1,也就是把雷克子波向右移动了一段距离,使主要部分都显示出来.(频率采用的是30hz) 从图中可以看出程序是正确的,符合理论上雷克子波的波形. 第二部分:主 ...
声波正演c语言程序,二维频率域声波方程正演模拟
1. 概述频率域波场正演相对于时间域数值模拟来说,有其自身的优势.首先,在多炮数值模拟情况下,频率域相对于时间域效率更高,每个频率成分的阻抗矩阵只需要计算一次,加入并行计算后可以极大地提高计算效率: ...
基于Python的二维有限元声波方程正演计算
基于Python的二维有限元声波方程正演计算一.基础理论与相关公式的导出什么是有限元方法? 有限元是计算复杂数学问题近似解的工具.当数学方程过于复杂,无法用正常的方法求解,并且一定程度的误差是可以 ...
【电磁】基于matlab求解瞬变电磁TEM层状介质正演【含Matlab源码 2164期】
⛄一.获取代码方式获取代码方式1: 完整代码已上传我的资源:[电磁]基于matlab求解瞬变电磁TEM层状介质正演[含Matlab源码 2164期] 点击上面蓝色字体,直接付费下载,即可. 获取代码 ...
ParaView绘制gprMax正演模拟的波场快照方法（1）
ParaView绘制gprMax正演模拟的波场快照方法(1) gprMax是一款优秀的基于时域有限差分方法(FDTD)的电磁波数值模拟脚本软件,其正演模拟的结果通过波场快照的形式可以直观的显示出来,通 ...
gprMax 正演模拟中Ex、Ey、Ez三个分量之间的关系分析
gprMax 正演模拟中Ex.Ey.Ez三个分量之间的关系分析在 GPR 应用中,电场分量通常是测得量.我们一般的正演模拟用哪个电场分量呢文章目录 gprMax 正演模拟中Ex.Ey.Ez三个分量 ...
关于Gprmax正演模拟结果显示空白的原因分析
Gprmax正演模拟结果显示空白 gprMax正演模拟结果显示空白,与发射天线极化方向有关,电场方向应垂直剖面方向,与参数设置关系不是很大.二维 gprMax 显示空白与天线收发距有关. 文章目录 G ...
gprMax电磁波正演模拟方法
文章首发于:https://blog.zhaoxuan.site/archives/37.html: 第一时间获取最新文章请关注博客个人站:https://blog.zhaoxuan.site. 目录 ...
重力勘探正演模拟matlab,裴雪林, 郭万松 (1995) 高精度重力勘探技术在国内外的应用. 断块油田, 5, 8-11....
ABSTRACT: 本文典型地质模型三维重力正演模拟研究是在二维重力正演方法总结的基础上,对地质模型作三维重力正演探究.其中所建的地质模型以球体为主.分别导出了二维与三维重力布格异常,重力异常各阶导数 ...

基于Madagascar的二维地震声波波动方程正演模拟

基于Madagascar的二维地震声波波动方程正演模拟相关推荐

最新文章

热门文章