差分格式的二维波动方程C语言

img


/*
 * @description:
 * @Author: PingMin Cheung
 * @Date: 2023-06-05 11:32:28
 * @LastEditors: PingMin Cheung
 * @LastEditTime: 2023-06-05 11:56:40
 * @FilePath: /modeling/fdmod2dacoustic.v1.0.c
 *
 * Copyright (c) 2023 by WaveTomo, All Rights Reserved.
 */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
//
#include <omp.h>
//
#include <cwp.h>
#include <par.h>
#include <su.h>
//
#include <segy.h>
/*********************** self documentation **********************/
char *sdoc[] = {
    /**/
    "               test                                      ",
    "write the manual of this program here line by line       ", NULL};
/**
 * Ricker wavelet
 * @description: Ricker wavelet
 * @param {float} *src: Ricker wavelet
 * @param {float} A: peak amplitude
 * @param {float} fpeak: peak frequency
 * @param {float} delay: delay time
 * @param {float} dt: the time step
 * @param {int} nt: umber of samples along t
 * @LastEditors: LaurentCheung
 * @Version: 1.0
 */
void riker(float *src, float A, float fpeak, float delay, float dt, int nt);
/**
 * padding 2D model
 * @description: padding 2D model
 * @param {float**} model: model array
 * @param {int} NX: the sample number along 2nd dimension
 * @param {int} NZ: the sample number along 1st dimension
 * @param {int} nlayer: padding width
 * @LastEditors: LaurentCheung
 * @Version: 1.0
 */
void paddingModel2D(float **model, int NX, int NZ, int nlayer);
/**
 * 2D MEAL boundary condition
 * @description: 2D MEAL boundary condition
 * @param {float**} D: the value of damping function
 * @param {float**} Beta: scaling factor
 * @param {float} vpmax: the max element of vp array
 * @param {float} fpeak: the dominant frequency of the source defined above,
 * in order to study the influence of that parameter on the solution
 * @param {float} eal_alpha: the order of Gaussian weight function
 * @param {int} NZ: The sample number along 1st dimension
 * @param {int} NX: The sample number along 2nd dimension
 * @param {float} dz: the sampling interval along z, the 1st dimension
 * @param {float} dx: the sampling interval along x, the 2nd dimension
 * @param {int} nlayer: number of EAL
 * @return {*}
 * @LastEditors: LaurentCheung
 * @Version: 1.0
 */
void meal_coefs2d(float **D, float **Beta, float vpmax, float fpeak,
                  float eal_alpha, int NZ, int NX, float dz, float dx,
                  int nlayer);
