fdmod.c

地震波正演和显示模块

      for (iz=0; iz<nz; ++iz) {
         w_r [ix][iz] = caz_r [ix][iz] * w_r [ix][iz] +
                        cbz_r [ix][iz] * (ux_r [ix-1][iz] + uz_r [ix-1][iz]
                                         -ux_r [ix  ][iz] - uz_r [ix  ][iz]);
      }
   }


   /* Calculate ux and uz in bottom pad */

   for (ix=0; ix<nx+pml_thickness-1; ++ix) {
      for (iz=0; iz<pml_thickness+2; ++iz) {
         ux_b [ix][iz] = dax_b [ix][iz] * ux_b [ix  ][iz] +
                         dbx_b [ix][iz] * (w_b [ix][iz  ] - w_b [ix+1][iz]);
      }
   }

   for (ix=nx+pml_thickness-1, iz=0; iz<pml_thickness+2; ++iz) {
      ux_b [ix][iz] = dax_b [ix][iz] * ux_b [ix  ][iz] + 
                      dbx_b [ix][iz] * (w_b [ix][iz  ] - w_b [   0][iz]);
   }


   for (ix=0; ix<nx+pml_thickness; ++ix) {
      for (iz=0; iz<pml_thickness+2; ++iz) {
         uz_b [ix][iz] = daz_b [ix][iz] * uz_b [ix][iz  ] +
                         dbz_b [ix][iz] * (v_b [ix][iz+1] - v_b [ix][iz  ]);
      }
   }


   /* Calculate ux and uz in right pad */

   for (ix=0; ix<pml_thickness+2; ++ix) {
      for (iz=0; iz<nz; ++iz) {
         ux_r [ix][iz] = dax_r [ix][iz] * ux_r [ix  ][iz] +
                         dbx_r [ix][iz] * (w_r [ix  ][iz] - w_r [ix+1][iz]);
      }
   }

   for (ix=0; ix<pml_thickness+2; ++ix) {
      for (iz=0; iz<nz-1; ++iz) {
         uz_r [ix][iz] = daz_r [ix][iz] * uz_r [ix][iz  ] +
                         dbz_r [ix][iz] * (v_r [ix][iz+1] - v_r [ix][iz  ]);
      }
   }

   for (ix=0, iz=nz-1, jx=nx-1; ix<pml_thickness+2; ++ix, ++jx) {
      if (jx == nx+pml_thickness) jx = 0;

      uz_r [ix][ 0] = uz_b [jx][pml_thickness+1];
      uz_r [ix][iz] = uz_b [jx][              0];
   }


   /* Update top and bottom edge of main grid with new field values */

   for (ix=0, jz=pml_thickness+1; ix<nx; ++ix) {
      if (abs [0] != 0) pp [ix][   0] = ux_b [ix][jz] + uz_b [ix][jz];
      if (abs [2] != 0) pp [ix][nz-1] = ux_b [ix][ 0] + uz_b [ix][ 0];
   }


   /* Update left and right edges of main grid with new field values */

   for (iz=1, jx=pml_thickness+1; iz<nz-1; ++iz) {
      if (abs [1] != 0) pp [   0][iz] = ux_r [jx][iz] + uz_r [jx][iz];
      if (abs [3] != 0) pp [nx-1][iz] = ux_r [ 0][iz] + uz_r [ 0][iz];
   }
}


