sufdmod2_pml.c

su 的源代码库

   }


   /* Calculate w for left edge of bottom pad */

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


   /* Calculate w for main part of bottom pad */

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


   /* Calculate v along top and bottom edge of right pad */

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

      v_r [ix][   0] = cax_r [ix][   0] * v_r [ix][   0] +
                       cbx_r [ix][   0] * (ux_r [ix][ 0] + uz_r [ix][ 0]
                                          -ux_b [jx][jz] - uz_b [jx][jz]);

      v_r [ix][nz-1] = cax_r [ix][nz-1] * v_r [ix][nz-1] +
                       cbx_r [ix][nz-1] * (ux_b [jx][   0] + uz_b [jx][   0]
                                          -ux_r [ix][nz-2] - uz_r [ix][nz-2]);
   }


   /* Calculate v in rest of right pad */

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


   /* Calculate w along left and right sides of right pad */

   for (iz=0, jx=pml_thickness+2; iz<nz; ++iz) {
      w_r [ 0][iz] = caz_r [ 0][iz] * w_r [ 0][iz] +
                     cbz_r [ 0][iz] * (((abs[3]!=0) ? p [nx-2][iz] : 0.0)
                                      -ux_r [   0][iz] - uz_r [   0][iz]);

      w_r [jx][iz] = caz_r [jx][iz] * w_r [jx][iz] +
                     cbz_r [jx][iz] * (ux_r [jx-1][iz] + uz_r [jx-1][iz]
                                      -((abs[1]!=0) ? p [1][iz] : 0.0));
   }


   /* Calculate w in main part of right pad */

   for (ix=1; ix<pml_thickness+2; ++ix) {
      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, float **od, int verbose)
{
   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 (od!=NULL) {
	    /*	    if (verbose) warn("initializing right pml with density");*/
	    if (ix == 0) {
	       dvv_0 = sqrt (dvv [nx-1][iz] * od [nx-1][iz]);
	    } else if (ix == pml_thickness+1) {
	       dvv_0 = sqrt (dvv [   0][iz] * od [   0][iz]);
	    } else {
	       dvv_0 = sqrt (dvv [nx-1][iz] * od [nx-1][iz]);
	       dvv_1 = sqrt (dvv [   0][iz] * od [   0][iz]);

	       dvv_0 = ((ix) * dvv_1 + (1+pml_thickness-ix)*dvv_0);
	       dvv_0 /= (pml_thickness+1);
	    }
	 }
	 else {
	    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 (od!=NULL) {
	    /*	    if (verbose) warn("initializing bottom pml with density");*/
	    if (ix < nx) {
	       if (iz == 0) {
		  dvv_0 = sqrt (dvv [ix][nz-1] * od [ix][nz-1]);
	       } else if (iz == pml_thickness+1) {
		  dvv_0 = sqrt (dvv [ix][   0] * od [ix][   0]);
	       } else {
		  dvv_0 = sqrt (dvv [ix][nz-1] * od [ix][nz-1]);
		  dvv_1 = sqrt (dvv [ix][   0] * od [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] * od [nx-1][nz-1]);
		  dvv_1 = sqrt (dvv [   0][nz-1] * od [   0][nz-1]);
	       } else if (iz == pml_thickness+1) {
		  dvv_0 = sqrt (dvv [nx-1][   0] * od [nx-1][   0]);
		  dvv_1 = sqrt (dvv [   0][   0] * od [   0][   0]);
	       } else {
		  dvv_2 = sqrt (dvv [nx-1][nz-1] * od [nx-1][nz-1]);
		  dvv_3 = sqrt (dvv [nx-1][   0] * od [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] * od [0][nz-1]);
		  dvv_3 = sqrt (dvv [0][0] * od [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);
	    }
	 }

	 else {
	    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));
      }
   }
}