int main(int argc, char **argv) {
  /* hook up getpar to handle the parameters */
  initargs(argc, argv);
  requestdoc(0);
  /* read parameters */
    int nz,nx,nt,pml_thick;
    float dx,dz,dt;
    char *vpfile=NULL;
    getparint("nz",&nz);
    getparint("nx",&nx);
    getparint("nt",&nt);
    getparint("pml_thick",&pml_thick);
    getparfloat("dz",&dz);
    getparfloat("dx",&dx);
    getparfloat("dt",&dt);
    getparstring("vpfile",&vpfile);
    warn("nz=%d,nx=%d,nt=%d.",nz,nx,nt);
    warn("dz=%f,dx=%f,dt=%f.",dz,dx,dt);
    warn("vpfile=%s.",vpfile);
  /* read model */
    float **vp=NULL;
    vp=alloc2float(nz,nx);
    FILE *fp;
    fp=fopen(vpfile,"r");
    fread(vp[0],FSIZE,nz*nx,fp);
    fclose(fp);

    fp=fopen("vp.bin","w");
    fwrite(vp[0].FSIZE,nz*nx,fp);
    fclose(fp);
  /* make source wavelet */
    float *src=NULL;
    src=alloc1float(nt);
    riker(src,1.0f,10.0f,0.15f,dt,nt);
    fp=fopen("src.bin","w");
    fwrite(src,FSIZE,nt,fp);
    fclose(fp);
  /* absorbing boundary condition */

  /* initial condition */
    float **p,**pn,**pn_1;
    p=alloc2float(nz,nx);
    pn=alloc2float(nz,nx);
    pn_1=alloc2float(nz,nx);
    int ix,iz;
    for(ix=0;ix<nx;ix++)
    {
        for(iz=0;iz<nz;iz++)
        {
            p[ix][iz]=0.0f;
            pn[ix][iz]=0.0f;
            pn_1[ix][iz]=0.0f;
        }
    }
  /* finite differential coefficient */

  /* finite differece scheme */
    int src_z,src_x;
    src_z=100;
    src_x=500;
    int it;
    float dpdxx,dpdzz;
    for(it=0;it<nt;it++)
    {
        for(ix=1;ix<nx-1;ix++)
        {
            for(iz=1;iz<nz-1;iz++)
            {
                dpdxx=(pn[ix+1][iz]-2.0f*pn[ix][iz]+pn[ix-1][iz])/dx/dx;
                dpdzz=(pn[ix][iz+1]-2.0f*pn[ix][iz]+pn[ix][iz-1])/dz/dz;
                p[ix][iz]=2.0f*pn[ix][iz]-pn_1[ix][iz]+dt*dt*vp[ix][iz]/vp[ix][iz]*(dpdxx+dpdzz);
            }
        }
        p[src_x][src_z]+=src[it];
        for(ix=0;ix<nx;ix++)
        {
            for(iz=0;iz<nz;iz++)
            {
                pn_1[ix][iz]=pn[ix][iz];
                pn[ix][iz]=p[ix][iz];
            }
        }
    }
  /* output snap */
    fp=fopen("snap.bin","w");
    fwrite(p[0],FSIZE,nz*nx,fp);
    fclose(fp);
  /* free memery */
    free2float(vp);
    free1float(src);
    free2float(p);
    free2float(pn);
    free2float(pn_1);
  return (CWP_Exit());
}
/**
 * Ricker wavelet
 * @description: Ricker wavelet
 * @param {float} *src: Ricker wavelet
 * @param {float} A: peak amplitude
 * @param {float} fpeak: peak frequency
 * @param {float} delay: delay time
 * @param {float} dt: the time step
 * @param {int} nt: umber of samples along t
 * @LastEditors: LaurentCheung
 * @Version: 1.0
 */
void riker(float *src, float A, float fpeak, float delay, float dt, int nt) {
  for (int i = 0; i < nt; i++) {
    if (i * dt <= 2.0f * delay) {
      src[i] =
          A *
          (1.0f - 2.0f * powf(3.1415926f * fpeak * (i * dt - delay), .0f2)) *
          expf(-1.0f * powf(3.1415926f * fpeak * (i * dt - delay), 2.0f));
    } else {
      src[i] = 0.0f;
    }
  }
}
/**
 * padding 2D model
 * @description: padding 2D model
 * @param {float**} model: model array
 * @param {int} NX: the sample number along 2nd dimension
 * @param {int} NZ: the sample number along 1st dimension
 * @param {int} nlayer: padding width
 * @LastEditors: LaurentCheung
 * @Version: 1.0
 */
void paddingModel2D(float **model, int NX, int NZ, int nlayer) {
  int i, j;
  /*padding top and bottom of modelelocity*/
  for (i = nlayer; i < NX - nlayer; i++) {
    for (j = 0; j < nlayer; j++) {
      model[i][j] = model[i][nlayer];
    }
    for (j = NZ - nlayer; j < NZ; j++) {
      model[i][j] = model[i][NZ - nlayer - 1];
    }
  }
  /*padding left and right of modelelocity*/
  for (i = 0; i < nlayer; i++) {
    for (j = 0; j < NZ; j++) {
      model[i][j] = model[nlayer][j];
      model[NX - i - 1][j] = model[NX - nlayer - 1][j];
    }
  }
}
float absorbWeight(float alpha, int l, int nlayer) {
  float w, d;
  d = (0.1f * l) / (0.1f * nlayer);
  w = expf(logf(2.0f) * powf(d, alpha)) - 1.0f;
  return w;
}
/**
 * 2D MEAL boundary condition
 * @description: 2D MEAL boundary condition
 * @param {float**} D: the value of damping function
 * @param {float**} Beta: scaling factor
 * @param {float} vpmax: the max element of vp array
 * @param {float} fpeak: the dominant frequency of the source defined above,
 * in order to study the influence of that parameter on the solution
 * @param {float} eal_alpha: the order of Gaussian weight function
 * @param {int} NZ: The sample number along 1st dimension
 * @param {int} NX: The sample number along 2nd dimension
 * @param {float} dz: the sampling interval along z, the 1st dimension
 * @param {float} dx: the sampling interval along x, the 2nd dimension
 * @param {int} nlayer: number of EAL
 * @return {*}
 * @LastEditors: LaurentCheung
 * @Version: 1.0
 */
