有限差分正演模拟

  • 0Z4wR7g38r
    了解作者
  • C/C++
    开发工具
  • 1.5MB
    文件大小
  • zip
    文件格式
  • 0
    收藏次数
  • 10 积分
    下载积分
  • 0
    下载次数
  • 2022-05-02 22:50
    上传日期
非常好用的有限差分正演程序,适用于新手学习使用。可以模拟地震波在介质中的传播规律,生成地震记录
有限差分程序.zip
  • Debug
  • 适合复杂的三角形网格剖分(直角三角形).obj
    21.8KB
  • 适合复杂的三角形网格剖分(直角三角形).exe
    228.1KB
  • vc60.pdb
    52KB
  • 适合复杂的三角形网格剖分(直角三角形)1.pch
    201.5KB
  • 适合复杂的三角形网格剖分(直角三角形).pch
    201.8KB
  • 适合复杂的三角形网格剖分(直角三角形) -1200.obj
    21.9KB
  • 适合复杂的三角形网格剖分(直角三角形) -1200.pdb
    497KB
  • 适合复杂的三角形网格剖分(直角三角形) -1200.pch
    201.5KB
  • file.pch
    228KB
  • vc60.idb
    65KB
  • 适合复杂的三角形网格剖分(直角三角形).pdb
    497KB
  • 适合复杂的三角形网格剖分(直角三角形)1.pdb
    497KB
  • 适合复杂的三角形网格剖分(直角三角形) -2400.pch
    201.8KB
  • 适合复杂的三角形网格剖分(等边三角形).ilk
    228.7KB
  • file.obj
    14.9KB
  • 适合复杂的三角形网格剖分(等边三角形).pch
    201.7KB
  • 适合复杂的三角形网格剖分(直角三角形) -2400.ilk
    227.5KB
  • file.exe
    220KB
  • 适合复杂的三角形网格剖分(等边三角形).obj
    21KB
  • 适合复杂的三角形网格剖分(直角三角形) -1200.ilk
    278KB
  • 适合复杂的三角形网格剖分(直角三角形) -2400.exe
    228.1KB
  • file.pdb
    497KB
  • 适合复杂的三角形网格剖分(直角三角形)1.exe
    228.1KB
  • 适合复杂的三角形网格剖分(等边三角形).exe
    228.1KB
  • 适合复杂的三角形网格剖分(直角三角形).ilk
    224.8KB
  • file.ilk
    260.6KB
  • 适合复杂的三角形网格剖分(直角三角形)1.ilk
    229.2KB
  • 适合复杂的三角形网格剖分(等边三角形).pdb
    497KB
  • 适合复杂的三角形网格剖分(直角三角形)1.obj
    23.8KB
  • 适合复杂的三角形网格剖分(直角三角形) -2400.obj
    21.9KB
  • 适合复杂的三角形网格剖分(直角三角形) -1200.exe
    228.1KB
  • 适合复杂的三角形网格剖分(直角三角形) -2400.pdb
    497KB
  • file.cpp
    6.6KB
  • int_alloc.h
    939B