static void pml_init (int nx, int nz, float dx, float dz, float dt, float **dvv)
{
   int ix, iz;

   /* Allocate arrays for pad on right */

   cax_r = alloc2float (nz, pml_thickness+2);
   cbx_r = alloc2float (nz, pml_thickness+2);
   caz_r = alloc2float (nz, pml_thickness+3);
   cbz_r = alloc2float (nz, pml_thickness+3);
   dax_r = alloc2float (nz, pml_thickness+2);
   dbx_r = alloc2float (nz, pml_thickness+2);
   daz_r = alloc2float (nz, pml_thickness+2);
   dbz_r = alloc2float (nz, pml_thickness+2);

   ux_r  = alloc2float (nz, pml_thickness+2);
   uz_r  = alloc2float (nz, pml_thickness+2);
   v_r   = alloc2float (nz, pml_thickness+2);
   w_r   = alloc2float (nz, pml_thickness+3);


   /* Zero out arrays for pad on right */

   for (ix=0; ix<pml_thickness+2; ++ix) {
      for (iz=0; iz<nz; ++iz) {
         ux_r  [ix][iz] = uz_r  [ix][iz] = 0.0;
         v_r   [ix][iz] = w_r   [ix][iz] = 0.0;

         cax_r [ix][iz] = cbx_r [ix][iz] = 0.0;
         caz_r [ix][iz] = cbz_r [ix][iz] = 0.0;
         dax_r [ix][iz] = dbx_r [ix][iz] = 0.0;
         daz_r [ix][iz] = dbz_r [ix][iz] = 0.0;
      }
   }


   /* Zero out extra bit on right pad */

   for (ix=pml_thickness+2, iz=0; iz<nz; ++iz) {
      caz_r [ix][iz] = cbz_r [ix][iz] = 0.0;
      w_r   [ix][iz] = 0.0;
   }


   /* Allocate arrays for pad on bottom */

   cax_b = alloc2float (pml_thickness+3, nx + pml_thickness);
   cbx_b = alloc2float (pml_thickness+3, nx + pml_thickness);
   caz_b = alloc2float (pml_thickness+2, nx + pml_thickness);
   cbz_b = alloc2float (pml_thickness+2, nx + pml_thickness);
   dax_b = alloc2float (pml_thickness+2, nx + pml_thickness);
   dbx_b = alloc2float (pml_thickness+2, nx + pml_thickness);
   daz_b = alloc2float (pml_thickness+2, nx + pml_thickness);
   dbz_b = alloc2float (pml_thickness+2, nx + pml_thickness);

   ux_b  = alloc2float (pml_thickness+2, nx + pml_thickness);
   uz_b  = alloc2float (pml_thickness+2, nx + pml_thickness);
   v_b   = alloc2float (pml_thickness+3, nx + pml_thickness);
   w_b   = alloc2float (pml_thickness+2, nx + pml_thickness);
	

   /* Zero out arrays for pad on bottom */

   for (ix=0; ix<nx+pml_thickness; ++ix) {
      for (iz=0; iz<pml_thickness+2; ++iz) {
         ux_b  [ix][iz] = uz_b  [ix][iz] = 0.0;
         v_b   [ix][iz] = w_b   [ix][iz] = 0.0;

         cax_b [ix][iz] = cbx_b [ix][iz] = 0.0;
         caz_b [ix][iz] = cbz_b [ix][iz] = 0.0;
         dax_b [ix][iz] = dbx_b [ix][iz] = 0.0;
         daz_b [ix][iz] = dbz_b [ix][iz] = 0.0;
      }
   }


   /* Zero out extra bit on bottom pad */

   for (ix=0, iz=pml_thickness+2; ix<nx+pml_thickness; ++ix) {
      cax_b [ix][iz] = cbx_b [ix][iz] = 0.0;
      v_b   [ix][iz] = 0.0;
   }


   /* Initialize cax & cbx arrays */

   for (ix=0; ix<pml_thickness+2; ++ix) {
      for (iz=0; iz<nz; ++iz) {
         sigma_ez = 0.0;

         cax_r [ix][iz] = (2.0 - (sigma_ez * dt)) / (2.0 + (sigma_ez * dt));
         cbx_r [ix][iz] = (2.0 * dt / dz)         / (2.0 + (sigma_ez * dt));
      }
   }

   for (ix=0; ix<nx+pml_thickness; ++ix) {
      for (iz=0; iz<pml_thickness+3; ++iz) {
         if ((iz == 0) || (iz == pml_thickness + 2)) {
            sigma_ez = 0.0;
         } else {
            sigma_ez = pml_max * 0.5 * (1.0 - cos (2*PI*(iz-0.5)/(pml_thickness+1)));
         }

         cax_b [ix][iz] = (2.0 - (sigma_ez * dt)) / (2.0 + (sigma_ez * dt));
         cbx_b [ix][iz] = (2.0 * dt / dz)         / (2.0 + (sigma_ez * dt));
      }
   }


   /* Initialize caz & cbz arrays */

   for (ix=0; ix<pml_thickness+3; ++ix) {
      for (iz=0; iz<nz; ++iz) {
         if ((ix == 0) || (ix == pml_thickness+2)) {
            sigma_ex = 0.0;
         } else {
            sigma_ex = pml_max * 0.5 * (1.0 - cos (2*PI*(ix-0.5)/(pml_thickness+1)));
         }

         caz_r [ix][iz] = (2.0 - (sigma_ex * dt)) / (2.0 + (sigma_ex * dt));
         cbz_r [ix][iz] = (2.0 * dt / dz)         / (2.0 + (sigma_ex * dt));
      }
   }


   for (ix=0; ix<nx+pml_thickness; ++ix) {
      for (iz=0; iz<pml_thickness+2; ++iz) {
         if (ix == 0) {
            sigma_ex = pml_max * 0.5 * (1.0 - cos (2*PI*(0.5)/(pml_thickness+1)));
         } else if (ix < nx) {
            sigma_ex = 0.0;
         } else {
            sigma_ex = pml_max * 0.5 * (1.0 - cos (2*PI*(ix-nx+0.5)/(pml_thickness+1)));
         }

         caz_b [ix][iz] = (2.0 - (sigma_ex * dt)) / (2.0 + (sigma_ex * dt));
         cbz_b [ix][iz] = (2.0 * dt / dz)         / (2.0 + (sigma_ex * dt));
      }
   }


   /* Initialize right pad's dax dbx, daz, & dbz arrays */

   for (ix=0; ix<pml_thickness+2; ++ix) {
      for (iz=0; iz<nz; ++iz) {

         /* Determine sigma_mx and sigma_mz */

         if ((ix == 0) || (ix == pml_thickness+1)) {
            sigma_mx = 0.0;
         } else {
            sigma_mx = pml_max * 0.5 * (1.0 - cos (2*PI*(ix)/(pml_thickness+1)));
         }

         sigma_mz = 0.0;


         /* Determine velocity, interpolate */

         if (ix == 0) {
            dvv_0 = sqrt (dvv [nx-1][iz]);
         } else if (ix == pml_thickness+1) {
            dvv_0 = sqrt (dvv [   0][iz]);
         } else {
            dvv_0 = sqrt (dvv [nx-1][iz]);
            dvv_1 = sqrt (dvv [   0][iz]);

            dvv_0 = ((ix) * dvv_1 + (1+pml_thickness-ix)*dvv_0);
            dvv_0 /= (pml_thickness+1);
         }

         dvv_0 = dvv_0 * dvv_0;


         dax_r [ix][iz] = (2.0 - (sigma_mx * dt)) / (2.0 + (sigma_mx * dt));
         dbx_r [ix][iz] = (2.0 * dt * dvv_0 / dx) / (2.0 + (sigma_mx * dt));

         daz_r [ix][iz] = (2.0 - (sigma_mz * dt)) / (2.0 + (sigma_mz * dt));
         dbz_r [ix][iz] = (2.0 * dt * dvv_0 / dz) / (2.0 + (sigma_mz * dt));
      }
   }


   /* Initialize bottom pad's dax, dbx, daz, & dbz arrays */

   for (ix=0; ix<nx+pml_thickness; ++ix) {
      for (iz=0; iz<pml_thickness+2; ++iz) {

         /* Determine sigma_mx and sigma_mz */

         if (ix < nx) {
            sigma_mx = 0.0;
         } else {
            sigma_mx = pml_max * 0.5 * (1.0 - cos (2*PI*(ix-nx+1)/(pml_thickness+1)));
         }

         if ((iz == 0) || (iz == pml_thickness+1)) {
            sigma_mz = 0.0;
         } else {
            sigma_mz = pml_max * 0.5 * (1.0 - cos (2*PI*(iz)/(pml_thickness+1)));
         }


         /* Determine velocity, interpolate */

         if (ix < nx) {
            if (iz == 0) {
               dvv_0 = sqrt (dvv [ix][nz-1]);
            } else if (iz == pml_thickness+1) {
               dvv_0 = sqrt (dvv [ix][   0]);
            } else {
               dvv_0 = sqrt (dvv [ix][nz-1]);
               dvv_1 = sqrt (dvv [ix][   0]);

               dvv_0 = ((iz) * dvv_1 + (1+pml_thickness-iz)*dvv_0);
               dvv_0 /= (pml_thickness+1);
            }
         } else {
            if (iz == 0) {
               dvv_0 = sqrt (dvv [nx-1][nz-1]);
               dvv_1 = sqrt (dvv [   0][nz-1]);
            } else if (iz == pml_thickness+1) {
               dvv_0 = sqrt (dvv [nx-1][   0]);
               dvv_1 = sqrt (dvv [   0][   0]);
            } else {
               dvv_2 = sqrt (dvv [nx-1][nz-1]);
               dvv_3 = sqrt (dvv [nx-1][   0]);

               dvv_0 = ((iz) * dvv_3 + (1+pml_thickness-iz)*dvv_2);
               dvv_0 /= (pml_thickness+1);

               dvv_2 = sqrt (dvv [0][nz-1]);
               dvv_3 = sqrt (dvv [0][   0]);

               dvv_1 = ((iz) * dvv_3 + (1+pml_thickness-iz)*dvv_2);
               dvv_1 /= (pml_thickness+1);
            }

            dvv_0 = ((ix-nx+1) * dvv_1 + (nx+pml_thickness-ix)*dvv_0);
            dvv_0 /= (pml_thickness+1);
         }

         dvv_0 = dvv_0 * dvv_0;

         dax_b [ix][iz] = (2.0 - (sigma_mx * dt)) / (2.0 + (sigma_mx * dt));
         dbx_b [ix][iz] = (2.0 * dt * dvv_0 / dx) / (2.0 + (sigma_mx * dt));

         daz_b [ix][iz] = (2.0 - (sigma_mz * dt)) / (2.0 + (sigma_mz * dt));
         dbz_b [ix][iz] = (2.0 * dt * dvv_0 / dz) / (2.0 + (sigma_mz * dt));
      }
   }
}