void meal_coefs2d(float **D, float **Beta, float vpmax, float fpeak,
                  float eal_alpha, int NZ, int NX, float dz, float dx,
                  int nlayer) {
  float *wx, *wz, *bx, *bz;
  float R, d0_x, d0_z, beta0_x, beta0_z;
  int iz, ix, l;
  wx = alloc1float(NX);
  wz = alloc1float(NZ);
  bx = alloc1float(NX);
  bz = alloc1float(NZ);
  R = powf(10.0f, -1.0f * ((log10f(nlayer) - 1.0f) / (log10f(2.0f))) - 3.0f);
  d0_x = (3.0f * vpmax * logf(1.0f / R)) / (2.0f * nlayer * dx);
  d0_z = (3.0f * vpmax * logf(1.0f / R)) / (2.0f * nlayer * dz);
  beta0_x = (vpmax / sqrtf(3.0f)) / (5.0f * dx * fpeak);
  beta0_z = (vpmax / sqrtf(3.0f)) / (5.0f * dz * fpeak);
  if (eal_alpha == -1.0f) {
    eal_alpha = 0.05f * nlayer + 1.5f;
  }
  float weight;
  for (ix = 0; ix < NX; ix++) {
    bx[ix] = 1.0f;
    wx[ix] = 0.0f;
    // 鍖哄煙2锛堝乏杈圭晫銆佸彸杈圭晫锛?
    if (ix <= nlayer) {
      l = nlayer - ix;
      weight = absorbWeight(eal_alpha, l, nlayer);
      bx[ix] = 1.0f + (beta0_x - 1.0f) * weight;
      wx[ix] = d0_x * weight;
    } else if (ix >= NX - nlayer - 1) {
      l = ix - (NX - nlayer - 1);
      weight = absorbWeight(eal_alpha, l, nlayer);
      bx[ix] = 1.0f + (beta0_x - 1.0f) * weight;
      wx[ix] = d0_x * weight;
    }
  }
  for (iz = 0; iz < NZ; iz++) {
    bz[iz] = 1.0f;
    wz[iz] = 0.0f;
    // 鍖哄煙1锛堜笂杈圭晫锛屼笅杈圭晫锛?
    if (iz <= nlayer) {
      l = nlayer - iz;
      weight = absorbWeight(eal_alpha, l, nlayer);
      bz[iz] = 1.0f + (beta0_z - 1.0f) * weight;
      wz[iz] = d0_z * weight;
    } else if (iz >= NZ - nlayer - 1) {
      l = iz - (NZ - nlayer - 1);
      weight = absorbWeight(eal_alpha, l, nlayer);
      bz[iz] = 1.0f + (beta0_z - 1.0f) * weight;
      wz[iz] = d0_z * weight;
    }
  }
  //
#pragma omp parallel for private(ix, iz)
  for (ix = 0; ix < NX; ix++) {
#pragma ivdep
    for (iz = 0; iz < NZ; iz++) {
      D[ix][iz] = sqrtf(powf(wz[iz], 2.0f) + powf(wx[ix], 2.0f));
      Beta[ix][iz] =
          sqrtf(powf(bz[iz] - 1.0f, 2.0f) + powf(bx[ix] - 1.0f, 2.0f)) + 1.0f;
    }
  }
  free1float(wx);
  free1float(wz);
  free1float(bx);
  free1float(bz);
}

