root/OptimizingMATLABCode/GPUWave/WaveCPU.m @ 11
| 10 | anderm8 | function vv = WaveCPU(N, Nsteps, hTopAxes, hBottomAxes, hStartStopBtn, ...
|
|
hIterationText, hElapsedTimeText, plotStep, ticT)
|
|||
%% 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 Chebysehv spectral method in space (using FFT).
|
|||
%
|
|||
% The code has been modified from an example in Spectral Methods in MATLAB
|
|||
% by Trefethen, Lloyd N.
|
|||
%
|
|||
% Copyright 2015 The MathWorks, Inc.
|
|||
if nargin == 2
|
|||
guiMode = false;
|
|||
else
|
|||
guiMode = true;
|
|||
on_state = get(hStartStopBtn,'Max');
|
|||
end
|
|||
% Points in X and Y
|
|||
x = cos(pi*(0:N)/N); % using Chebyshev points
|
|||
y = x';
|
|||
% Calculating time step
|
|||
dt = 6/N^2;
|
|||
% Setting up grid
|
|||
[xx,yy] = meshgrid(x,y);
|
|||
% Calculate initial values
|
|||
vv = exp(-40*((xx-.4).^2 + yy.^2));
|
|||
vvold = vv;
|
|||
ii = 2:N;
|
|||
index1 = 1i*[0:N-1 0 1-N:-1];
|
|||
index2 = -[0:N 1-N:-1].^2;
|
|||
% Weights used for spectral differentiation via FFT
|
|||
W1T = repmat(index1,N-1,1);
|
|||
W2T = repmat(index2,N-1,1);
|
|||
W3T = repmat(index1.',1,N-1);
|
|||
W4T = repmat(index2.',1,N-1);
|
|||
WuxxT1 = repmat((1./(1-x(ii).^2)),N-1,1);
|
|||
WuxxT2 = repmat(x(ii)./(1-x(ii).^2).^(3/2),N-1,1);
|
|||
WuyyT1 = repmat(1./(1-y(ii).^2),1,N-1);
|
|||
WuyyT2 = repmat(y(ii)./(1-y(ii).^2).^(3/2),1,N-1);
|
|||
% Start time-stepping
|
|||
n = 1;
|
|||
while n <= Nsteps
|
|||
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))).';
|
|||
% Calculating 2nd derivative in x
|
|||
uxx(ii,ii) = W2(:,ii).* WuxxT1 - W1(:,ii).*WuxxT2;
|
|||
uxx([1,N+1],[1,N+1]) = 0;
|
|||
V = [vv(:,ii); vv((N:-1:2),ii)];
|
|||
U = real(fft(V));
|
|||
W1 = real(ifft(U.*W3T));
|
|||
W2 = real(ifft(U.*W4T));
|
|||
% Calculating 2nd derivative in y
|
|||
uyy(ii,ii) = W2(ii,:).* WuyyT1 - W1(ii,:).*WuyyT2;
|
|||
uyy([1,N+1],[1,N+1]) = 0;
|
|||
% Computing new value using 2nd order central finite difference in
|
|||
% time
|
|||
vvnew = 2*vv - vvold + dt*dt*(uxx+uyy);
|
|||
vvold = vv;
|
|||
vv = vvnew;
|
|||
if guiMode
|
|||
if ~isequal(get(hStartStopBtn, 'Value'), on_state)
|
|||
break
|
|||
end
|
|||
if n == 1 || mod(n, plotStep) == 0
|
|||
updateUI(xx, yy, vv, vvold, n, N, hTopAxes, ...
|
|||
hBottomAxes, hIterationText, hElapsedTimeText, ticT);
|
|||
end
|
|||
end
|
|||
n = n + 1;
|
|||
end
|