内容介绍
#include<stdio.h> #include<math.h> #include <stdlib.h> #include <time.h> #define T_Number 1301 /////波场快照时间的数值 #define L_Number 601 /////横向网格数 #define D_Number 601 /////纵向网格数 #define PI 3.141593 #define PML 40 /////////////给动态二维数组分配空间//////////// float **dmatrix(int a,int b) { float **v; int i; v=(float **)malloc(a*sizeof(float *)); for(i=0;i<a;i++) { v[i]=(float *)malloc(b*sizeof(float *)); } return v; } //////////////初始化处理///////////////// void zero(int a,int b,float **zero) { int i,j; for(i=0;i<a;i++) { for(j=0;j<b;j++) { zero[i][j]=0.0; } } } /////求横向空间导数/////////////////// void spatial_derivative_xx(float **Dx,float **U) { float a0,a1,a2,a3,a4,a5; int i,j; a0=-2.92722; a1=1.66667; a2=-0.238095; a3=0.0396825; a4=-0.00496032; a5=0.000317460; for(i=5;i<L_Number-5;i++) { for(j=5;j<D_Number-5;j++) { Dx[i][j]=2.0*(a0*U[i][j]+a1*(U[i+1][j]+U[i-1][j])+ a2*(U[i+2][j]+U[i-2][j])+a3*(U[i+3][j] +U[i-3][j])+a4*(U[i+4][j]+U[i-4][j]) +a5*(U[i+5][j]+U[i-5][j])); } } } /////求横向空间导数/////////////////// void spatial_derivative_zz(float **Dz,float **U) { float a0,a1,a2,a3,a4,a5; int i,j; a0=-2.92722; a1=1.66667; a2=-0.238095; a3=0.0396825; a4=-0.00496032; a5=0.000317460; for(i=5;i<L_Number-5;i++) { for(j=5;j<D_Number-5;j++) { Dz[i][j]=2.0*(a0*U[i][j]+a1*(U[i][j+1]+U[i][j-1]) +a2*(U[i][j+2]+U[i][j-2])+a3*(U[i][j+3] +U[i][j-3])+a4*(U[i][j+4]+U[i][j-4]) +a5*(U[i][j+5]+U[i][j-5])); } } } ///////////////////////读取速度///////////////////////// void read_file(int vnx,int vnz,float**vp0) { FILE *fp=NULL; int i,j; fp=fopen("Velocity_model.bin","rb"); if(fp==NULL) printf("NULL\n"); else{ for(i=0;i<vnx;i++) for(j=0;j<vnz;j++) { fread(&vp0[i][j],sizeof(float),1,fp); } } fclose(fp); } void absortion(float **vx) { int i,j; float qp=0.01; /////上边界///// for(i=PML;i<L_Number-PML;i++) for(j=0;j<PML;j++) { vx[i][j]=exp(-qp*(PML-j)*qp*(PML-j))*vx[i][j]; } ////下边界//// for(i=PML;i<L_Number-PML;i++) for(j=D_Number-PML;j<D_Number;j++) { vx[i][j]=exp(-qp*(j-D_Number+PML)*qp*(j-D_Number+PML))*vx[i][j]; } /////左边界///// for(i=0;i<PML;i++) for(j=0;j<D_Number;j++) { vx[i][j]=exp(-qp*(PML-i)*qp*(PML-i))*vx[i][j]; } ///////右边界//////// for(i=L_Number-PML;i<L_Number;i++) for(j=0;j<D_Number;j++) { vx[i][j]=exp(-qp*(i-L_Number+PML)*qp*(i-L_Number+PML))*vx[i][j]; } } //////////////主函数///////////////// void main() { clock_t start,finish; double total_time; FILE *fp=NULL,*fp_temp=NULL,*fp_temp1=NULL; float dt,dx,dz,depth,length,time,f0,vp,R; int i,j,k,m,n;////循环变量 int izs,ixs; int node_number; float f[T_Number]; int temp1,temp2; FILE *fp5; fp5=fopen("Velocity_FD_model.bin","wb"); float **U,**U0,**U1,**uxx,**uzz,**alphy,**p,**vp0,**ddx,**ddz; start = clock(); ////////采样间隔////////// dt=0.0002; dx=2.0; dz=2.0; f0=40.0; ///////////震源位置点///////// izs= floor( L_Number/2 ); ixs= floor (D_Number/2 ); U0 = (float**)dmatrix(L_Number,D_Number); U1 = (float**)dmatrix(L_Number,D_Number); U = (float**)dmatrix(L_Number,D_Number); alphy =(float**)dmatrix(L_Number,D_Number); uxx = (float**)dmatrix(L_Number,D_Number); uzz = (float**)dmatrix(L_Number,D_Number); p = (float**)dmatrix(L_Number,T_Number); vp0 = (float**)dmatrix(L_Number,D_Number); ddx = (float**)dmatrix(L_Number,D_Number); ddz = (float**)dmatrix(L_Number,D_Number); vp0 = (float**)dmatrix(L_Number,D_Number); zero(L_Number,D_Number,vp0); zero(L_Number,D_Number,uxx); zero(L_Number,T_Number,p); zero(L_Number,D_Number,uzz); zero(L_Number,D_Number,ddx); zero(L_Number,D_Number,ddz); zero(L_Number,D_Number,U0); zero(L_Number,D_Number,U1); zero(L_Number,D_Number,U); zero(L_Number,D_Number,alphy); zero(L_Number,D_Number,vp0); alphy[ixs][izs]=1.0; ////////////点震源处理///////////// fp=fopen("wavelet.xls","w"); for (k=0;k<T_Number;k++) { f[k]=(1-2*PI*f0*(dt*k-1.0/f0)*PI*f0*(dt*k-1.0/f0))*exp(-PI*f0*(dt*k-1.0/f0)*PI*f0*(dt*k-1.0/f0)); fprintf(fp,"%f\t%f\n",k*dt,f[k]); } fclose(fp); /////////////////////////////////////////// // read_file(L_Number,D_Number,vp0); for(j=0;j<=300;j++) for(i=0;i<L_Number;i++) { vp0[i][j]=2000; } for(j=301;j<D_Number;j++) for(i=0;i<L_Number;i++) { vp0[i][j]=2500; } temp1=100; temp2=500; for(j=301;j<=400;j++) { for(i=temp1;i<=temp2;i++) { vp0[i][j]=2000; } temp1=temp1+1; temp2=temp2-1; } printf("%f,%f,%f", vp0[1][0],vp0[1][3],vp0[2][2]); for(i=0;i<L_Number;i++) for(j=0;j<D_Number;j++) { fwrite(&vp0[i][j],sizeof(float),1,fp5); } for (k=0;k<T_Number;k++) { printf("k=%d\n",k); spatial_derivative_xx( uxx, U); spatial_derivative_zz( uzz, U); ////////////////应力计算/////////// for(i=0;i<L_Number;i++) { for(j=0;j<D_Number;j++) { U[i][j]=2.0*U1[i][j]-U0[i][j]+0.5*pow((vp0[i][j]*dt)/dx,2)*uxx[i][j]+0.5*pow((vp0[i][j]*dt)/dz,2)*uzz[i][j]+f[k]*alphy[i][j]; } } for(i=0;i<L_Number;i++) { p[i][k]=U[i][izs]; } absortion(U1);absortion(U);absortion(U0); ////////////////迭代替换/////////// for(i=0;i<L_Number;i++) { for(j=0;j<D_Number;j++) { float temp; temp=U1[i][j]; U1[i][j]=U[i][j]; U0[i][j]=temp; } } } ///////////////地震记录计算////////////// fp_temp1=fopen("record1.dat","wb"); for(i=0;i<L_Number;i++) { for(k=0;k<T_Number;k++) { fwrite(&p[i][k],sizeof(float),1,fp_temp1); } } fclose(fp_temp1); ////////////////声波应力场输出/////////// // fp_temp=fopen("Snapshot3.dat","wb"); for(i=0;i<L_Number;i++) { for(j=0;j<D_Number;j++) { fwrite(&U[i][j],sizeof(float),1,fp_temp); } } fclose(fp_temp); finish=clock(); printf("the start time is %f,the finish time is %f\n",(double)start/CLOCKS_PER_SEC,(double)finish/CLOCKS_PER_SEC); total_time=(double)(finish-start)/CLOCKS_PER_SEC; printf("the computed total time is %f seconds\n",total_time); }
评论
    相关推荐