Code covered by the BSD License  

Highlights from
Demo files for "Parallel Computing with MATLAB on Multicore Desktops and GPUs" Webinar

image thumbnail

Demo files for "Parallel Computing with MATLAB on Multicore Desktops and GPUs" Webinar

by

 

05 May 2011 (Updated )

Demo files for "Parallel Computing with MATLAB on Multicore Desktops and GPUs" Webinar (May 4, 2011)

WaveEqn_CPU(h,h2,hObject,hText,N,maxIter)
function vv = WaveEqn_CPU(h,h2,hObject,hText,N,maxIter)
%% Solving 2nd Order Wave Equation Using Spectral Methods
% This example solves a 2nd order wave equation: utt = uxx + uyy, with u =
% 0 on the boundaries. It uses a 2nd order central finite difference in
% time and a Chebyshev spectral method in space (using FFT). The code has
% been modified from an example in Spectral Methods in MATLAB by Trefethen,
% Lloyd N.

% Copyright 2011 The MathWorks, Inc.

if nargin <= 2
   guiMode = false;
   N = h;
   if nargin == 1
      maxIter = 4000;
   else
      maxIter = h2;
   end
else
   guiMode = true;
   
   on_state = get(hObject,'Max');
end

x = cos(pi*(0:N)/N); % using Chebyshev points
y = x';
dt = 6/N^2;

[xx,yy] = meshgrid(x,y);

vv = exp(-40*((xx-.4).^2 + yy.^2));
vvold = vv;
ii = 2:N;

if guiMode
   plotstep = 10;
   
   [xxx,yyy] = meshgrid(-1:1/16:1,-1:1/16:1); % grid for plotting
   ylim = 5e4*10^-(log2(N));
   
   sh  = surf(h , xxx, yyy, zeros(size(xxx)));
   sh2 = surf(h2, xxx, yyy, zeros(size(xxx)));
   axis(h , [-1 1 -1 1 -0.1 1]);
   axis(h2, [-1 1 -1 1 -ylim ylim]);
   
   plotWave(sh,sh2,xx,yy,vv,vvold,xxx,yyy)
   set(hText, 'String', '0');
end

% Weights used for spectral differentiation via FFT
index1 = 1i*[0:N-1 0 1-N:-1];
W1T = repmat(index1,length(ii),1);
index2 = -[0:N 1-N:-1].^2;
W2T = repmat(index2,length(ii),1);
index3 = 1i*[0:N-1 0 1-N:-1]';
W3T = repmat(index3,1,length(ii));
index4 = -[0:N 1-N:-1]'.^2;
W4T = repmat(index4,1,length(ii));

WuxxT1 = repmat((1./(1-x(ii).^2)),length(ii),1);
WuxxT2 = repmat(x(ii)./(1-x(ii).^2).^(3/2),length(ii),1);
WuyyT1 = repmat(1./(1-y(ii).^2),1,length(ii));
WuyyT2 = repmat(y(ii)./(1-y(ii).^2).^(3/2),1,length(ii));

uxx = zeros(N+1,N+1);
uyy = zeros(N+1,N+1);

% Start time-stepping
for n = 1:maxIter
   
   if guiMode && ~isequal(get(hObject, 'Value'), on_state)
      break;
   end
   
   V = [vv(ii,:) vv(ii,N:-1:2)];
   U = real(fft(V.')).';
   
   W1test = (U.*W1T).';
   W2test = (U.*W2T).';
   W1 = (real(ifft(W1test))).';
   W2 = (real(ifft(W2test))).';
   
   uxx(ii,ii) = W2(:,ii).* WuxxT1 - W1(:,ii).*WuxxT2;
   
   V = [vv(:,ii); vv((N:-1:2),ii)];
   U = real(fft(V));
   
   W1 = real(ifft(U.*W3T));
   W2 = real(ifft(U.*W4T));
   
   uyy(ii,ii) = W2(ii,:).* WuyyT1 - W1(ii,:).*WuyyT2;
   vvnew = 2*vv - vvold + dt*dt*(uxx+uyy);
   vvold = vv; vv = vvnew;
   
   % plot every plotstep iterations
   if guiMode && mod(n,plotstep) == 0
      
      
      plotWave(sh,sh2,xx,yy,vv,vvold,xxx,yyy)
      set(hText,'String',num2str(n));
   end
end

Contact us