上面的代码是利用差分方法得到的二维声波方程，可不可以按照图片中的题目在这个的基础上把它改成差分格式的二维波动方程

  • 你可以参考下这个问题的回答, 看看是否对你有帮助, 链接: https://ask.csdn.net/questions/7516585
  • 这篇博客你也可以参考下：如何用C语言获取一个数二进制序列中所有的偶数位和奇数位，分别输出二进制序列
  • 你还可以看下c语言参考手册中的 c语言-内存模型与数据竞争
  • 除此之外, 这篇博客: 【C语言进阶】③探究浮点数在内存中的存储方式中的 二、浮点数在内存中存储方式 部分也许能够解决你的问题, 你可以仔细阅读以下内容或跳转源博客中阅读:
    1. IEEE二进制浮点数算术标准（IEEE 754）：是20世纪80年代以来最广泛使用的浮点数运算标准，为许多CPU与浮点运算器所采用。
    2. IEEE 754规定了四种表示浮点数值的方式：单精确度（32位）、双精确度（64位）、延伸单精确度（43比特以上，很少使用）与延伸双精确度（79比特以上，通常以80位实现）。
    • 以上内容来自于百度百科。
    • 简单点来说：一个任意的二进制浮点数A可以表示成下面的形式：
    1. （-1）^S * M * 2^E
      2.s=（-1)^S，s表示的是符号位，当s=0时，A为正数；当s=1时，A为负数；
    2. M表示有效数字，1 <= M < 2；
    3. 2^E表示指数位；

    举个例子：

    • 十进制的11.0，写成二进制为：1011.0，相当于1.011x2^3;
    • 那么按照上面的形式：s=0，M=1.011，E=3；
    • 而十进制的-11.0，写出二进制为：-1011.0，相当于-1.011x2^3;此时，s=1，M=1.011，E=3;
    • IEEE 754规定：
      对于32位的单精度浮点数而言，最高的1位是符号位s，接着的8位是指数E，剩下的23位为有效数字M。
    • 即：
      在这里插入图片描述

    实际上:

    • 前面说过，1<= M < 2，也就是说，M可以写成：1.xxxx的形式，其中xxxx表示小数部分。
      IEEE754规定：计算机内部保存M时，默认这个数的第一位是1，因此可以被舍去，只保留后面的小数部分；比如对二进制1.011来说，只保存.011，等到读取是，再在小数点前加1；这样做的目的是节省了1位有效数字。以单精度浮点数（32bit）来说，留给M的只有23位，将第一位1舍去后，等于可以保存24位有效数字；

    • 对于指数E：
      （1）E作为一个无符号整数，这意味着，如果E为8位（单精度浮点数），它的取值范围为：0-255；如果E为11位（双精度浮点数），它的取值为：0-2047；
      （2）但是，科学记数法中，指数位是可以为负数的，所以**IEEE 754规定:E存入内存时，必须在加上一个中间数，对于单精度浮点数来说，这个数是127，对于双精度浮点数来说，这个数是1023；
      例如：2^10，E为10，保存在内存中时为：10+127=137,即写成二进制为：1000 1001；

    • 然后指数E从内存中取出还可以分为三种情况：

      1.E不全为0或者全为1：

      这时，浮点数就采用下面规则表示，即指数E的计算值减去127（或者1023）得到真实值，在将真实值得有效数字M前加上第一位1，即得到二进制数；
      例如：
      十进制：0.5(1/2) 的二进制形式为0.1，由于规定正数部分必须为1，即将小数点右移动1位，为：1.0*2^(-1)，其E为：-1+127=126，表示为：0111 1110，而M为：1.0去掉整数部分1，为0，补齐到23位为：00000000000000000000000；
      则这个数的二进制在内存中保存形式为：
      0 0111110 0000000000000000000000

    2. E全为0
    这时，浮点数的指数E等于1-127（或者1-1023）即为真实值，
    有效数字M不再加上第一位的1，而是还原为0.xxxxxx的小数。这样做是为了表示±0，以及接近于
    0的很小的数字。

    3. E全为1
    这时，如果有效数字M全为0，表示±无穷大（正负取决于符号位s）;

  • 您还可以看一下 李飞老师的C语言开发之数据结构与算法二课程中的 树的插入、查找操作小节, 巩固相关知识点