ori-attn.rar_Psuedospectral_The Given_harmonic generation_solve_

function solve_sps2d(med,iters,excit,attn,record) % solve_sps2d(med,iters,excit,attn,record) % % MATLAB function to implement 2-D solution of % the nonlinear wave equation using a psuedospectral % / staggered A-B method, including attenuation modeled % with two relaxation mechanisms. % % This version accepts an arbitrary excitation (forcing) % function(s) at arbitrary points on the 2-D grid. % % Martin Anderson % Created: 1-26-00 Last modified: 11-1-00 % % Based on: % Wojcik G. L., B. Fornberg, R. Waag, L. Carcione, J. Mould, % L. Nikodym, and T. Driscoll, "Pseudospectral methods for % large-scale bioacoustic models," Proc. IEEE-UFFC % Ultrasonics Symposium, pp. 1501-1506, 1997. % % Nachman, A. I., J. F. Smith, and R. C. Waag, "An equation for % acoustic propagation in inhomogeneous media with relaxation % losses," JASA 88, pp. 1584-1595, 1990. % % Yuan, X., D. Borup, J. W. Wiskin, M. Berggren, R. Eidens, and % S. Johnson, "Formulation and validation of Berenger's PML % absorbing boundary for the FDTD simulation of acoustic % scattering," IEEE Trans. UFFC 44, 816-822, 1997. % % Yuan, X., D. Borup, J. W. Wiskin, M. Berggren, and % S. Johnson, "Simulation of acoustic wave propagation in % dispersive media with relaxation losses by using FDTD method % with PML absorbing boundary conditions," IEEE Trans. % UFFC 46, 14-23, 1999. % % G. L. Wojcik, J. Mould, S. Ayter and L. M. Carcione, "A study of % second harmonic generation by focused medical transducer pulses," % Proc. IEEE-UFFC Ultrasonics Symposium, pp. 1583-1588, 1998. % % Ghrist, M., B. Fornberg, and T. A. Driscoll, "Staggered time % integrators for wave equations," SIAM Journal on Numerical Analysis % 38(3), pp. 718-741, 2000. % % solve_sps2d.m accepts structures containing various variables % decribed in the accompanying READ_ME_solve_sps2d.txt file. % --------------------------------- nonlin_flag = (sum(med.ba(:)) > 0); % Calculate time increment [s] dt = med.stabf*min([med.dx med.dz])/max(med.c(:)); % Calculation of pressure and velocity are "staggered", % i.e. velocity at integer time steps [0 1 2...] * dt, % and pressure at half-integer time steps [0.5 1.5 2.5 ...] * dt; dt2 = dt/2; % Calculate bulk modulus matrix: K = med.rho.*med.c.^2.*ones(med.num_axpts,med.num_latpts); % Define axes for output image: xax = [-med.num_latpts/2+0.5:med.num_latpts/2-0.5]*med.dx; zax = [-med.num_pml:med.num_axpts-med.num_pml-1]*med.dz; % Define Perfectly Matched Boundary Layer (PML) % Only one taper is used for both x and y axes. % After Section IV of Yuan, et al., 1997: alpha_max = -log(med.pml_ampred).*med.c ./ (K.*med.dx); quadratic = ([med.num_pml:-1:1].'/med.num_pml).^2; % After equation 4 in Wojcik, et al., 1997: % Note the factor of -K actually cancels out. % (Included here for clarity w.r.t. references.) left_alpha = -K(:,1:med.num_pml) .*... alpha_max(:,1:med.num_pml) .*... repmat(quadratic.',med.num_axpts,1); right_alpha = -K(:,med.num_latpts-med.num_pml+1:med.num_latpts) .*... alpha_max(:,med.num_latpts-med.num_pml+1:med.num_latpts) .*... repmat(fliplr(quadratic.'),med.num_axpts,1); top_alpha = -K(1:med.num_pml,:) .*... alpha_max(1:med.num_pml,:) .*... repmat(quadratic,1,med.num_latpts); bottom_alpha = -K(med.num_axpts-med.num_pml+1:med.num_axpts,:) .*... alpha_max(med.num_axpts-med.num_pml+1:med.num_axpts,:).*... repmat(flipud(quadratic),1,med.num_latpts); % Define spatial frequency multipliers used to calculate % spectral derivatives: ps_axshifts = [[0:1:med.num_axpts/2-1] [-med.num_axpts/2:1:-1]].'*i*pi; ps_axshifts = repmat(ps_axshifts,1,med.num_latpts); ps_latshifts = [[0:1:med.num_latpts/2-1] [-med.num_latpts/2:1:-1]].'*i*pi; ps_latshifts = repmat(ps_latshifts,1,med.num_axpts); fft_axscale = 1/(med.num_axpts*med.dz/2); fft_latscale = 1/(med.num_latpts*med.dx/2); % Calculate Adams-Bashforth coefficients for ABS4 % integrator in Ghrist, et al. (2000): ab1 = (dt/24) * 26; ab2 = (dt/24) * -5; ab3 = (dt/24) * 4 ; ab4 = (dt/24) * -1; % --------------------------------- % Initialize pressure and velocity vectors: p = zeros(med.num_axpts,med.num_latpts); px = zeros(med.num_axpts,med.num_latpts); pz = zeros(med.num_axpts,med.num_latpts); vx = zeros(med.num_axpts,med.num_latpts); vz = zeros(med.num_axpts,med.num_latpts); dpx = zeros(med.num_axpts,med.num_latpts); dpz = zeros(med.num_axpts,med.num_latpts); dvx = zeros(med.num_axpts,med.num_latpts); dvz = zeros(med.num_axpts,med.num_latpts); div_u = zeros(med.num_axpts,med.num_latpts); % Initialize relaxation mechanism state variables: if attn.flag ~= 0 S = zeros(med.num_axpts,med.num_latpts,4,5); dS = zeros(med.num_axpts,med.num_latpts,4,5); kt1 = attn.kappas1./attn.taus1; kt2 = attn.kappas2./attn.taus2; tk_term = kt1 + kt2; end % Initialize W and dW matrices: % (These are 4-D work matrices that are % the size of the spatial domain in the % first two dimensions stacked as [vx vz px pz] % and d[vx vz px pz]/dt % across the third dimension for linear conditions, % [vx vz px pz div_ux div_uy] and % d[vx vz px pz div_ux div_uy]/dt for nonlinear, % with the 5 timesteps of these used in the % forward time integration stacked across the % fourth dimension.) if nonlin_flag ~= 0 w = zeros(med.num_axpts,med.num_latpts,6,5); dw = zeros(med.num_axpts,med.num_latpts,6,5); else w = zeros(med.num_axpts,med.num_latpts,4,5); dw = zeros(med.num_axpts,med.num_latpts,4,5); end % Initialize matrix for storing time ("hydrophone") record at focus: psf = zeros(iters,med.num_latpts); % Start clocks and initialize time counters: tic t0 = clock; start_time = cputime; forc_time = 0; fft_time = 0; relax_time = 0; pml_time = 0; adbash_time = 0; update_time = 0; output_time = 0; % Initialize time-step pointer vector: % (When the work matrices W and dW are updated % with each time step, their contents from % previous four time-steps are actually left % "in place", and a set of pointers to the % time-steps are updated. This saves computation % time and overhead.) tsp = [1:5]; % ------------------------------------ % Main calculation loop clf rec_count = 0; for cr = 1:iters % Integer time step: calculation of velocity % Calculate spatial partial derivatives: % Apply forcing function at points beyond the attenuation taper: t_stage = cputime; if cr<=size(excit.p,1)/2 % The next loop is not efficient, but allows for completely % arbitrary excitation locations throughout the medium: % (Usually very little time is spent on excitation, anyway.) for cpt = 1:length(excit.ax) p(excit.ax(cpt),excit.lat(cpt)) = excit.p(cr*2-1,cpt); end end forc_time = forc_time + (cputime-t_stage); t_stage = cputime; dpx = real(ifft(fft(p.').*ps_latshifts)).' *fft_latscale; dpz = real(ifft(fft(p).*ps_axshifts)) *fft_axscale; fft_time = fft_time + (cputime-t_stage); % --------------------------------- % Calculate relaxation mechanism state variables: % (Based on equations 17 and 18' in Yuan, et al., 1999) t_stage = cputime; if attn.flag ~= 0 sum1x = S(:,:,1,tsp(4))./attn.taus1 + S(:,:,2,tsp(4))./attn.taus2; sum2x = px*tk_term; dS(:,:,1,tsp(4)) = -S(:,:,1,tsp(4))./attn.taus1 + px.*kt1; dS(:,:,2,tsp(4)) = -S(:,:,2,tsp(4))./attn.taus2 + px.*kt2; sum1z = S(:,:,3,tsp(4))./attn.taus1 + S(:,:,4,tsp(4))./attn.taus2; sum2z = pz*tk_term; dS(:,:,3,tsp(4)) = -S(:,:,3,tsp(4))./attn.taus1 + pz.*kt1; dS(:,:,4,tsp(4)) = -S(:,:,4,tsp(4))./attn.taus2 + pz.*kt2; end relax_time = relax_time + (cputime-t_stage); % --------------------------------- % Calculate PML boundary region masks: t_stage = cputime; pml_vxleft = left_alpha .* vx(:,1:med.num_pml); pml_vxright = right_alpha .* ... vx(:,med.num_latpts-med.num_pml+1:med.num_latpts); pml_vztop = top_alpha .* vz(1:med.num_pml,:); pml_vzbottom = bottom_alpha .* ... vz(med.num_axpts-med.num_pml+1:med.num_axpts,:); % Calculate ODE co