Version 1.1 (see change log)

ucl-bug · Apr 21, 2020 · 3c7fa1e · 3c7fa1e
1 parent 35ca9d6
commit 3c7fa1e
Show file tree

Hide file tree

Showing 9 changed files with 1,132 additions and 228 deletions.
diff --git 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: [email protected]
@@ -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 <https://www.gnu.org/licenses/>.
+You should have received a copy of the GNU General Public License along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+## 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
diff --git 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
@@ -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 → examples/example_wave_eq_pstd_1D.m b/examples/example_wave_eq_pstd.m → 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