sjxishou1.cpp

二维声波方程正演

  //有限差分时间域二阶、空间域四阶二维声波方程加吸收边界常速模拟研究 (F=30  C=2500  时好)   震源不可放在边界
  //
  //      u(i,k,n+1)=(p*p/12)*{16[u(i+1,k,n)+u(i-1,k,n)+u(i,k+1,n)+u(i,k-1,n)]-[u(i-2,k,n)+u(i+2,k,n)+u(i,k+2,n)+u(i,k-2,n)]}
  //                 +(2-5*p*p)*u(i,k,n)-u(i,k,n-1)+C*C*DT*DT*S(m,n,k)
  //
  //       p=C*DT/h
  //
  //  //eponge边界吸收函数
  //    xx=0.305/IABMAX;
  //    for(i=0;i<=IABMAX-1;i++)
  //      {  aaa=pow((IABMAX-i),2);
  //         eponge[i]=exp(-xx*aaa);
  //      }
  //
  // 不要求:  DX=DZ      25hz ,吸收条件好 
  //            　　                  

         
#include <stdio.h>
#include <math.h>
#include <iostream.h>
#include <iomanip.h>       //输入输出格式控制
#include <fstream.h>
#define F 31.0            //雷克子波主频，单位:hz 
#define N0 2            //震源起爆时间 
#define NZ 91             //纵向最大采样点号 
#define NX 91             //横向最大采样点号 
#define DT 0.002          //时间采样间隔，单位：s 
#define NT 100             //时间最大样点号 
#define DZ 9             //纵向采样间隔 
#define DX 9             //横向采样间隔 
#define XS 45             //震源位置横坐标样点号 
#define ZS 45             //震源位置纵坐标样点号 
#define C 2500            //波传播速度，单位：m/s 
#define T1 49          //输出t=T1*DT时的切片 
#define IABMAX 20           //边界吸收宽度 

double u0[NX][NZ],u1[NX][NZ],u2[NX][NZ];
double wavelet(double f,int t0,int t,double dt)
{double wav,a,b;
 b=dt*(t-t0);
 a=pow(3.1415,2)*pow(f,2)*pow(b,2)*(-1);
 wav=(1+2*a)*exp(a);
 return wav;
 }
