wavefield.zip_PML_PML边界_wavefield_模拟二维波场_边界波场

function p = defwave2d_pml(src,vel,at,ax,az,... ixsrc,izsrc,nxext,nzext,vm) ns = length(src); nx = length(ax); nz = length(az); nt = length(at); dx = ax(2) - ax(1); dz = az(2) - az(1); dt = at(2) - at(1); % For PML in x and z ---------------------------- [sigpmlxa,sigpmlxb] = defparamcpml1d(ax,vm,nxext); [sigpmlza,sigpmlzb] = defparamcpml1d(az,vm,nzext); % Extend the model ------------------------------ nxe = nx + 2*nxext; nze = nz + 2*nzext; % corners vele = zeros(nze,nxe); vele(1:nzext,1:nxext) = vel(1,1); vele(1:nzext,nxe-nxext+1:nxe) = vel(1,nx); vele(nze-nzext+1:nze,1:nxext) = vel(nz,1); vele(nze-nzext+1:nze,nxe-nxext+1:nxe) = vel(nz,nx); % edges for iz = 1:nzext vele(iz,1+nxext:nxe-nxext) = vel(1,:); vele(nze-nzext+iz,1+nxext:nxe-nxext) = vel(nz,:); end for ix = 1:nxext vele(1+nzext:nze-nzext,ix) = vel(:,1); vele(1+nzext:nze-nzext,nxe-nxext+ix) = vel(:,nx); end % inside vele(nzext+1:nze-nzext,nxext+1:nxe-nxext) = vel; %% Propagate p = zeros(nz,nx,nt); px = zeros(nze,nxe); pz = zeros(nze,nxe); ux = zeros(nze,nxe-1); uz = zeros(nze-1,nxe); pxn = zeros(nze,nxe); pzn = zeros(nze,nxe); uxn = zeros(nze,nxe-1); uzn = zeros(nze-1,nxe); for it = 1:nt-1 % % Derivative in space % dpdx = zeros(nze,nxe-1); dpdz = zeros(nze-1,nxe); % for ix = 1:nxe-1 dpdx(:,ix) = (px(:,ix+1)+pz(:,ix+1)-px(:,ix)-pz(:,ix))/dx; end for iz = 1:nze-1 dpdz(iz,:) = (px(iz+1,:)+pz(iz+1,:)-px(iz,:)-pz(iz,:))/dz; end % % Integration in time % for ix = 1:nxe-1 uxn(:,ix) = ux(:,ix) + dpdx(:,ix)*dt - sigpmlxb(ix) * ux(:,ix) * dt; end for iz = 1:nze-1 uzn(iz,:) = uz(iz,:) + dpdz(iz,:)*dt - sigpmlzb(iz) * uz(iz,:) * dt; end % % Derivative in space % duxdx = zeros(nze,nxe-1); duzdz = zeros(nze-1,nxe); % for ix = 2:nxe-1 duxdx(:,ix) = (uxn(:,ix) - uxn(:,ix-1))/dx; end for iz = 2:nze-1 duzdz(iz,:) = (uzn(iz,:) - uzn(iz-1,:))/dz; end % % Integration in time % for ix = 2:nxe-1 for iz = 2:nze-1 % if(ix == ixsrc + nxext) && (iz == izsrc + nzext) dirac = 1; else dirac = 0; end if(it <= ns) addsrc = src(it); else addsrc = 0; end % pxn(iz,ix) = px(iz,ix) + (duxdx(iz,ix) + dirac*addsrc)*dt*vele(iz,ix)^2 - ... sigpmlxa(ix)*px(iz,ix)*dt; pzn(iz,ix) = pz(iz,ix) + (duzdz(iz,ix) + dirac*addsrc)*dt*vele(iz,ix)^2 - ... sigpmlza(iz)*pz(iz,ix)*dt; end end % % Save the wavefield % p(:,:,it+1) = pxn(1+nzext:nze-nzext,1+nxext:nxe-nxext) + ... pzn(1+nzext:nze-nzext,1+nxext:nxe-nxext); % % Update % px = pxn; pz = pzn; ux = uxn; uz = uzn; % end