From 3c7fa1ea20242e883028cba9a5eb7010d01279a0 Mon Sep 17 00:00:00 2001 From: Bradley Treeby Date: Tue, 21 Apr 2020 09:14:12 +0100 Subject: [PATCH] Version 1.1 (see change log) --- README.md | 17 +- examples/example_gradient_calculation.m | 22 +- examples/example_transform_wsws_sequence.m | 4 +- ...ve_eq_pstd.m => example_wave_eq_pstd_1D.m} | 90 ++-- .../example_wave_eq_pstd_1D_non_reflecting.m | 284 +++++++++++++ examples/example_wave_eq_pstd_2D_dirichlet.m | 217 ++++++++++ examples/example_wave_eq_pstd_2D_neumann.m | 216 ++++++++++ examples/example_wsws_gradient.m | 113 +++++ gradientDTT1D.m => gradientDtt1D.m | 397 +++++++++--------- 9 files changed, 1132 insertions(+), 228 deletions(-) rename examples/{example_wave_eq_pstd.m => example_wave_eq_pstd_1D.m} (73%) create mode 100644 examples/example_wave_eq_pstd_1D_non_reflecting.m create mode 100644 examples/example_wave_eq_pstd_2D_dirichlet.m create mode 100644 examples/example_wave_eq_pstd_2D_neumann.m create mode 100644 examples/example_wsws_gradient.m rename gradientDTT1D.m => gradientDtt1D.m (53%) diff --git a/README.md b/README.md index 7d2c196..8deabd4 100644 --- a/README.md +++ b/README.md @@ -1,4 +1,6 @@ # MATLAB Discrete Trigonometric Transform Library +[![View MATLAB Discrete Trigonometric Transform Library on File Exchange](https://www.mathworks.com/matlabcentral/images/matlab-file-exchange.svg)](https://www.mathworks.com/matlabcentral/fileexchange/75071-matlab-discrete-trigonometric-transform-library) + ## Author This library is written by Bradley Treeby, University College London. Contact: b.treeby@ucl.ac.uk @@ -31,4 +33,17 @@ This program is free software: you can redistribute it and/or modify it under th This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -You should have received a copy of the GNU General Public License along with this program. If not, see . \ No newline at end of file +You should have received a copy of the GNU General Public License along with this program. If not, see . + +## Change Log + +* v1.1 (21 April 2020): + * Fixed bug in `gradientDtt1D` (missing `numDim` function) + * Added `align_output` option to `gradientDtt1D` + * Updated `example_wave_eq_pstd_1D` to enforce correct boundary position + * Added `example_wave_eq_pstd_1D_non_reflecting` to simulate general boundary conditions + * Added `example_wsws_gradient` to illustrate grid shifting with and without `align_output` +* Added `example_wave_eq_pstd_2D_neumann` and `example_wave_eq_pstd_2D_dirichlet` + +* v1.0 (17 April 2020): + * Initial release \ No newline at end of file diff --git a/examples/example_gradient_calculation.m b/examples/example_gradient_calculation.m index c7a0ec8..66b5439 100644 --- a/examples/example_gradient_calculation.m +++ b/examples/example_gradient_calculation.m @@ -1,9 +1,9 @@ % DESCRIPTION: % This example script creates different periodic sequences of known % symmetry, and then numerically calculates the gradients using -% gradientDTT1D. The numerical gradient is then compared with the known +% gradientDtt1D. The numerical gradient is then compared with the known % analytical gradient. A test sequence is created for all supported DTT -% symmetries. The grid shifting functionality of gradientDTT1D can also +% symmetries. The grid shifting functionality of gradientDtt1D can also % be tested by changing the value of shift. % % ABOUT: @@ -13,7 +13,7 @@ % % Copyright (C) 2013-2020 Bradley Treeby % -% See also dtt1D, gradientDTT1D +% See also dtt1D, gradientDtt1D % This program is free software: you can redistribute it and/or modify it % under the terms of the GNU General Public License as published by the @@ -68,7 +68,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = -sin(x_out); % plot functions @@ -114,7 +114,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = -sin(x_out); % plot functions @@ -160,7 +160,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = -sin(x_out); % plot functions @@ -206,7 +206,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = -sin(x_out); % plot functions @@ -252,7 +252,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = cos(x_out); % plot functions @@ -298,7 +298,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = cos(x_out); % plot functions @@ -344,7 +344,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = cos(x_out); % plot functions @@ -390,7 +390,7 @@ end % compute gradient - dfdx_num = gradientDTT1D(f, dx, test_case, shift); + dfdx_num = gradientDtt1D(f, dx, test_case, shift); dfdx_analytical = cos(x_out); % plot functions diff --git a/examples/example_transform_wsws_sequence.m b/examples/example_transform_wsws_sequence.m index 5735e12..827b82e 100644 --- a/examples/example_transform_wsws_sequence.m +++ b/examples/example_transform_wsws_sequence.m @@ -16,9 +16,9 @@ % ABOUT: % author - Bradley Treeby % date - 31 March 2020 -% last update - 2 April 2020 +% last update - 20 April 2020 % -% See also dtt1D, gradientDTT1D +% See also dtt1D, gradientDtt1D % This program is free software: you can redistribute it and/or modify it % under the terms of the GNU General Public License as published by the diff --git a/examples/example_wave_eq_pstd.m b/examples/example_wave_eq_pstd_1D.m similarity index 73% rename from examples/example_wave_eq_pstd.m rename to examples/example_wave_eq_pstd_1D.m index e2beabd..55ef445 100644 --- a/examples/example_wave_eq_pstd.m +++ b/examples/example_wave_eq_pstd_1D.m @@ -1,16 +1,22 @@ % DESCRIPTION: % This example script solves the 1D wave equation (written as two % coupled first-order equations) using a DTT-based PSTD method subject -% to different combinations of Dirichlet and Neumann boundary -% conditions at each end of the domain. The spatial gradients are -% calculated using discrete cosine transforms (DCTs) and discrete sine -% transforms (DST) chosen to enfore the required boundary condition. -% Time integration is performed using a first-order accurate foward -% difference scheme. A classical Fourier PSTD method is also used with -% periodic bounday conditions for comparison. Further details are given -% in [1]. +% to four different combinations of Dirichlet and Neumann boundary +% conditions at each end of the domain. % -% This script reproduces Figure 5 of [1]. +% The spatial gradients are calculated using discrete cosine transforms +% (DCTs) and discrete sine transforms (DST) chosen to enfore the +% required boundary condition. Time integration is performed using a +% first-order accurate foward difference scheme. A classical Fourier +% PSTD method is also used with periodic bounday conditions for +% comparison. +% +% For the DTT-based simulations, the position of the boundary is +% assumed to be at the first and last grid points. To account for the +% different implied symmetries for the different transform types, the +% calculations are performed on grids of different sizes. +% +% Further details are given in [1]. % % [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral % time-domain (PSTD) methods for the wave equation: Realising boundary @@ -19,11 +25,11 @@ % ABOUT: % author - Bradley Treeby % date - 31 March 2020 -% last update - 1 April 2020 +% last update - 20 April 2020 % % Copyright (C) 2020 Bradley Treeby % -% See also dtt1D, gradientDTT1D +% See also dtt1D, gradientDtt1D % This program is free software: you can redistribute it and/or modify it % under the terms of the GNU General Public License as published by the @@ -44,11 +50,12 @@ % DEFINE LITERALS % ========================================================================= -% shift variables used by gradientDTT1D -shift_pos = 1; -shift_neg = 2; +% shift and align variables used by gradientDtt1D (shift = 1 and shift = 2 +% return the same output if align_output = false) +shift_output = 1; +align_output = false; -% transform types used by gradientDTT1D +% transform types used by gradientDtt1D DCT1 = 1; % WSWS DCT2 = 2; % HSHS DCT3 = 3; % WSWA @@ -118,6 +125,7 @@ % select figure figure(plot_fig); + hold on; % preallocate matrix to store data for waterfall plot p_waterfall = zeros(waterfall_snapshots, Nx); @@ -132,57 +140,86 @@ % Neumann-Neumann / WSWS T1 = DCT1; T2 = DST2; + + % set indices for representative sample + ind1 = 1; + ind2 = Nx; + + % set title waterfall_title = 'Neumann-Neumann (WSWS)'; + title(waterfall_title); case 2 % Neumann-Dirichlet / WSWA T1 = DCT3; T2 = DST4; + + % set indices for representative sample + ind1 = 1; + ind2 = Nx - 1; + + % set title waterfall_title = 'Neumann-Dirichlet (WSWA)'; + title(waterfall_title); case 3 % Dirichlet-Neumann / WAWS T1 = DST3; T2 = DCT4; + + % set indices for representative sample + ind1 = 2; + ind2 = Nx; + + % set title waterfall_title = 'Dirichlet-Neumann (WAWS)'; + title(waterfall_title); case 4 - % Neumann-Neumann / WSWS + % Dirichlet-Dirichlet / WAWA T1 = DST1; T2 = DCT2; + + % set indices for representative sample + ind1 = 2; + ind2 = Nx - 1; + + % set title waterfall_title = 'Dirichlet-Dirichlet (WAWA)'; + title(waterfall_title); end - % assign initial condition for the pressure - p = p0; + % assign initial condition for p, just selecting representative sample + p = p0(ind1:ind2); % run the model backwards for dt/2 to calculate the initial condition % for u, to take into account the time staggering between p and u - u = u0 + (dt / 2) / rho0 * gradientDTT1D(p0, dx, T1, shift_pos); + u = u0(1:end - 1) + (dt / 2) / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output); % calculate pressure in a loop for time_ind = 1:Nt % update u - u = u - dt / rho0 * gradientDTT1D(p, dx, T1, shift_pos); + u = u - dt / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output); % update p - p = p - dt * rho0 * c0^2 * gradientDTT1D(u, dx, T2, shift_neg); + p = p - dt * rho0 * c0^2 * gradientDtt1D(u, dx, T2, shift_output, align_output); % plot pressure field if ~rem(time_ind, plot_freq) - plot(x, p, 'k-'); + cla; + plot(x(ind1:ind2), p, 'k-'); set(gca, 'YLim', [-1, 1]); drawnow; end % save the pressure field if ~rem(time_ind, Nt / waterfall_snapshots) - p_waterfall(waterfall_ind, :) = p; + p_waterfall(waterfall_ind, ind1:ind2) = p; waterfall_ind = waterfall_ind + 1; end @@ -209,6 +246,10 @@ % select figure figure(plot_fig); +% set title +waterfall_title = 'Periodic'; +title(waterfall_title); + % preallocate matrix to store data for waterfall plot p_waterfall = zeros(waterfall_snapshots, Nx); @@ -242,6 +283,7 @@ % plot pressure field if ~rem(time_ind, plot_freq) + cla; plot(x, p, 'k-'); set(gca, 'YLim', [-1, 1]); drawnow; @@ -264,7 +306,7 @@ ylabel('time'); zlabel('pressure'); xlabel('position'); -title('Periodic'); +title(waterfall_title); grid off % close animation figure diff --git a/examples/example_wave_eq_pstd_1D_non_reflecting.m b/examples/example_wave_eq_pstd_1D_non_reflecting.m new file mode 100644 index 0000000..d94b433 --- /dev/null +++ b/examples/example_wave_eq_pstd_1D_non_reflecting.m @@ -0,0 +1,284 @@ +% DESCRIPTION: +% This example script solves the 1D wave equation (written as two +% coupled first-order equations) using a DTT-based PSTD method subject +% to four different combinations of Dirichlet and Neumann boundary +% conditions at each end of the domain. +% +% The solutions for different boundary conditions are then summed to +% give different reflection coefficients at each end of the domain, +% including non-reflecting boundaries. This builds on +% example_wave_eq_pstd.m. +% +% Further details are given in [1]. +% +% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral +% time-domain (PSTD) methods for the wave equation: Realising boundary +% conditions with discrete sine and cosine transforms", 2020. +% +% ABOUT: +% author - Bradley Treeby +% date - 19 April 2020 +% last update - 20 April 2020 +% +% Copyright (C) 2020 Bradley Treeby +% +% See also dtt1D, gradientDtt1D + +% This program is free software: you can redistribute it and/or modify it +% under the terms of the GNU General Public License as published by the +% Free Software Foundation, either version 3 of the License, or (at your +% option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +% Public License for more details. +% +% You should have received a copy of the GNU General Public License along +% with this program. If not, see . + +clearvars; + +%% ======================================================================== +% DEFINE LITERALS +% ========================================================================= + +% shift and align variables used by gradientDtt1D (shift = 1 and shift = 2 +% return the same output if align_output = false) +shift_output = 1; +align_output = false; + +% transform types used by gradientDtt1D +DCT1 = 1; % WSWS +DCT2 = 2; % HSHS +DCT3 = 3; % WSWA +DCT4 = 4; % HSHA +DST1 = 5; % WAWA +DST2 = 6; % HAHA +DST3 = 7; % WAWS +DST4 = 8; % HAHS + +% number of time snapshots to include in waterfall plot +waterfall_snapshots = 80; + +% number of boundary conditions +bc_num = 4; + +%% ======================================================================== +% DEFINE SIMULATION SETTINGS +% ========================================================================= + +% set the grid size +Nx = 256; % grid size [m] +dx = 1/Nx; % grid spacing [m] + +% set the medium properties +c0 = 1500; % sound speed [m/s] +rho0 = 1000; % density [kg/m^3] + +% set CFL and number of time steps +CFL = 0.2; +Nt = 2720; + +%% ======================================================================== +% DEFINE INITIAL CONDITIONS +% ========================================================================= + +% spatial grid +x = (0:Nx - 1) * dx; + +% properties of Gaussian initial condition +width = Nx * dx / 28; +offset = Nx * dx / 3; + +% define initial pressure distribution on regular grid as a Gaussian +p0 = exp(-((x - offset) / width).^2); + +%% ======================================================================== +% RUN SIMULATIONS USING DTT-BASED PSTD METHOD +% ========================================================================= + +% calculate the time step +dt = CFL * dx / c0; + +% normalised plot axes +x_ax = (0:Nx - 1) / (Nx - 1); +t_ax = (0:waterfall_snapshots - 1) / (waterfall_snapshots - 1); + +% preallocate matrix to store simulation data +p_four_bc = zeros(bc_num, Nx, waterfall_snapshots); + +% loop over boundary conditions +for bc_ind = 1:bc_num + + % initialise index for storing waterfall data + waterfall_ind = 1; + + % define which transforms to use + switch bc_ind + case 1 + + % Neumann-Neumann / WSWS + T1 = DCT1; + T2 = DST2; + + % set indices for representative sample + ind1 = 1; + ind2 = Nx; + + case 2 + + % Neumann-Dirichlet / WSWA + T1 = DCT3; + T2 = DST4; + + % set indices for representative sample + ind1 = 1; + ind2 = Nx - 1; + + case 3 + + % Dirichlet-Neumann / WAWS + T1 = DST3; + T2 = DCT4; + + % set indices for representative sample + ind1 = 2; + ind2 = Nx; + + case 4 + + % Dirichlet-Dirichlet / WAWA + T1 = DST1; + T2 = DCT2; + + % set indices for representative sample + ind1 = 2; + ind2 = Nx - 1; + + end + + % assign initial condition for p, just selecting representative sample + p = p0(ind1:ind2); + + % run the model backwards for dt/2 to calculate the initial condition + % for u, to take into account the time staggering between p and u, in + % this case setting u0 = 0 + u = (dt / 2) / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output); + + % calculate pressure in a loop + for time_ind = 1:Nt + + % update u + u = u - dt / rho0 * gradientDtt1D(p, dx, T1, shift_output, align_output); + + % update p + p = p - dt * rho0 * c0^2 * gradientDtt1D(u, dx, T2, shift_output, align_output); + + % save the pressure field + if ~rem(time_ind, floor(Nt / waterfall_snapshots)) + p_four_bc(bc_ind, ind1:ind2, waterfall_ind) = p; + waterfall_ind = waterfall_ind + 1; + end + + end + +end + +%% ======================================================================== +% COMBINE SIMULATIONS WITH DIFFERENT BC +% ========================================================================= + +% form matrix of the different possible boundary conditions, where +1 +% corresponds to a positive image source, and -1 to a negative image +% source, thus the four columns correspond to Neumann-Neumann, +% Neumann-Dirichlet, Dirichlet-Neumann, and Dirichlet-Dirichlet +M = [1 1 -1 -1; + 1 1 1 1; + 1 -1 1 -1; + 1 -1 -1 1]; + + % open a new figure window +plot_fig = figure; + +% loop over different reflection coefficients +for ind = 1:2 + + % set the desired reflection coefficient + switch ind + case 1 + + % non-reflecting boundaries + R_l = 0; + R_r = 0; + + case 2 + + % partially reflecting boundary + R_l = 0; + R_r = 0.5; + + end + + % form r vector (this corresponds to the image source amplitudes) + r = [R_l, 1, R_r, R_r * R_l].'; + + % solve for the weights for each of the pre-calculated fields + w = 1/4 * M.' * r %#ok + + % preallocate matrix + p_waterfall = zeros(Nx, waterfall_snapshots); + + % loop over time (waterfall) index + for waterfall_ind = 1:waterfall_snapshots + + % form the pressure field by summing weighted fields + p_waterfall(:, waterfall_ind) = w.' * p_four_bc(:, :, waterfall_ind); + + % plot the fields and the summation + figure(plot_fig); + + subplot(5, 1, 1); + plot(x_ax, p_four_bc(1, :, waterfall_ind), 'k-'); + set(gca, 'YLim', [-1, 1]); + title('Neumann-Neumann (WSWS)'); + + subplot(5, 1, 2); + plot(x_ax, p_four_bc(2, :, waterfall_ind), 'k-'); + set(gca, 'YLim', [-1, 1]); + title('Neumann-Dirichlet (WSWA)'); + + subplot(5, 1, 3); + plot(x_ax, p_four_bc(3, :, waterfall_ind), 'k-'); + set(gca, 'YLim', [-1, 1]); + title('Dirichlet-Neumann (WAWS)'); + + subplot(5, 1, 4); + plot(x_ax, p_four_bc(4, :, waterfall_ind), 'k-'); + set(gca, 'YLim', [-1, 1]); + title('Dirichlet-Dirichlet (WAWA)'); + + subplot(5, 1, 5); + plot(x_ax, p_waterfall(:, waterfall_ind), 'k-'); + set(gca, 'YLim', [-1, 1]); + title(['Sum r = [' num2str(r(1)) ', ' num2str(r(2)) ', ' num2str(r(3)) ', ' num2str(r(4)) ']']); + drawnow; + + end + + % waterfall plot of evolution of field + figure; + waterfall(x_ax, t_ax, p_waterfall.'); + view(10, 70); + colormap([0, 0, 0]); + set(gca, 'ZLim', [-1, 1], 'FontSize', 12); + ylabel('time'); + zlabel('pressure'); + xlabel('position'); + title(['r = [' num2str(r(1)) ', ' num2str(r(2)) ', ' num2str(r(3)) ', ' num2str(r(4)) ']']); + grid off; + +end + +% close animation figure +close(plot_fig); \ No newline at end of file diff --git a/examples/example_wave_eq_pstd_2D_dirichlet.m b/examples/example_wave_eq_pstd_2D_dirichlet.m new file mode 100644 index 0000000..30d7565 --- /dev/null +++ b/examples/example_wave_eq_pstd_2D_dirichlet.m @@ -0,0 +1,217 @@ +% DESCRIPTION: +% This example script solves the 2D wave equation (written as two +% coupled first-order equations) using a DTT-based PSTD method subject +% to Dirichlet boundary conditions on each side of the domain. +% +% The spatial gradients are calculated using discrete cosine transforms +% (DCTs) and discrete sine transforms (DST) chosen to enfore the +% required boundary condition. Time integration is performed using a +% first-order accurate foward difference scheme. +% +% Due to the symmetry of the DTTs used (WAWA), the position of the +% boundary is one point outside the first and last grid points in the +% domain. +% +% Further details are given in [1]. +% +% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral +% time-domain (PSTD) methods for the wave equation: Realising boundary +% conditions with discrete sine and cosine transforms", 2020. +% +% ABOUT: +% author - Bradley Treeby +% date - 21 April 2020 +% last update - 21 April 2020 +% +% Copyright (C) 2020 Bradley Treeby +% +% See also dtt1D, gradientDtt1D + +% This program is free software: you can redistribute it and/or modify it +% under the terms of the GNU General Public License as published by the +% Free Software Foundation, either version 3 of the License, or (at your +% option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +% Public License for more details. +% +% You should have received a copy of the GNU General Public License along +% with this program. If not, see . + +clearvars; + +% ========================================================================= +% DEFINE LITERALS +% ========================================================================= + +% transform types used by dtt1D +DCT1 = 1; % WSWS +DCT2 = 2; % HSHS +DCT3 = 3; % WSWA +DCT4 = 4; % HSHA +DST1 = 5; % WAWA +DST2 = 6; % HAHA +DST3 = 7; % WAWS +DST4 = 8; % HAHS + +% plot frequency for snapshots of the field (time steps) +plot_freq = 250; + +% ========================================================================= +% DEFINE SIMULATION SETTINGS +% ========================================================================= + +% set the grid size (assuming grid is square) +Nx = 256; % grid size [m] +dx = 1/Nx; % grid spacing [m] + +% set the medium properties +c0 = 1500; % sound speed [m/s] +rho0 = 1000; % density [kg/m^3] + +% set CFL and number of time steps +CFL = 0.2; +Nt = 1000; + +% ========================================================================= +% DEFINE INITIAL CONDITIONS +% ========================================================================= + +% spatial grid +x = (0:Nx - 1) * dx; + +% properties of Gaussian initial condition +width = Nx * dx / 28; +offset = Nx * dx / 3; + +% define initial pressure distribution on regular grid as a Gaussian +p = 0.5 * exp(-((x - offset) / width).^2); +p = p.' * p; + +% define initial velocity to be zero +ux = 0; +uy = 0; + +% ========================================================================= +% DEFINE DTT VARIABLES +% ========================================================================= + +% compute the implied period of the input function +M = 2 .* (Nx + 1); + +% define wavenumber vectors for WAWA / DST-I +k_wawa = 2 .* pi .* (1:1:(M/2 - 1)).' ./ (M .* dx); + +% define wavenumber vectors for HSHS / DCT-II +k_hshs = 2 .* pi .* (0:1:(M/2 - 1)).' ./ (M .* dx); + +% ========================================================================= +% RUN SIMULATIONS USING DTT-BASED PSTD METHOD +% ========================================================================= + +% calculate the time step +dt = CFL * dx / c0; + +% create custom colour map +cm = [bone(256); flipud(hot(256))]; + +% calculate pressure in a loop +for time_ind = 1:Nt + + % --------------- + % calculate dpdx + % --------------- + + % forward transform over x + dpdx = dtt1D(p, DST1, 1); + + % spectral derivative + dpdx = bsxfun(@times, k_wawa, dpdx); + + % WAWA -> HSHS, prepend zero + dpdx = [zeros(1, Nx); dpdx]; %#ok + + % inverse transform over x, C2^-1 = C3 + dpdx = dtt1D(dpdx, DCT3, 1) ./ M; + + % --------------- + % calculate dpdy + % --------------- + + % forward transform over y + dpdy = dtt1D(p, DST1, 2); + + % spectral derivative + dpdy = bsxfun(@times, k_wawa.', dpdy); + + % WAWA -> HSHS, prepend zero + dpdy = [zeros(Nx, 1), dpdy]; %#ok + + % inverse transform over y, C2^-1 = C3 + dpdy = dtt1D(dpdy, DCT3, 2) ./ M; + + % --------------- + % update u + % --------------- + + % update components of particle velocity + ux = ux - dt / rho0 * dpdx; + uy = uy - dt / rho0 * dpdy; + + % --------------- + % calculate duxdx + % --------------- + + % forward transform over x + duxdx = dtt1D(ux, DCT2, 1); + + % spectral derivative + duxdx = bsxfun(@times, -k_hshs, duxdx); + + % HSHS -> WAWA, remove left endpoint + duxdx = duxdx(2:end, :); + + % inverse transform over x, S1^-1 = S1 + duxdx = dtt1D(duxdx, DST1, 1) ./ M; + + % --------------- + % calculate duydy + % --------------- + + % forward transform over y + duydy = dtt1D(uy, DCT2, 2); + + % spectral derivative + duydy = bsxfun(@times, -k_hshs.', duydy); + + % HSHS -> WAWA, remove left endpoint + duydy = duydy(:, 2:end); + + % inverse transform over y, S1^-1 = S1 + duydy = dtt1D(duydy, DST1, 2) ./ M; + + % --------------- + % update p + % --------------- + + % update p + p = p - dt * rho0 * c0^2 * (duxdx + duydy); + + % --------------- + % plot + % --------------- + + % plot snapshots of pressure field + if (time_ind == 1) || ~rem(time_ind, plot_freq) + figure; + imagesc(x, x, p, 0.075 * [-1, 1]); + axis image; + colormap(cm); + box on; + set(gca, 'XTick', 0:0.25:1, 'YTick', 0:0.25:1, 'XTickLabel', {}, 'YTickLabel', {}); + drawnow; + end + +end \ No newline at end of file diff --git a/examples/example_wave_eq_pstd_2D_neumann.m b/examples/example_wave_eq_pstd_2D_neumann.m new file mode 100644 index 0000000..b6f987f --- /dev/null +++ b/examples/example_wave_eq_pstd_2D_neumann.m @@ -0,0 +1,216 @@ +% DESCRIPTION: +% This example script solves the 2D wave equation (written as two +% coupled first-order equations) using a DTT-based PSTD method subject +% to Neumann boundary conditions on each side of the domain. +% +% The spatial gradients are calculated using discrete cosine transforms +% (DCTs) and discrete sine transforms (DST) chosen to enfore the +% required boundary condition. Time integration is performed using a +% first-order accurate foward difference scheme. +% +% Due to the symmetry of the DTTs used (WSWS), the position of the +% boundary is at the first and last grid points in the domain. +% +% Further details are given in [1]. +% +% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral +% time-domain (PSTD) methods for the wave equation: Realising boundary +% conditions with discrete sine and cosine transforms", 2020. +% +% ABOUT: +% author - Bradley Treeby +% date - 20 April 2020 +% last update - 21 April 2020 +% +% Copyright (C) 2020 Bradley Treeby +% +% See also dtt1D, gradientDtt1D + +% This program is free software: you can redistribute it and/or modify it +% under the terms of the GNU General Public License as published by the +% Free Software Foundation, either version 3 of the License, or (at your +% option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +% Public License for more details. +% +% You should have received a copy of the GNU General Public License along +% with this program. If not, see . + +clearvars; + +% ========================================================================= +% DEFINE LITERALS +% ========================================================================= + +% transform types used by dtt1D +DCT1 = 1; % WSWS +DCT2 = 2; % HSHS +DCT3 = 3; % WSWA +DCT4 = 4; % HSHA +DST1 = 5; % WAWA +DST2 = 6; % HAHA +DST3 = 7; % WAWS +DST4 = 8; % HAHS + +% plot frequency for snapshots of the field (time steps) +plot_freq = 250; + +% ========================================================================= +% DEFINE SIMULATION SETTINGS +% ========================================================================= + +% set the grid size (assuming grid is square) +Nx = 256; % grid size [m] +dx = 1/Nx; % grid spacing [m] + +% set the medium properties +c0 = 1500; % sound speed [m/s] +rho0 = 1000; % density [kg/m^3] + +% set CFL and number of time steps +CFL = 0.2; +Nt = 1000; + +% ========================================================================= +% DEFINE INITIAL CONDITIONS +% ========================================================================= + +% spatial grid +x = (0:Nx - 1) * dx; + +% properties of Gaussian initial condition +width = Nx * dx / 28; +offset = Nx * dx / 3; + +% define initial pressure distribution on regular grid as a Gaussian +p = 0.5 * exp(-((x - offset) / width).^2); +p = p.' * p; + +% define initial velocity to be zero +ux = 0; +uy = 0; + +% ========================================================================= +% DEFINE DTT VARIABLES +% ========================================================================= + +% compute the implied period of the input function +M = 2 .* (Nx - 1); + +% define wavenumber vectors for WSWS / DCT-I +k_wsws = 2 .* pi .* (0:1:M/2).' ./ (M .* dx); + +% define wavenumber vectors for HAHA / DST-II +k_haha = 2 .* pi .* (1:1:M/2).' ./ (M .* dx); + +% ========================================================================= +% RUN SIMULATIONS USING DTT-BASED PSTD METHOD +% ========================================================================= + +% calculate the time step +dt = CFL * dx / c0; + +% create custom colour map +cm = [bone(256); flipud(hot(256))]; + +% calculate pressure in a loop +for time_ind = 1:Nt + + % --------------- + % calculate dpdx + % --------------- + + % forward transform over x + dpdx = dtt1D(p, DCT1, 1); + + % spectral derivative + dpdx = bsxfun(@times, -k_wsws, dpdx); + + % WSWS -> HAHA, remove left endpoint + dpdx = dpdx(2:end, :); + + % inverse transform over x, S2^-1 = S3 + dpdx = dtt1D(dpdx, DST3, 1) ./ M; + + % --------------- + % calculate dpdy + % --------------- + + % forward transform over y + dpdy = dtt1D(p, DCT1, 2); + + % spectral derivative + dpdy = bsxfun(@times, -k_wsws.', dpdy); + + % WSWS -> HAHA, remove left endpoint + dpdy = dpdy(:, 2:end); + + % inverse transform over y, S2^-1 = S3 + dpdy = dtt1D(dpdy, DST3, 2) ./ M; + + % --------------- + % update u + % --------------- + + % update components of particle velocity + ux = ux - dt / rho0 * dpdx; + uy = uy - dt / rho0 * dpdy; + + % --------------- + % calculate duxdx + % --------------- + + % forward transform over x + duxdx = dtt1D(ux, DST2, 1); + + % spectral derivative + duxdx = bsxfun(@times, k_haha, duxdx); + + % HAHA -> WSWS, prepend zero + duxdx = [zeros(1, Nx); duxdx]; %#ok + + % inverse transform over x, C1^-1 = C1 + duxdx = dtt1D(duxdx, DCT1, 1) ./ M; + + % --------------- + % calculate duydy + % --------------- + + % forward transform over y + duydy = dtt1D(uy, DST2, 2); + + % spectral derivative + duydy = bsxfun(@times, k_haha.', duydy); + + % HAHA -> WSWS, prepend zero + duydy = [zeros(Nx, 1), duydy]; %#ok + + % inverse transform over y, C1^-1 = C1 + duydy = dtt1D(duydy, DCT1, 2) ./ M; + + % --------------- + % update p + % --------------- + + % update p + p = p - dt * rho0 * c0^2 * (duxdx + duydy); + + % --------------- + % plot + % --------------- + + % plot snapshots of pressure field + if (time_ind == 1) || ~rem(time_ind, plot_freq) + figure; + imagesc(x, x, p, 0.075 * [-1, 1]); + axis image; + colormap(cm); + box on; + set(gca, 'XTick', 0:0.25:1, 'YTick', 0:0.25:1, 'XTickLabel', {}, 'YTickLabel', {}); + drawnow; + end + +end \ No newline at end of file diff --git a/examples/example_wsws_gradient.m b/examples/example_wsws_gradient.m new file mode 100644 index 0000000..e7bf125 --- /dev/null +++ b/examples/example_wsws_gradient.m @@ -0,0 +1,113 @@ +% DESCRIPTION: +% This example script shows a WSWS gradient calculated on regular and +% staggered grids, both with and without alignment. +% +% Further details are given in [1]. +% +% [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral +% time-domain (PSTD) methods for the wave equation: Realising boundary +% conditions with discrete sine and cosine transforms", 2020. +% +% ABOUT: +% author - Bradley Treeby +% date - 20 April 2020 +% last update - 20 April 2020 +% +% Copyright (C) 2020-2020 Bradley Treeby +% +% See also dtt1D, gradientDtt1D + +% This program is free software: you can redistribute it and/or modify it +% under the terms of the GNU General Public License as published by the +% Free Software Foundation, either version 3 of the License, or (at your +% option) any later version. +% +% This program is distributed in the hope that it will be useful, but +% WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General +% Public License for more details. +% +% You should have received a copy of the GNU General Public License along +% with this program. If not, see . + +clearvars; + +% create function with WSWS symmetry +Nx = 10; +dx = pi/(Nx-1); +x_in = 0:dx:pi; +f = cos(x_in); + +% ========================================================================= +% GRADIENT WITH ALIGNMENT +% ========================================================================= + +% define x vectors +x_out_def = x_in; +x_out_pos = x_in + dx/2; +x_out_neg = x_in - dx/2; + +% compute gradient using DCT-I (WSWS symmetry) +dfdx_def = gradientDtt1D(f, dx, 1, 0); +dfdx_pos = gradientDtt1D(f, dx, 1, 1); +dfdx_neg = gradientDtt1D(f, dx, 1, 2); + +% analytical derivative +x_an = -dx:dx/5:2*pi - dx; +dfdx_an = -sin(x_an); + +% plot and compare +figure; +plot(... + x_an, dfdx_an, 'k-', ... + x_out_def, dfdx_def, 'r.', ... + x_out_pos, dfdx_pos, 'bo', ... + x_out_neg, dfdx_neg, 'g.'); +legend(... + 'Analytical Derivative', ... + 'DTT shift = 0', ... + 'DTT shift = 1', ... + 'DTT shift = 2', ... + 'Location', 'Best'); +xlabel('x'); +ylabel('f''(x)'); +title('DCT-I WSWS (align output = true)'); +set(gca, 'XLim', [-dx, 2*pi - dx]); +grid on; + +% ========================================================================= +% GRADIENT WITHOUT ALIGNMENT +% ========================================================================= + +% define x vectors +x_out_def = x_in(2:end-1); +x_out_pos = x_in(1:end-1) + dx/2; +x_out_neg = x_in(1:end-1) + dx/2; + +% compute gradient using DCT-I (WSWS symmetry) +dfdx_def = gradientDtt1D(f, dx, 1, 0, false); +dfdx_pos = gradientDtt1D(f, dx, 1, 1, false); +dfdx_neg = gradientDtt1D(f, dx, 1, 2, false); + +% analytical derivative +x_an = -dx:dx/5:2*pi - dx; +dfdx_an = -sin(x_an); + +% plot and compare +figure; +plot(... + x_an, dfdx_an, 'k-', ... + x_out_def, dfdx_def, 'r.', ... + x_out_pos, dfdx_pos, 'bo', ... + x_out_neg, dfdx_neg, 'g.'); +legend(... + 'Analytical Derivative', ... + 'DTT shift = 0', ... + 'DTT shift = 1', ... + 'DTT shift = 2', ... + 'Location', 'Best'); +xlabel('x'); +ylabel('f''(x)'); +title('DCT-I WSWS (align output = false)'); +set(gca, 'XLim', [-dx, 2*pi - dx]); +grid on; \ No newline at end of file diff --git a/gradientDTT1D.m b/gradientDtt1D.m similarity index 53% rename from gradientDTT1D.m rename to gradientDtt1D.m index c06d2aa..49fa19e 100644 --- a/gradientDTT1D.m +++ b/gradientDtt1D.m @@ -1,8 +1,8 @@ -function deriv = gradientDTT1D(f, dx, dtt_type, shift) +function deriv = gradientDtt1D(f, dx, dtt_type, shift, align_output) %GRADIENTDTT1D Calculate gradient using discrete trigonometric transforms. % % DESCRIPTION: -% gradientDTT1D computes the spectral gradient of a 1D input function +% gradientDtt1D computes the spectral gradient of a 1D input function % using discrete trigonometric transforms (DTTs). The DTTs used to % transform the input to and from the frequency domain are chosen based % on the symmetry of the input f, defined using dtt_type. The gradient @@ -14,47 +14,55 @@ % need to be defined). This means that in some cases the basis % functions weights are trimmed or appended after the forward % transform. To ensure the output of gradientDTT1D is the same length -% as the input f, additional values are appended or removed from the -% calculated gradient after the inverse transform is calculated. These -% values are known from the symmetry of the output sequence. +% as the input f, additional values can be appended or removed from the +% calculated gradient after the inverse transform is calculated by +% setting the optional input pad to true (the default). These values +% are known from the symmetry of the output sequence. % -% For additional details, see [1]. +% For additional details on gradient calculation using DTTs, see [1]. % % [1] E. Wise, J. Jaros, B. Cox, and B. Treeby, "Pseudospectral % time-domain (PSTD) methods for the wave equation: Realising boundary % conditions with discrete sine and cosine transforms", 2020. % % INPUTS: -% f - matrix or vector to find the gradient of -% dx - grid point spacing -% dtt_type - type of discrete trigonometric transform, this should -% correspond to the assumed input symmetry of the input -% function, where: +% f - Vector to find the gradient of. +% dx - Grid point spacing. +% dtt_type - Type of discrete trigonometric transform. This should +% correspond to the assumed input symmetry of the input +% function, where: % -% 1: DCT-I WSWS -% 2: DCT-II HSHS -% 3: DCT-III WSWA -% 4: DCT-IV HSHA -% 5: DST-I WAWA -% 6: DST-II HAHA -% 7: DST-III WAWS -% 8: DST-IV HAHS +% 1: DCT-I WSWS +% 2: DCT-II HSHS +% 3: DCT-III WSWA +% 4: DCT-IV HSHA +% 5: DST-I WAWA +% 6: DST-II HAHA +% 7: DST-III WAWS +% 8: DST-IV HAHS % % OPTIONAL INPUTS: -% shift - integer controlling whether derivative is shifted to a -% staggered grid (default = 0), where +% shift - Integer controlling whether derivative is shifted to a +% staggered grid (default = 0), where % -% 0: no shift -% 1: shift by + dx/2 -% 2: shift by - dx/2 +% 0: no shift +% 1: shift by + dx/2 +% 2: shift by - dx/2 +% +% align_output - Boolean controlling whether the returned values are +% padded and trimmed based on the implied symmetry so +% the output is the same length as the input (default = +% true). Note, if align_output is false, then the +% gradient calculated for shift = 1 and shift = 2 will +% be the same. % % OUTPUTS: -% dfdx - gradient of the input function +% dfdx - Gradient of the input function. % % ABOUT: -% author - Bradley Treeby -% date - 23 April 2013 -% last update - 2 April 2020 +% author - Bradley Treeby +% date - 23 April 2013 +% last update - 20 April 2020 % % Copyright (C) 2013-2020 Bradley Treeby % @@ -84,15 +92,22 @@ DST4 = 8; % HAHS % check for shift input -if nargin < 4 - shift = false; +if (nargin < 4) || isempty(shift) + shift = 0; end -% make sure the input is 1D -if numDim(f) > 1 - error('Input function must be 1D.'); +% check for pad input +if (nargin < 5) || isempty(align_output) + align_output = true; end +% check inputs +validateattributes(f, {'numeric'}, {'real', 'vector', 'nonsparse'}, 'gradientDTT1D', 'f', 1); +validateattributes(dx, {'numeric'}, {'real', 'scalar', 'nonsparse'}, 'gradientDTT1D', 'dx', 2); +validateattributes(dtt_type, {'numeric'}, {'integer', 'scalar', 'nonsparse', '>=', 1, '<=', 8}, 'gradientDTT1D', 'dtt_type', 3); +validateattributes(shift, {'numeric'}, {'integer', 'scalar', 'nonsparse', '>=', 0, '<=', 2}, 'gradientDTT1D', 'shift', 4); +validateattributes(align_output, {'logical'}, {'scalar'}, 'gradientDTT1D', 'align_output', 5); + % reshape to be a row vector f = reshape(f, 1, []); @@ -341,161 +356,163 @@ % add back in the implied values so the output is the same length as the % input -switch dtt_type - case 1 - - switch shift - case 0 - - % WSWS -> WAWA, add both endpoints - deriv = [0, deriv, 0]; - - case 1 - - % WSWS -> HAHA, shift right, mirror right endpoint - deriv = [deriv, -deriv(end)]; - - case 2 - - % WSWS -> HAHA, shift left, mirror left endpoint - deriv = [-deriv(1), deriv]; - - end - - case 2 - - switch shift - case 0 - - % HSHS -> HAHA, no change - - case 1 - - % HSHS -> WAWA, shift right, append zero - deriv = [deriv, 0]; - - case 2 - - % HSHS -> WAWA, shift left, prepend zero - deriv = [0, deriv]; - - end - - case 3 - - switch shift - case 0 - - % WSWA -> WAWS, prepend zero, remove right endpoint - deriv = [0, deriv(1:end - 1)]; - - case 1 - - % WSWA -> HAHS, shift right, no change - - case 2 - - % WSWA -> HAHS, shift left, mirror left endpoint, remove - % right endpoint - deriv = [-deriv(1), deriv(1:end - 1)]; - - end - - case 4 - - switch shift - case 0 - - % HSHA -> HAHS, no change - - case 1 - - % HSHA -> WAWS, shift right, no change - - case 2 - - % HSHA -> WAWS, shift left, prepend zero, remove right - % endpoint - deriv = [0, deriv(1:end - 1)]; - - end - - case 5 - - switch shift - case 0 - - % WAWA -> WSWS, remove both endpoints - deriv = deriv(2:end - 1); - - case 1 - - % WAWA -> HSHS, shift right, remove left endpoint - deriv = deriv(2:end); - - case 2 - - % WAWA -> HSHS, shift left, remove right endpoint - deriv = deriv(1:end - 1); - - end - - case 6 - - switch shift - case 0 - - % HAHA -> HSHS, no change - - case 1 - - % HAHA -> WSWS, shift right, remove left endpoint - deriv = deriv(2:end); - - case 2 - - % HAHA -> WSWS, shift left, remove right endpoint - deriv = deriv(1:end - 1); - - end - - case 7 - - switch shift - case 0 - - % WAWS -> WSWA, remove left endpoint, append zero - deriv = [deriv(2:end), 0]; - - case 1 - - % WAWS -> HSHA, shift right, remove left endpoint, mirror - % right endpoint - deriv = [deriv(2:end), -deriv(end)]; - - case 2 - - % WAWS -> HSHA, shift left, no change - - end - - case 8 - - switch shift - case 0 - - % HAHS -> HSHA, no change - - case 1 - - % HAHS -> WSWA, shift right, remove left endpoint, append - % zero - deriv = [deriv(2:end), 0]; - - case 2 - - % HAHS -> WSWA, shift left, no change - - end - +if align_output + switch dtt_type + case 1 + + switch shift + case 0 + + % WSWS -> WAWA, add both endpoints + deriv = [0, deriv, 0]; + + case 1 + + % WSWS -> HAHA, shift right, mirror right endpoint + deriv = [deriv, -deriv(end)]; + + case 2 + + % WSWS -> HAHA, shift left, mirror left endpoint + deriv = [-deriv(1), deriv]; + + end + + case 2 + + switch shift + case 0 + + % HSHS -> HAHA, no change + + case 1 + + % HSHS -> WAWA, shift right, append zero + deriv = [deriv, 0]; + + case 2 + + % HSHS -> WAWA, shift left, prepend zero + deriv = [0, deriv]; + + end + + case 3 + + switch shift + case 0 + + % WSWA -> WAWS, prepend zero, remove right endpoint + deriv = [0, deriv(1:end - 1)]; + + case 1 + + % WSWA -> HAHS, shift right, no change + + case 2 + + % WSWA -> HAHS, shift left, mirror left endpoint, remove + % right endpoint + deriv = [-deriv(1), deriv(1:end - 1)]; + + end + + case 4 + + switch shift + case 0 + + % HSHA -> HAHS, no change + + case 1 + + % HSHA -> WAWS, shift right, no change + + case 2 + + % HSHA -> WAWS, shift left, prepend zero, remove right + % endpoint + deriv = [0, deriv(1:end - 1)]; + + end + + case 5 + + switch shift + case 0 + + % WAWA -> WSWS, remove both endpoints + deriv = deriv(2:end - 1); + + case 1 + + % WAWA -> HSHS, shift right, remove left endpoint + deriv = deriv(2:end); + + case 2 + + % WAWA -> HSHS, shift left, remove right endpoint + deriv = deriv(1:end - 1); + + end + + case 6 + + switch shift + case 0 + + % HAHA -> HSHS, no change + + case 1 + + % HAHA -> WSWS, shift right, remove left endpoint + deriv = deriv(2:end); + + case 2 + + % HAHA -> WSWS, shift left, remove right endpoint + deriv = deriv(1:end - 1); + + end + + case 7 + + switch shift + case 0 + + % WAWS -> WSWA, remove left endpoint, append zero + deriv = [deriv(2:end), 0]; + + case 1 + + % WAWS -> HSHA, shift right, remove left endpoint, mirror + % right endpoint + deriv = [deriv(2:end), -deriv(end)]; + + case 2 + + % WAWS -> HSHA, shift left, no change + + end + + case 8 + + switch shift + case 0 + + % HAHS -> HSHA, no change + + case 1 + + % HAHS -> WSWA, shift right, remove left endpoint, append + % zero + deriv = [deriv(2:end), 0]; + + case 2 + + % HAHS -> WSWA, shift left, no change + + end + + end end \ No newline at end of file