static void pml_free(){

   free2float( cax_r);
   free2float(cbx_r  );
  free2float( caz_r );
  free2float( cbz_r );
  free2float( dax_r );
  free2float( dbx_r );
  free2float( daz_r );
  free2float( dbz_r );

  free2float(  ux_r  );
   free2float(uz_r  );
   free2float(v_r   );
  free2float( w_r   );
  free2float( cax_b );
  free2float( cbx_b );
  free2float( caz_b );
  free2float( cbz_b );
  free2float( dax_b );
  free2float( dbx_b );
  free2float( daz_b );
  free2float( dbz_b );

  free2float( ux_b );
  free2float( uz_b );
  free2float( v_b   );
  free2float( w_b   );
}

void intcub (int ideriv, int nin, float xin[], float ydin[][4],
        int nout, float xout[], float yout[])
/*****************************************************************************
evaluate y(x), y'(x), y''(x), ... via piecewise cubic interpolation
******************************************************************************
Input:
ideriv          =0 if y(x) desired; =1 if y'(x) desired, ...
nin             length of xin and ydin arrays
xin             array[nin] of monotonically increasing or decreasing x values
ydin            array[nin][4] of y(x), y'(x), y''(x), and y'''(x)
nout            length of xout and yout arrays
xout            array[nout] of x values at which to evaluate y(x), y'(x), ...

Output:
yout            array[nout] of y(x), y'(x), ... values
******************************************************************************
Notes:
xin values must be monotonically increasing or decreasing.

Extrapolation of the function y(x) for xout values outside the range
spanned by the xin values is performed using the derivatives in
ydin[0][0:3] or ydin[nin-1][0:3], depending on whether xout is closest
to xin[0] or xin[nin-1], respectively.
******************************************************************************
Author:  Dave Hale, Colorado School of Mines, 06/02/89
*****************************************************************************/
{
        static int idx;
        int iout;
        float delx;

        /* y(x) is desired, then */
        if (ideriv==0) {
                for (iout=0; iout<nout; iout++) {
                        xindex(nin,xin,xout[iout],&idx);
                        delx = xout[iout]-xin[idx];
                        yout[iout] = (ydin[idx][0]+delx*
                                (ydin[idx][1]+delx*
                                (ydin[idx][2]*O2+delx*
                                (ydin[idx][3]*O6))));
                }

        /* else, if y'(x) is desired, then */
        } else if (ideriv==1) {
                for (iout=0; iout<nout; iout++) {
                        xindex(nin,xin,xout[iout],&idx);
                        delx = xout[iout]-xin[idx];
                        yout[iout] = (ydin[idx][1]+delx*
                                (ydin[idx][2]+delx*
                                (ydin[idx][3]*O2)));
                }
  /* else, if y''(x) is desired, then */
        } else if (ideriv==2) {
                for (iout=0; iout<nout; iout++) {
                        xindex(nin,xin,xout[iout],&idx);
                        delx = xout[iout]-xin[idx];
                        yout[iout] = (ydin[idx][2]+delx*
                                (ydin[idx][3]));
                }

        /* else, if y'''(x) is desired, then */
        } else if (ideriv==3) {
                for (iout=0; iout<nout; iout++) {
                        xindex(nin,xin,xout[iout],&idx);