updated readme file and cleaned up code
20 changed files+450−326
.DS_Storemodified+0−0View file
Binary file not shown.
README.mdmodified+14−1View file
@@ -3,4 +3,17 @@
33
44 ## Dependencies
55
6-This project requires [finufft](https://github.com/flatironinstitute/finufft.git).
6+This project requires [Matlab](https://www.mathworks.com/help/install/ug/install-products-with-internet-connection.html)
7+and [finufft](https://github.com/flatironinstitute/finufft.git).
8+
9+## Example runs
10+Enter experiment number (expNum) in the driver file main.m to run built-in examples.
11+Parameters for example runs are found in get_experimentParameters(expNum):
12+expNum = 1, interpolation order 3 manufactured solution convergence study with 10 randomly located sources in [-2,2] up to time 3*pi
13+expNum = 2, interpolation order 8 scattering from 150 springs located randomly on [-2,2] up to time 10*pi
14+expNum = 3, interpolation order 8 scattering from 1000 springs located randomly on [-2,2] up to time 3*pi with narrow incident wave
15+expNum = 4, interpolation order 8 scattering from 10000 springs located randomly on [-2,2] up to time 3*pi with narrow incident wave
16+expNum = 5, interpolation order 8 scattering from 2 stiff springs located on -0.5 and 0.5 up to time 30*pi with wide incident wave to simulate Fabry-Perot interferometer
17+expNum = 6, interpolation order 8 scattering from 10 stiff springs located uniformly between -2 to 2 up to time 40*pi final time with wide incident wave
18+expNum = 7, interpolation order 8 scattering from 200 springs located uniformly between -2 to 2 up to time 40*pi final time with wide incident wave
19+expNum = 999, experiment number reserved for running tests and debugging
wfp_1Dspringscattering/.DS_Storemodified+0−0View file
Binary file not shown.
wfp_1Dspringscattering/checkMatchingGrids.mmodified+1−0View file
@@ -8,6 +8,7 @@ end
88
99 tol = 1e-12;
1010 if(tGridErr>tol || xGridErr>tol)
11+ fprintf('tGridErr = %1.2e, xGridErr = %1.2e\n',tGridErr,xGridErr);
1112 warning('Time or space grids do not match to report self-convergence');
1213 end
1314
wfp_1Dspringscattering/get_experimentParameters.asvadded+310−0View file
@@ -0,0 +1,310 @@
1+function [plotOpt,saveWorkspaceOpt,addErrResults,savePlot,...
2+ logScale,printTimeStep,evalSol,order,numResolutions,solnType,tFinal,...
3+ Nx,Nt_sol,M,maxNumNeighbors,src_dmn,ds,...
4+ uniform_sgrid,betaMax,dataParam,dt0_min,uniform_beta,keepdt]...
5+ =get_experimentParameters(expNum)
6+
7+dataParam = struct;
8+
9+switch expNum
10+ case 1
11+ % testing buttons
12+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
13+ saveWorkspaceOpt = 0;
14+ addErrResults = 0;
15+ savePlot = 0;
16+ logScale = 0; % '1' to use log scale in heat maps
17+ printTimeStep = 0; % '1' to print out each time step
18+ evalSol = 1;
19+
20+ % parameters
21+ order = 2; % order of accuracy = interpolation
22+ numResolutions = 1; % num of grid resolutions
23+ solnType = 'ms'; % type of solution 'ms' or 'true'
24+
25+ tFinal = 3*pi; % final time
26+
27+ %%% Choose number of spatial and temporal set for solution computation
28+ Nx = 5; % size of spatial grid for testing and plotting
29+ Nt_sol = 5; % size of the time grid for testing and plotting
30+
31+ %%% sources
32+ M = 10; % number of sources
33+ maxNumNeighbors = NaN; % set max number of neighbors
34+ src_dmn = [-2,2];
35+ ds = 1e-4; % set min distance between sources
36+ uniform_sgrid = 1;
37+
38+ %%% spring constants
39+ betaMax = 3; uniform_beta = 1;
40+
41+ %%% incident pulse parameter
42+ keepdt = 0;
43+
44+ dt0_min = 0.01; % RBC constraint
45+
46+ case 2 % localization plots
47+ % testing buttons
48+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
49+ saveWorkspaceOpt = 1;
50+ addErrResults = 0;
51+ savePlot = 0;
52+ logScale = 0; % '1' to use log scale in heat maps
53+ printTimeStep = 1; % '1' to print out each time step
54+ evalSol = 1;
55+
56+ % parameters
57+ order = 8; % order of accuracy = interpolation
58+ numResolutions = 2; % num of grid resolutions
59+ solnType = 'true'; % type of solution 'ms' or 'true'
60+
61+ tFinal = 10*pi; % final time
62+
63+ %%% Choose number of spatial and temporal set for solution computation
64+ Nx = 400; % size of spatial grid for testing and plotting
65+ Nt_sol = 400; % size of the time grid for testing and plotting
66+
67+ %%% sources
68+ M = 150; % number of sources
69+ maxNumNeighbors = 10; % set max number of neighbors
70+ src_dmn = [-2,2];
71+ ds = 1e-4; % set min distance between sources
72+ uniform_sgrid = 0;
73+
74+ %%% spring constants
75+ betaMax = 10; uniform_beta = 0;
76+
77+ %%% incident pulse parameter
78+ dataParam.mu = 30; dataParam.t0 = -3;
79+ keepdt = 0;
80+
81+ dt0_min = 0.1;
82+
83+ case 3 % M=1000, mu = M^2
84+
85+ % testing buttons
86+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
87+ saveWorkspaceOpt = 1;
88+ addErrResults = 0;
89+ savePlot = 0;
90+ logScale = 0; % '1' to use log scale in heat maps
91+ printTimeStep = 1; % '1' to print out each time step
92+ evalSol = 1;
93+
94+ % parameters
95+ order = 8; % order of accuracy = interpolation
96+ numResolutions = 2; % num of grid resolutions
97+ solnType = 'true'; % type of solution 'ms' or 'true'
98+
99+ tFinal = 3*pi; % final time
100+
101+ %%% Choose number of spatial and temporal set for solution computation
102+ Nx = 800; % size of spatial grid for testing and plotting
103+ Nt_sol = 800; % size of the time grid for testing and plotting
104+
105+ %%% sources
106+ M = 1000; % number of sources
107+ maxNumNeighbors = 50; % set max number of neighbors
108+ src_dmn = [-2,2];
109+ ds = 1e-4; % set min distance between sources
110+ uniform_sgrid = 0;
111+
112+ %%% spring constants
113+ betaMax = 3; uniform_beta = 0;
114+
115+ %%% incident pulse parameter
116+ dataParam.mu = M^2; dataParam.t0 = -3;
117+ keepdt = 0;
118+
119+ dt0_min = 0.1;
120+ case 4 % M=10000, mu = M^2
121+ % testing buttons
122+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
123+ saveWorkspaceOpt = 1;
124+ addErrResults = 0;
125+ savePlot = 0;
126+ logScale = 0; % '1' to use log scale in heat maps
127+ printTimeStep = 1; % '1' to print out each time step
128+ evalSol = 1;
129+
130+ % parameters
131+ order = 8; % order of accuracy = interpolation
132+ numResolutions = 2; % num of grid resolutions
133+ solnType = 'true'; % type of solution 'ms' or 'true'
134+
135+ tFinal = 3*pi; % final time
136+
137+ %%% Choose number of spatial and temporal set for solution computation
138+ Nx = 400; % size of spatial grid for testing and plotting
139+ Nt_sol = 800; % size of the time grid for testing and plotting
140+
141+ %%% sources
142+ M = 10000; % number of sources
143+ maxNumNeighbors = 500; % set max number of neighbors
144+ src_dmn = [-2,2];
145+ ds = 1e-4; % set min distance between sources
146+ uniform_sgrid = 0;
147+
148+ %%% spring constants
149+ betaMax = 3; uniform_beta = 0;
150+
151+ %%% incident pulse parameter
152+ dataParam.mu = M^2; dataParam.t0 = -3;
153+ dataParam.doubleTimeStep = 1;
154+ keepdt = 0;
155+
156+ dt0_min = 0.1;
157+ case 5 % FP
158+ % testing buttons
159+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
160+ saveWorkspaceOpt = 1;
161+ addErrResults = 0;
162+ savePlot = 0;
163+ logScale = 0; % '1' to use log scale in heat maps
164+ printTimeStep = 0; % '1' to print out each time step
165+ evalSol = 1;
166+
167+ % parameters
168+ order = 8; % order of accuracy = interpolation
169+ numResolutions = 2; % num of grid resolutions
170+ solnType = 'true'; % type of solution 'ms' or 'true'
171+
172+ tFinal = 30*pi; % final time
173+
174+ %%% Choose number of spatial and temporal set for solution computation
175+ Nx = 400; % size of spatial grid for testing and plotting
176+ Nt_sol = 500; % size of the time grid for testing and plotting
177+
178+ %%% sources
179+ M = 2; % number of sources
180+ maxNumNeighbors = 2; % set max number of neighbors
181+ src_dmn = [-0.5,0.5];
182+ ds = 1e-4; % set min distance between sources
183+ uniform_sgrid = 1;
184+
185+ %%% spring constants
186+ betaMax = 100; uniform_beta = 1;
187+
188+ %%% incident pulse parameter
189+ dataParam.mu = 5; dataParam.t0 = -3;
190+ keepdt = 0;
191+
192+ dt0_min = 0.02; % RBC constraint
193+ dt0_min = min(0.01,dt0_min); % stability constraint
194+
195+ case 6 % FP M = 10
196+ % testing buttons
197+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
198+ saveWorkspaceOpt = 1;
199+ addErrResults = 0;
200+ savePlot = 0;
201+ logScale = 0; % '1' to use log scale in heat maps
202+ printTimeStep = 1; % '1' to print out each time step
203+ evalSol = 1;
204+
205+ % parameters
206+ order = 8; % order of accuracy = interpolation
207+ numResolutions = 2; % num of grid resolutions
208+ solnType = 'true'; % type of solution 'ms' or 'true'
209+
210+ tFinal = 40*pi; % final time
211+
212+ %%% Choose number of spatial and temporal set for solution computation
213+ Nx = 800; % size of spatial grid for testing and plotting
214+ Nt_sol = 800; % size of the time grid for testing and plotting
215+
216+ %%% sources
217+ M = 10; % number of sources
218+ maxNumNeighbors = 10; % set max number of neighbors
219+ src_dmn = [-2,2];
220+ ds = 1e-4; % set min distance between sources
221+ uniform_sgrid = 1;
222+
223+ %%% spring constants
224+ betaMax = 100; uniform_beta = 1;
225+
226+ %%% incident pulse parameter
227+ dataParam.mu = 5; dataParam.t0 = -4;
228+ keepdt = 0;
229+
230+ dt0_min = 0.02; % RBC constraint
231+ dt0_min = min(0.005,dt0_min); % stability constraint
232+ case 7 % FP M = 200
233+ % testing buttons
234+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
235+ saveWorkspaceOpt = 1;
236+ addErrResults = 0;
237+ savePlot = 0;
238+ logScale = 0; % '1' to use log scale in heat maps
239+ printTimeStep = 0; % '1' to print out each time step
240+ evalSol = 1;
241+
242+ % parameters
243+ order = 8; % order of accuracy = interpolation
244+ numResolutions = 2; % num of grid resolutions
245+ solnType = 'true'; % type of solution 'ms' or 'true'
246+
247+ tFinal = 40*pi; % final time
248+
249+ %%% Choose number of spatial and temporal set for solution computation
250+ Nx = 800; % size of spatial grid for testing and plotting
251+ Nt_sol = 800; % size of the time grid for testing and plotting
252+
253+ %%% sources
254+ M = 200; % number of sources
255+ maxNumNeighbors = 20; % set max number of neighbors
256+ src_dmn = [-2,2];
257+ ds = 1e-4; % set min distance between sources
258+ uniform_sgrid = 1;
259+
260+ %%% spring constants
261+ betaMax = 1; uniform_beta = 1;
262+
263+ %%% incident pulse parameter
264+ dataParam.mu = 5; dataParam.t0 = -4;
265+ keepdt = 0;
266+
267+ dt0_min = 0.02; % RBC constraint
268+ dt0_min = min(0.005,dt0_min); % stability constraint
269+ case 999
270+ % testing buttons
271+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
272+ saveWorkspaceOpt = 0;
273+ addErrResults = 0;
274+ savePlot = 0;
275+ logScale = 0; % '1' to use log scale in heat maps
276+ printTimeStep = 0; % '1' to print out each time step
277+ evalSol = 1;
278+
279+ % parameters
280+ order = 2; % order of accuracy = interpolation
281+ numResolutions = 3; % num of grid resolutions
282+ solnType = 'true'; % type of solution 'ms' or 'true'
283+
284+ tFinal = 3*pi; % final time
285+
286+ %%% Choose number of spatial and temporal set for solution computation
287+ Nx = 5; % size of spatial grid for testing and plotting
288+ Nt_sol = 5; % size of the time grid for testing and plotting
289+
290+ %%% sources
291+ M = 10; % number of sources
292+ maxNumNeighbors = NaN; % set max number of neighbors
293+ src_dmn = [-2,2];
294+ ds = 1e-4; % set min distance between sources
295+ uniform_sgrid = 1;
296+
297+ %%% spring constants
298+ betaMax = 3; uniform_beta = 1;
299+
300+ %%% incident pulse parameter
301+ dataParam.mu = 30; dataParam.t0 = -3;
302+ keepdt = 0;
303+
304+ dt0_min = 0.01; % RBC constraint
305+ otherwise
306+ disp('experiment number inapplicable');
307+end
308+
309+end
310+
wfp_1Dspringscattering/get_experimentParameters.mmodified+69−32View file
@@ -1,12 +1,13 @@
11 function [plotOpt,saveWorkspaceOpt,addErrResults,savePlot,...
22 logScale,printTimeStep,evalSol,order,numResolutions,solnType,tFinal,...
33 Nx,Nt_sol,M,maxNumNeighbors,src_dmn,ds,...
4- uniform_sgrid,betaMax,mu,dt0_min,uniform_beta,keepdt]...
4+ uniform_sgrid,betaMax,dataParam,dt0_min,uniform_beta,keepdt]...
55 =get_experimentParameters(expNum)
66
7+dataParam = struct;
8+
79 switch expNum
810 case 1
9-
1011 % testing buttons
1112 plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
1213 saveWorkspaceOpt = 0;
@@ -17,33 +18,32 @@ switch expNum
1718 evalSol = 1;
1819
1920 % parameters
20- order = 2; % order of accuracy = interpolation
21+ order = 3; % order of accuracy = interpolation
2122 numResolutions = 3; % num of grid resolutions
2223 solnType = 'ms'; % type of solution 'ms' or 'true'
2324
2425 tFinal = 3*pi; % final time
2526
2627 %%% Choose number of spatial and temporal set for solution computation
27- Nx = 10; % size of spatial grid for testing and plotting
28- Nt_sol = 10; % size of the time grid for testing and plotting
28+ Nx = 5; % size of spatial grid for testing and plotting
29+ Nt_sol = 5; % size of the time grid for testing and plotting
2930
3031 %%% sources
31- M = 2; % number of sources
32+ M = 10; % number of sources
3233 maxNumNeighbors = NaN; % set max number of neighbors
33- src_dmn = [-1,1];
34+ src_dmn = [-2,2];
3435 ds = 1e-4; % set min distance between sources
35- uniform_sgrid = 1;
36+ uniform_sgrid = 0;
3637
3738 %%% spring constants
38- betaMax = 3; uniform_beta = 1;
39+ betaMax = 3; uniform_beta = 0;
3940
4041 %%% incident pulse parameter
41- mu = 50;
42- keepdt = 1;
42+ keepdt = 0;
4343
44- dt0_min = 0.1;
45- case 2
44+ dt0_min = 0.01; % RBC constraint
4645
46+ case 2 % localization plots
4747 % testing buttons
4848 plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
4949 saveWorkspaceOpt = 1;
@@ -75,12 +75,12 @@ switch expNum
7575 betaMax = 10; uniform_beta = 0;
7676
7777 %%% incident pulse parameter
78- mu = 30;
78+ dataParam.mu = 30; dataParam.t0 = -3;
7979 keepdt = 0;
8080
8181 dt0_min = 0.1;
8282
83- case 3
83+ case 3 % M=1000, mu = M^2
8484
8585 % testing buttons
8686 plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
@@ -113,18 +113,18 @@ switch expNum
113113 betaMax = 3; uniform_beta = 0;
114114
115115 %%% incident pulse parameter
116- mu = M^2;
116+ dataParam.mu = M^2; dataParam.t0 = -3;
117117 keepdt = 0;
118118
119119 dt0_min = 0.1;
120- case 4
120+ case 4 % M=10000, mu = M^2
121121 % testing buttons
122122 plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
123123 saveWorkspaceOpt = 1;
124124 addErrResults = 0;
125125 savePlot = 0;
126126 logScale = 0; % '1' to use log scale in heat maps
127- printTimeStep = 0; % '1' to print out each time step
127+ printTimeStep = 1; % '1' to print out each time step
128128 evalSol = 1;
129129
130130 % parameters
@@ -135,12 +135,12 @@ switch expNum
135135 tFinal = 3*pi; % final time
136136
137137 %%% Choose number of spatial and temporal set for solution computation
138- Nx = 800; % size of spatial grid for testing and plotting
138+ Nx = 400; % size of spatial grid for testing and plotting
139139 Nt_sol = 800; % size of the time grid for testing and plotting
140140
141141 %%% sources
142142 M = 10000; % number of sources
143- maxNumNeighbors = 50; % set max number of neighbors
143+ maxNumNeighbors = 500; % set max number of neighbors
144144 src_dmn = [-2,2];
145145 ds = 1e-4; % set min distance between sources
146146 uniform_sgrid = 0;
@@ -149,7 +149,8 @@ switch expNum
149149 betaMax = 3; uniform_beta = 0;
150150
151151 %%% incident pulse parameter
152- mu = M^2;
152+ dataParam.mu = M^2; dataParam.t0 = -3;
153+ dataParam.doubleTimeStep = 1;
153154 keepdt = 0;
154155
155156 dt0_min = 0.1;
@@ -185,7 +186,7 @@ switch expNum
185186 betaMax = 100; uniform_beta = 1;
186187
187188 %%% incident pulse parameter
188- mu = 5;
189+ dataParam.mu = 5; dataParam.t0 = -3;
189190 keepdt = 0;
190191
191192 dt0_min = 0.02; % RBC constraint
@@ -194,7 +195,7 @@ switch expNum
194195 case 6 % FP M = 10
195196 % testing buttons
196197 plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
197- saveWorkspaceOpt = 0;
198+ saveWorkspaceOpt = 1;
198199 addErrResults = 0;
199200 savePlot = 0;
200201 logScale = 0; % '1' to use log scale in heat maps
@@ -203,17 +204,17 @@ switch expNum
203204
204205 % parameters
205206 order = 8; % order of accuracy = interpolation
206- numResolutions = 3; % num of grid resolutions
207+ numResolutions = 2; % num of grid resolutions
207208 solnType = 'true'; % type of solution 'ms' or 'true'
208209
209- tFinal = 3*pi; % final time
210+ tFinal = 40*pi; % final time
210211
211212 %%% Choose number of spatial and temporal set for solution computation
212- Nx = 10; % size of spatial grid for testing and plotting
213- Nt_sol = 10; % size of the time grid for testing and plotting
213+ Nx = 800; % size of spatial grid for testing and plotting
214+ Nt_sol = 800; % size of the time grid for testing and plotting
214215
215216 %%% sources
216- M = 2; % number of sources
217+ M = 10; % number of sources
217218 maxNumNeighbors = 10; % set max number of neighbors
218219 src_dmn = [-2,2];
219220 ds = 1e-4; % set min distance between sources
@@ -223,11 +224,11 @@ switch expNum
223224 betaMax = 100; uniform_beta = 1;
224225
225226 %%% incident pulse parameter
226- mu = 5;
227+ dataParam.mu = 5; dataParam.t0 = -4;
227228 keepdt = 0;
228229
229230 dt0_min = 0.02; % RBC constraint
230- dt0_min = min(0.001,dt0_min); % stability constraint
231+ dt0_min = min(0.005,dt0_min); % stability constraint
231232 case 7 % FP M = 200
232233 % testing buttons
233234 plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
@@ -260,11 +261,47 @@ switch expNum
260261 betaMax = 1; uniform_beta = 1;
261262
262263 %%% incident pulse parameter
263- mu = 5;
264+ dataParam.mu = 5; dataParam.t0 = -4;
264265 keepdt = 0;
265266
266267 dt0_min = 0.02; % RBC constraint
267- dt0_min = min(0.01,dt0_min); % stability constraint
268+ dt0_min = min(0.005,dt0_min); % stability constraint
269+ case 999
270+ % testing buttons
271+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
272+ saveWorkspaceOpt = 0;
273+ addErrResults = 0;
274+ savePlot = 0;
275+ logScale = 0; % '1' to use log scale in heat maps
276+ printTimeStep = 0; % '1' to print out each time step
277+ evalSol = 1;
278+
279+ % parameters
280+ order = 2; % order of accuracy = interpolation
281+ numResolutions = 3; % num of grid resolutions
282+ solnType = 'true'; % type of solution 'ms' or 'true'
283+
284+ tFinal = 3*pi; % final time
285+
286+ %%% Choose number of spatial and temporal set for solution computation
287+ Nx = 5; % size of spatial grid for testing and plotting
288+ Nt_sol = 5; % size of the time grid for testing and plotting
289+
290+ %%% sources
291+ M = 10; % number of sources
292+ maxNumNeighbors = NaN; % set max number of neighbors
293+ src_dmn = [-2,2];
294+ ds = 1e-4; % set min distance between sources
295+ uniform_sgrid = 1;
296+
297+ %%% spring constants
298+ betaMax = 3; uniform_beta = 1;
299+
300+ %%% incident pulse parameter
301+ dataParam.mu = 30; dataParam.t0 = -3;
302+ keepdt = 0;
303+
304+ dt0_min = 0.01; % RBC constraint
268305 otherwise
269306 disp('experiment number inapplicable');
270307 end
wfp_1Dspringscattering/get_gvalue.mdeleted+0−15View file
@@ -1,15 +0,0 @@
1-function g = get_gvalue(t,s,j,sig,beta,M,mu,t0,addProj)
2-% GET_GVALUE get the data at a given time t
3-%
4-% g = get_gvalue(t,s,j,sig,beta,M,mu,t0) returns the data g for each spring
5-% scatterer condition at a time 't'. 's' is a vector containing the source
6-% locations s = [s(1),...,s(M)], where 'M' is the number of sources. 'sig'
7-% is a cell that holds manufactured density values (Gaussians) where
8-% sig{j} = @(t) exp(mu(j)*(t - t0(j))). 'beta' is a factor of the spring
9-% constant, and 'j' represents the index of the current spring scatterer
10-
11-xval = s(j); % location at the current source
12-U = get_ExactSol(xval,t,mu,t0,M,s,addProj); % get the manufactured solution at xval
13-g = -sig(t,j) - beta*U;
14-
15-end
wfp_1Dspringscattering/get_gvalueFromIncomingField.mdeleted+0−18View file
@@ -1,18 +0,0 @@
1-function g = get_gvalueFromIncomingField(t,s,j,beta,mu,t0)
2-% GET_GVALUEFROMINCOMINGFIELD get data from incoming field
3-%
4-% g = get_gvalueFromIncomingField(t,s,j,beta) returns the data g for each spring
5-% scatterer condition at a time 't'. 's' is a vector containing the source
6-% locations s = [s(1),...,s(M)], where 'M' is the number of sources.
7-% 'beta' is a factor of the spring constant, and 'j' represents the index
8-% of the current spring scatterer. 'mu' and 't0' are parameters needed for
9-% the incident wave.
10-
11-x = s(j);
12-uin = get_uIncident(x,t,mu,t0);
13-g = beta*uin;
14-
15-% testing
16-% g = beta*cos(10*(s(j) - t - t0));
17-
18-end
wfp_1Dspringscattering/get_msOrIncidentField.mdeleted+0−40View file
@@ -1,40 +0,0 @@
1-function get_msOrIncidentField(M,tFinal,mu,solnType)
2-
3-
4-if(strcmp(solnType,'ms'))
5- override = 0; % turn on '1' to avoid restrictions on chosen test function
6-
7-mu = linspace(40,50,M)';
8-t0 = linspace(1,3,M)';
9-
10-% choose the initial time step
11-dt0_vec = pi./(4*sqrt(mu*log(1/eps)));
12-dt0 = min(dt0_vec);
13-
14-% ensure that the forcing function is bandlimited within the computational time domain
15-if(any(2*t0>tFinal)&&(override == 0))
16- error('Forcing inappropriate for testing. For tFinal=%1.2f\nChoose other mu or tFinal',tFinal);
17-end
18-
19-% Define the density functions to be Gaussians
20-sig = @(t) exp(-mu.*((t - t0).^2));
21-
22-dataParam.sig = sig;
23-dataParam.t0 = t0;
24-dataParam.mu = mu;
25-
26-else
27-
28-t0 = -3;
29-
30-% choose the initial time step
31-dt0 = pi/(4*sqrt(mu*log(1/eps)));
32-
33-dataParam.mu = mu;
34-dataParam.t0 = t0;
35-end
36-
37-
38-
39-
40-end
wfp_1Dspringscattering/get_uIncidentInfo.mmodified+29−4View file
@@ -1,4 +1,4 @@
1-function [dt0,dataParam] = get_uIncidentInfo(mu)
1+function [dt0,dataParam] = get_uIncidentInfo(dataParam)
22 % GET_UINCIDENTINFO prepares information needed for an incident wave of the
33 % form uin = exp(-mu*(x - t - t0).^2);
44 %
@@ -14,12 +14,37 @@ function [dt0,dataParam] = get_uIncidentInfo(mu)
1414 % dataParam: is struct to hold mu and t0 (for ease of switching between
1515 % types of test solutions
1616
17-t0 = -3;
17+if nargin == 0, test_get_uIncidentInfo(); return; end
18+
19+mu = dataParam.mu;
20+t0 = dataParam.t0;
21+
22+if(isnan(t0))
23+ t0 = -sqrt(log(1/eps)./mu);
24+end
1825
1926 % choose the initial time step
2027 dt0 = pi/(4*sqrt(mu*log(1/eps)));
2128
22-dataParam.mu = mu;
23-dataParam.t0 = t0;
29+if(isfield(dataParam,'doubleTimeStep'))
30+ if(dataParam.doubleTimeStep == 1)
31+ dt0 = 2*dt0;
32+ end
33+end
34+
35+dataParam.t0 = t0;
36+
37+end
38+
39+function test_get_uIncidentInfo()
40+
41+mu = 5;
42+[dt0,dataParam] = get_uIncidentInfo(mu);
43+
44+x = linspace(-pi,pi,1000); t = linspace(0,3*pi,1000); t0 = dataParam.t0;
45+uin = get_uIncident(x,t,mu,t0);
46+
47+figure(1)
48+imagesc(x,t,uin);
2449
2550 end
wfp_1Dspringscattering/get_xgrid.mmodified+1−1View file
@@ -16,7 +16,7 @@ h_temp = h0;
1616
1717 xmin = (ax + 2*h_temp*W); xmax = (bx - 2*h_temp*W);
1818 x = linspace(xmin,xmax,Nx)';
19-if(xmin>as || xmax<bs)
19+if(xmin>=as || xmax<=bs)
2020 xmin,xmax,h0
2121 error('Sources are outside allowed region; reduce the time-step or narrow down the domain where sources live');
2222 end
wfp_1Dspringscattering/main.mmodified+6−6View file
@@ -2,16 +2,16 @@ startMatlabFile
22
33 %% Set file directories to save data, figures, tables and movies
44 addpath ./utils;
5-dataFile = '/Users/nalhassanieh/Desktop/WFP_results/data';
5+dataFile = './data';
66 figFile = './fig';
77
88 %% Problem parameters
9-expNum = 6; % experiment number
9+expNum = 1; % experiment number
1010
1111 [plotOpt,saveWorkspaceOpt,addErrResults,savePlot,...
1212 logScale,printTimeStep,evalSol,order,numResolutions,solnType,tFinal,...
1313 Nx,Nt_sol,M,maxNumNeighbors,src_dmn,ds,...
14- uniform_sgrid,betaMax,mu,dt0_min,uniform_beta,keepdt]...
14+ uniform_sgrid,betaMax,dataParam,dt0_min,uniform_beta,keepdt]...
1515 = get_experimentParameters(expNum);
1616
1717 %%% domain
@@ -31,7 +31,7 @@ tol = 1e-12; % error tolerance
3131
3232 %% Manufactured density values for testing or incident wave for true scattering
3333 if(strcmp(solnType,'ms'));[dt0,dataParam] = manufacturedSolution(M,tFinal);
34-else; [dt0,dataParam] = get_uIncidentInfo(mu); end
34+else; [dt0,dataParam] = get_uIncidentInfo(dataParam); end
3535 dt0 = min(dt0_min,dt0); % domain RBC constraint
3636
3737 % fix initial time step
@@ -39,7 +39,7 @@ dt0 = min(dt0_min,dt0); % domain RBC constraint
3939 Nt0 = 2*ceil(0.5*tFinal/dt0);
4040 dt0 = tFinal/Nt0;
4141
42-%% Compute some values for the window function and GL
42+%% Compute some values for the window function and GL*
4343 P = W; % GL nodes
4444 gl = glwt_prep(P); % GL nodes and weights on [-1,1]
4545
@@ -163,7 +163,7 @@ for m = 1:numResolutions
163163 % update values for the next time step
164164 [an,bn,sn,sn_hat] = prepForNextStep(anp1,bnp1,sn,sn_hat,snIdxNow,tgh,s,N,tol);
165165
166- if(printTimeStep == 1)
166+ if(printTimeStep == 1 && mod(totalSteps,100) == 0)
167167 fprintf('t = %1.4e\n',t);
168168 end
169169
wfp_1Dspringscattering/manufacturedSolution.mmodified+19−2View file
@@ -22,10 +22,15 @@ function [dt0,dataParam] = manufacturedSolution(M,tFinal)
2222 % Notes:
2323 % Take the densities to be Gaussians for testing.
2424
25+if nargin == 0, test_manufacturedSolution(); return; end
26+
2527 override = 0; % turn on '1' to avoid restrictions on chosen test function
2628
27-mu = linspace(40,50,M)';
28-t0 = linspace(1,3,M)';
29+% mu = linspace(40,50,M)';
30+mu = linspace(2,5,M)';
31+
32+% choose t0 such that sig(0) = eps
33+t0 = sqrt(log(1/eps)./mu);
2934
3035 % choose the initial time step
3136 dt0_vec = pi./(4*sqrt(mu*log(1/eps)));
@@ -44,3 +49,15 @@ dataParam.t0 = t0;
4449 dataParam.mu = mu;
4550
4651 end
52+
53+%%%%
54+function test_manufacturedSolution()
55+
56+M = 1; tFinal = 3*pi;
57+[dt0,dataParam] = manufacturedSolution(M,tFinal);
58+
59+t = linspace(0,tFinal,1000);
60+figure(1)
61+plot(t,dataParam.sig(t)); grid on;
62+
63+end
wfp_1Dspringscattering/utils/antiderivmat_1d.mdeleted+0−37View file
@@ -1,37 +0,0 @@
1-function L = antiderivmat_1d(s)
2-% ANTIDERIVMAT_1D matrix from nodes in 1D to antiderivatives at nodes
3-%
4-% L = antiderivmat_1d(s) returns (N-1)*N matrix taking values on s,
5-% a list of N nodes, to their integrals from the first node to each of the
6-% other N-1 nodes in turn. Ie, the antiderivative with constant chosen so that
7-% its value at the first node would be zero. Nodes must be spaced in a sensible
8-% way. length(s) should not exceed around 30 for stability reasons.
9-%
10-% Notes: 1) should match Leslie's get_integratemat().
11-% 2) Computed in Helsing style, with centering and scaling for stability.
12-
13-% Barnett 8/13/24.
14-
15-if nargin==0, test_antiderivmat_1d; return; end
16-
17-cen = (max(s)+min(s))/2; hwid = (max(s)-min(s))/2; % affine map s to [-1,1]
18-s = (s-cen)/hwid;
19-n = numel(s); s = s(:); % col vec
20-V = ones(n); for j=2:n, V(:,j) = V(:,j-1).*s; end % Vandermonde (polyval) mat
21-U = diag(s)*V*diag(1./(1:n)); % mat evaluating an antideriv of poly
22-L = (V'\U')'; % backwards-stable way to solve for it (Helsing)
23-L = L(2:end,:) - L(1,:); % adjust const to zero at first node
24-L = L*hwid; % unscale
25-
26-%%%%%%
27-function test_antiderivmat_1d
28-off = 4.3; sc = 1.7; % test centering and scaling
29-x = sc*linspace(-1,1,16)' + off; % nodes
30-%x = sc*gauss(16) + off; % nodes
31-f = @(x) sin(0.8*x + 0.7); % the antiderivative
32-fp = @(x) 0.8*cos(0.8*x + 0.7); % the input func
33-Fex = f(x(2:end))-f(x(1)); % exact ans at nodes 2...N, col vec
34-L = antiderivmat_1d(x);
35-F = L * fp(x); % hit L against vec of func values
36-fprintf('max abs err for antideriv on nodes : %.3g\n',max(abs(F - Fex)))
37-fprintf('[mat inf-norm = %.3g; max element size = %.3g]\n',norm(L,inf),max(abs(L(:))))
wfp_1Dspringscattering/utils/chebEval.mmodified+1−0View file
@@ -1,3 +1,4 @@
1+
12 function y = chebEval(x,wts,chebApproxInfo)
23 % CHEBEval retruns the approximation of a function at x
34 %
wfp_1Dspringscattering/utils/interpmat_1d.mdeleted+0−38View file
@@ -1,38 +0,0 @@
1-function L = interpmat_1d(t,s)
2-% INTERPMAT_1D interpolation matrix from nodes in 1D to any target nodes
3-%
4-% L = interpmat_1d(t,s) returns interpolation matrix taking values on nodes s
5-% (a list of nodes) to target nodes t. It assumes smooth functions.
6-% length(s) should be kept small, eg, 30 or less.
7-%
8-% Run without arguments does a self test (see test code for usage example).
9-%
10-% Notes: Computed in Helsing style, with centering and scaling for stability.
11-
12-% Barnett 7/17/16. Auto-centering & scaling for stability 12/23/21.
13-
14-if nargin==0, test_interpmat_1d; return; end
15-
16-cen = (max(s)+min(s))/2; hwid = (max(s)-min(s))/2; % affine s to [-1,1]
17-s = (s-cen)/hwid;
18-t = (t-cen)/hwid;
19-
20-p = numel(s); q = numel(t); s = s(:); t = t(:); % all col vecs
21-n = p; % set the polynomial order we go up to
22-V = ones(p,n); for j=2:n, V(:,j) = V(:,j-1).*s; end % polyval matrix on nodes
23-R = ones(q,n); for j=2:n, R(:,j) = R(:,j-1).*t; end % polyval matrix on targs
24-L = (V'\R')'; % backwards-stable way to do it (Helsing) See corners/interpdemo.m
25-
26-%%%%%%
27-function test_interpmat_1d
28-off = 4.3; % test centering and scaling
29-sc = 1.7;
30-x = sc*linspace(-1,1,16)' + off;
31-f = @(x) sin(x + 0.7);
32-data = f(x); % func on smooth (src) nodes
33-t = sc*(2*rand(1000,1) - 1) + off; % cover same interval as the x lie
34-uex = f(t); % col vec
35-L = interpmat_1d(t,x);
36-u = L * data;
37-fprintf('max abs err for interp in [a,b] : %.3g\n',max(abs(u - uex)))
38-fprintf('interp mat inf-norm = %.3g; max element size = %.3g\n',norm(L,inf),max(abs(L(:))))
wfp_1Dspringscattering/maxNorm.m →wfp_1Dspringscattering/utils/maxNorm.mrenamed+0−0View file
No changes to the file's content.
wfp_1Dspringscattering/mergeWeights.m →wfp_1Dspringscattering/utils/mergeWeights.mrenamed+0−0View file
No changes to the file's content.
wfp_1Dspringscattering/utils/perispecdiff.mdeleted+0−25View file
@@ -1,25 +0,0 @@
1-function g = perispecdiff(f)
2-% PERISPECDIFF - use FFT to take periodic spectral differentiation of vector
3-%
4-% g = perispecdiff(f) returns g the derivative of the spectral interpolant
5-% of f, which is assumed to be the values of a smooth 2pi-periodic function
6-% at the N gridpoints 2.pi.j/N, for j=1,..,N (or any translation of such
7-% points). Can be row or col vec, and output is same shape.
8-%
9-% Without arguments, does a self-test.
10-
11-% Barnett 2/18/14
12-if nargin==0, test_perispecdiff; return; end
13-N = numel(f);
14-if mod(N,2)==0 % even
15- g = ifft(fft(f(:)).*[0 1i*(1:N/2-1) 0 1i*(-N/2+1:-1)].');
16-else
17- g = ifft(fft(f(:)).*[0 1i*(1:(N-1)/2) 1i*((1-N)/2:-1)].');
18-end
19-g = reshape(g,size(f));
20-
21-%%%%%%
22-function test_perispecdiff
23-N = 50; tj = 2*pi/N*(1:N)';
24-f = sin(3*tj); fp = 3*cos(3*tj); % trial periodic function & its deriv
25-norm(fp-perispecdiff(f))
wfp_1Dspringscattering/utils/vline.mdeleted+0−107View file
@@ -1,107 +0,0 @@
1-function hhh=vline(x,in1,in2)
2-% function h=vline(x, linetype, label)
3-%
4-% Draws a vertical line on the current axes at the location specified by 'x'. Optional arguments are
5-% 'linetype' (default is 'r:') and 'label', which applies a text label to the graph near the line. The
6-% label appears in the same color as the line.
7-%
8-% The line is held on the current axes, and after plotting the line, the function returns the axes to
9-% its prior hold state.
10-%
11-% The HandleVisibility property of the line object is set to "off", so not only does it not appear on
12-% legends, but it is not findable by using findobj. Specifying an output argument causes the function to
13-% return a handle to the line, so it can be manipulated or deleted. Also, the HandleVisibility can be
14-% overridden by setting the root's ShowHiddenHandles property to on.
15-%
16-% h = vline(42,'g','The Answer')
17-%
18-% returns a handle to a green vertical line on the current axes at x=42, and creates a text object on
19-% the current axes, close to the line, which reads "The Answer".
20-%
21-% vline also supports vector inputs to draw multiple lines at once. For example,
22-%
23-% vline([4 8 12],{'g','r','b'},{'l1','lab2','LABELC'})
24-%
25-% draws three lines with the appropriate labels and colors.
26-%
27-% By Brandon Kuczenski for Kensington Labs.
28-% brandon_kuczenski@kensingtonlabs.com
29-% 8 November 2001
30-
31-if length(x)>1 % vector input
32- for I=1:length(x)
33- switch nargin
34- case 1
35- linetype='r:';
36- label='';
37- case 2
38- if ~iscell(in1)
39- in1={in1};
40- end
41- if I>length(in1)
42- linetype=in1{end};
43- else
44- linetype=in1{I};
45- end
46- label='';
47- case 3
48- if ~iscell(in1)
49- in1={in1};
50- end
51- if ~iscell(in2)
52- in2={in2};
53- end
54- if I>length(in1)
55- linetype=in1{end};
56- else
57- linetype=in1{I};
58- end
59- if I>length(in2)
60- label=in2{end};
61- else
62- label=in2{I};
63- end
64- end
65- h(I)=vline(x(I),linetype,label);
66- end
67-else
68- switch nargin
69- case 1
70- linetype='r:';
71- label='';
72- case 2
73- linetype=in1;
74- label='';
75- case 3
76- linetype=in1;
77- label=in2;
78- end
79-
80-
81-
82-
83- g=ishold(gca);
84- hold on
85-
86- y=get(gca,'ylim');
87- h=plot([x x],y,linetype);
88- if length(label)
89- xx=get(gca,'xlim');
90- xrange=xx(2)-xx(1);
91- xunit=(x-xx(1))/xrange;
92- if xunit<0.8
93- text(x+0.01*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
94- else
95- text(x-.05*xrange,y(1)+0.1*(y(2)-y(1)),label,'color',get(h,'color'))
96- end
97- end
98-
99- if g==0
100- hold off
101- end
102- set(h,'tag','vline','handlevisibility','off')
103-end % else
104-
105-if nargout
106- hhh=h;
107-end