double shot(int a,int b,int e)
{double m,n,mn,s[NT];
  if(a==XS&&b==ZS)
   { m=pow(C,2)*pow(DT,2);
     n=DX*DZ;
     s[e]=wavelet(F,N0,e,DT);
     mn=m*s[e]/n;
    }
  else mn=0.0;
 return mn;
}
void printu2()
{
   int i,k;
   ofstream outfile("sjxishou.out");
   
   outfile<<setiosflags(ios::left)<<setiosflags(ios::scientific);
   for(i=IABMAX;i<=NX-IABMAX-1;i++)
   {   for(k=IABMAX;k<=NZ-IABMAX-1;k++)
        outfile<<setw(20)<<u2[i][k]<<endl;
        
   }
    cout<<"write the file success!"<<"\n";
    outfile.close();     // 调用成员函数close 
}
void main()
{  int i,k,n,g;
   double rx,rz,a,b,w;
   double eponge[IABMAX];
   double aaa,xx;
   //第一步：速度文件的载入及相关整理（替换）  
   rx=DT/DX;
   rz=DT/DZ;
   a=pow(C,2)*pow(rx,2);
   b=pow(C,2)*pow(rz,2);
   //eponge边界吸收函数
   xx=0.305/IABMAX;
   for(i=0;i<=IABMAX-1;i++)
      {  aaa=pow((IABMAX-i),2);
         eponge[i]=exp(-xx*aaa)*1e-10;
      }

   for(n=0;n<=NT-1;n++)
     {  //第二步：赋初值;初始时刻的全波场值均为零，P(i, j, 0)=0） 
         //  时刻dt时，在炮点S (XS, ZS)给出一个脉冲震源S(t)，其它点波场P =0； 

       if(n==0)
          {  for(i=0;i<=NX-1;i++)
             for(k=0;k<=NZ-1;k++)
                u0[i][k]=0.0;

           }
       else if(n==1)
          {  for(i=0;i<=NX-1;i++)
             for(k=0;k<=NZ-1;k++)
                u1[i][k]=0.0;

           }
       //第三步：边界条件处理及11点差分计算波场延拓  

       else if(n>=2)
           {  g=n-1;
               for(i=2;i<=NX-3;i++)
               for(k=2;k<=NZ-3;k++)
                 {  w=(2-2.5*a-2.5*b)*u1[i][k]-u0[i][k];
                    if(i==XS&&k==ZS)
                       u2[i][k]=a/12*(16*(u1[i-1][k]+u1[i+1][k])-(u1[i-2][k]+u1[i+2][k]))+b/12*(16*(u1[i][k-1]+u1[i][k+1])-(u1[i][k-2]+u1[i][k+2]))+w+shot(i,k,g);
                    else
                       u2[i][k]=a/12*(16*(u1[i-1][k]+u1[i+1][k])-(u1[i-2][k]+u1[i+2][k]))+b/12*(16*(u1[i][k-1]+u1[i][k+1])-(u1[i][k-2]+u1[i][k+2]))+w;
                 }
               u2[0][0]=a/12*(16*u1[1][0]-u1[2][0])+b/12*(16*u1[0][1]-u1[0][2])+(2-2.5*a-2.5*b)*u1[0][0]-u0[0][0];

               u2[1][0]=a/12*(16*(u1[0][0]+u1[2][0])-u1[3][0])+b/12*(16*u1[1][1]-u1[1][2])+(2-2.5*a-2.5*b)*u1[1][0]-u0[1][0];

               u2[0][1]=a/12*(16*u1[1][1]-u1[2][1])+b/12*(16*(u1[0][0]+u1[0][2])-u1[0][3])+(2-2.5*a-2.5*b)*u1[0][1]-u0[0][1];

               u2[1][1]=a/12*(16*(u1[0][1]+u1[2][1])-u1[3][1])+b/12*(16*(u1[1][0]+u1[1][2])-u1[1][3])+(2-2.5*a-2.5*b)*u1[1][1]-u0[1][1];

               u2[NX-1][0]=a/12*(16*u1[NX-2][0]-u1[NX-3][0])+b/12*(16*u1[NX-1][1]-u1[NX-1][2])+(2-2.5*a-2.5*b)*u1[NX-1][0]-u0[NX-1][0];

               u2[NX-2][0]=a/12*(16*(u1[NX-3][0]+u1[NX-1][0])-u1[NX-4][0])+b/12*(16*u1[NX-2][1]-u1[NX-2][2])+(2-2.5*a-2.5*b)*u1[NX-2][0]-u0[NX-2][0];

               u2[NX-1][1]=a/12*(16*u1[NX-2][1]-u1[NX-3][1])+b/12*(16*(u1[NX-1][0]+u1[NX-1][2])-u1[NX-1][3])+(2-2.5*a-2.5*b)*u1[NX-1][1]-u0[NX-1][1];

               u2[NX-2][1]=a/12*(16*(u1[NX-3][1]+u1[NX-1][1])-u1[NX-4][1])+b/12*(16*(u1[NX-2][0]+u1[NX-2][2])-u1[NX-2][3])+(2-2.5*a-2.5*b)*u1[NX-2][1]-u0[NX-2][1];

               u2[0][NZ-1]=a/12*(16*u1[1][NZ-1]-u1[2][NZ-1])+b/12*(16*u1[0][NZ-2]-u1[0][NZ-3])+(2-2.5*a-2.5*b)*u1[0][NZ-1]-u0[0][NZ-1];

               u2[0][NZ-2]=a/12*(16*u1[1][NZ-2]-u1[2][NZ-2])+b/12*(16*(u1[0][NZ-3]+u1[0][NZ-1])-u1[0][NZ-4])+(2-2.5*a-2.5*b)*u1[0][NZ-2]-u0[0][NZ-2];

               u2[1][NZ-1]=a/12*(16*(u1[0][NZ-1]+u1[2][NZ-1])-u1[3][NZ-1])+b/12*(16*u1[1][NZ-2]-u1[1][NZ-3])+(2-2.5*a-2.5*b)*u1[1][NZ-1]-u0[1][NZ-1];

               u2[1][NZ-2]=a/12*(16*(u1[0][NZ-2]+u1[2][NZ-2])-u1[3][NZ-2])+b/12*(16*(u1[1][NZ-3]+u1[1][NZ-1])-u1[1][NZ-4])+(2-2.5*a-2.5*b)*u1[1][NZ-2]-u0[1][NZ-2];

               u2[NX-1][NZ-1]=a/12*(16*u1[NX-2][NZ-1]-u1[NX-3][NZ-1])+b/12*(16*u1[NX-1][NZ-2]-u1[NX-1][NZ-3])+(2-2.5*a-2.5*b)*u1[NX-1][NZ-1]-u0[NX-1][NZ-1];

               u2[NX-2][NZ-1]=a/12*(16*(u1[NX-3][NZ-1]+u1[NX-1][NZ-1])-u1[NX-4][NZ-1])+b/12*(16*u1[NX-2][NZ-2]-u1[NX-2][NZ-3])+(2-2.5*a-2.5*b)*u1[NX-2][NZ-1]-u0[NX-2][NZ-1];

               u2[NX-1][NZ-2]=a/12*(16*u1[NX-2][NZ-2]-u1[NX-3][NZ-2])+b/12*(16*(u1[NX-1][NZ-3]+u1[NX-1][NZ-1])-u1[NX-1][NZ-4])+(2-2.5*a-2.5*b)*u1[NX-1][NZ-2]-u0[NX-1][NZ-2];

               u2[NX-2][NZ-2]=a/12*(16*(u1[NX-3][NZ-2]+u1[NX-1][NZ-2])-u1[NX-4][NZ-2])+b/12*(16*(u1[NX-2][NZ-3]+u1[NX-2][NZ-1])-u1[NX-2][NZ-4])+(2-2.5*a-2.5*b)*u1[NX-2][NZ-2]-u0[NX-2][NZ-2];

               for(i=2;i<=NX-3;i++)
                  u2[i][0]=a/12*(16*(u1[i-1][0]+u1[i+1][0])-(u1[i-2][0]+u1[i+2][0]))+b/12*(16*u1[i][1]-u1[i][2])+(2-2.5*a-2.5*b)*u1[i][0]-u0[i][0];

               for(i=2;i<=NX-3;i++)
                  u2[i][NZ-1]=a/12*(16*(u1[i-1][NZ-1]+u1[i+1][NZ-1])-(u1[i-2][NZ-1]+u1[i+2][NZ-1]))+b/12*(16*u1[i][NZ-2]-u1[i][NZ-3])+(2-2.5*a-2.5*b)*u1[i][NZ-1]-u0[i][NZ-1];

               for(k=2;k<=NZ-3;k++)
                  u2[0][k]=a/12*(16*u1[1][k]-u1[2][k])+b/12*(16*(u1[0][k-1]+u1[0][k+1])-(u1[0][k-2]+u1[0][k+2]))+(2-2.5*a-2.5*b)*u1[0][k]-u0[0][k];

               for(k=2;k<=NZ-3;k++)
                  u2[NX-1][k]=a/12*(16*u1[NX-2][k]-u1[NX-3][k])+b/12*(16*(u1[NX-1][k-1]+u1[NX-1][k+1])-(u1[NX-1][k-2]+u1[NX-1][k+2]))+(2-2.5*a-2.5*b)*u1[NX-1][k]-u0[NX-1][k];

               for(i=2;i<=NX-3;i++)
                  u2[i][1]=a/12*(16*(u1[i-1][1]+u1[i+1][1])-(u1[i-2][1]+u1[i+2][1]))+b/12*(16*(u1[i][0]+u1[i][2])-u1[i][3])+(2-2.5*a-2.5*b)*u1[i][1]-u0[i][1];

               for(i=2;i<=NX-3;i++)
                  u2[i][NZ-2]=a/12*(16*(u1[i-1][NZ-2]+u1[i+1][NZ-2])-(u1[i-2][NZ-2]+u1[i+2][NZ-2]))+b/12*(16*(u1[i][NZ-3]+u1[i][NZ-1])-u1[i][NZ-4])+(2-2.5*a-2.5*b)*u1[i][NZ-2]-u0[i][NZ-2];

               for(k=2;k<=NZ-3;k++)
                  u2[1][k]=a/12*(16*(u1[0][k]+u1[2][k])-u1[3][k])+b/12*(16*(u1[1][k-1]+u1[1][k+1])-(u1[1][k-2]+u1[1][k+2]))+(2-2.5*a-2.5*b)*u1[1][k]-u0[1][k];

               for(k=2;k<=NZ-3;k++)
                  u2[NX-2][k]=a/12*(16*(u1[NX-3][k]+u1[NX-1][k])-u1[NX-4][k])+b/12*(16*(u1[NX-2][k-1]+u1[NX-2][k+1])-(u1[NX-2][k-2]+u1[NX-2][k+2]))+(2-2.5*a-2.5*b)*u1[NX-2][k]-u0[NX-2][k];
               //左右边界吸收 
               for(k=0;k<=NZ-1;k++)
               for(i=0;i<=IABMAX-1;i++)
			   {  
				  u2[i][k]=u2[i][k]*eponge[i];
                  u2[i][k]=u2[i][k]*eponge[i];
                  u2[NX-1-i][k]=u2[NX-1-i][k]*eponge[i];
                  u2[NX-1-i][k]=u2[NX-1-i][k]*eponge[i];
			   }
               //上下边界吸收 
               for(k=0;k<=IABMAX-1;k++)
               for(i=0;i<=NX-1;i++)
               {  
                  u2[i][k]=u2[i][k]*eponge[k];
                  u2[i][k]=u2[i][k]*eponge[k];
                  u2[i][NZ-1-k]=u2[i][NZ-1-k]*eponge[k];
                  u2[i][NZ-1-k]=u2[i][NZ-1-k]*eponge[k];
               }
               if(n==T1)
                  printu2();
               // 数组循环覆盖  
               for(i=0;i<=NX-1;i++)
               for(k=0;k<=NZ-1;k++)
                  u0[i][k]=u1[i][k];
               for(i=0;i<=NX-1;i++)
               for(k=0;k<=NZ-1;k++)
                  u1[i][k]=u2[i][k];
           }

      }

    

}