updated code
10 changed files+356−69
.DS_Storemodified+0−0View file
Binary file not shown.
wfp_1Dspringscattering/.DS_Storemodified+0−0View file
Binary file not shown.
wfp_1Dspringscattering/fixDtBasedOnTypNeighborsNum.mmodified+2−2View file
@@ -1,4 +1,4 @@
1-function [dt,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,dt)
1+function [dt,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,dt,keepdt)
22 % FIXDTBASEDONTYPNEIGHBORSNUM adjusts dt to make sure num of local
33 % neighbors is less than a maximum value
44 %
@@ -25,7 +25,7 @@ end
2525
2626 typNumOfNeighbors = ceil(2*W*dt*M/(bs - as));
2727
28-if(typNumOfNeighbors>maxNumNeighbors)
28+if(typNumOfNeighbors>maxNumNeighbors && keepdt == 0)
2929 dt = (bs - as)*maxNumNeighbors/(2*W*M);
3030 fprintf('time step changed to dt = %1.2e to ensure maximum number of neighbors %d\n',dt,maxNumNeighbors);
3131 typNumOfNeighbors = maxNumNeighbors;
wfp_1Dspringscattering/get_experimentParameters.madded+273−0View file
@@ -0,0 +1,273 @@
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,mu,dt0_min,uniform_beta,keepdt]...
5+ =get_experimentParameters(expNum)
6+
7+switch expNum
8+ case 1
9+
10+ % testing buttons
11+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
12+ saveWorkspaceOpt = 0;
13+ addErrResults = 0;
14+ savePlot = 0;
15+ logScale = 0; % '1' to use log scale in heat maps
16+ printTimeStep = 0; % '1' to print out each time step
17+ evalSol = 1;
18+
19+ % parameters
20+ order = 2; % order of accuracy = interpolation
21+ numResolutions = 3; % num of grid resolutions
22+ solnType = 'ms'; % type of solution 'ms' or 'true'
23+
24+ tFinal = 3*pi; % final time
25+
26+ %%% 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
29+
30+ %%% sources
31+ M = 2; % number of sources
32+ maxNumNeighbors = NaN; % set max number of neighbors
33+ src_dmn = [-1,1];
34+ ds = 1e-4; % set min distance between sources
35+ uniform_sgrid = 1;
36+
37+ %%% spring constants
38+ betaMax = 3; uniform_beta = 1;
39+
40+ %%% incident pulse parameter
41+ mu = 50;
42+ keepdt = 1;
43+
44+ dt0_min = 0.1;
45+ case 2
46+
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+ mu = 30;
79+ keepdt = 0;
80+
81+ dt0_min = 0.1;
82+
83+ case 3
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+ mu = M^2;
117+ keepdt = 0;
118+
119+ dt0_min = 0.1;
120+ case 4
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 = 0; % '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 = 800; % 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 = 50; % 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+ mu = M^2;
153+ keepdt = 0;
154+
155+ dt0_min = 0.1;
156+ case 5 % FP
157+ % testing buttons
158+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
159+ saveWorkspaceOpt = 1;
160+ addErrResults = 0;
161+ savePlot = 0;
162+ logScale = 0; % '1' to use log scale in heat maps
163+ printTimeStep = 0; % '1' to print out each time step
164+ evalSol = 1;
165+
166+ % parameters
167+ order = 8; % order of accuracy = interpolation
168+ numResolutions = 2; % num of grid resolutions
169+ solnType = 'true'; % type of solution 'ms' or 'true'
170+
171+ tFinal = 30*pi; % final time
172+
173+ %%% Choose number of spatial and temporal set for solution computation
174+ Nx = 400; % size of spatial grid for testing and plotting
175+ Nt_sol = 500; % size of the time grid for testing and plotting
176+
177+ %%% sources
178+ M = 2; % number of sources
179+ maxNumNeighbors = 2; % set max number of neighbors
180+ src_dmn = [-0.5,0.5];
181+ ds = 1e-4; % set min distance between sources
182+ uniform_sgrid = 1;
183+
184+ %%% spring constants
185+ betaMax = 100; uniform_beta = 1;
186+
187+ %%% incident pulse parameter
188+ mu = 5;
189+ keepdt = 0;
190+
191+ dt0_min = 0.02; % RBC constraint
192+ dt0_min = min(0.01,dt0_min); % stability constraint
193+
194+ case 6 % FP M = 10
195+ % testing buttons
196+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
197+ saveWorkspaceOpt = 0;
198+ addErrResults = 0;
199+ savePlot = 0;
200+ logScale = 0; % '1' to use log scale in heat maps
201+ printTimeStep = 1; % '1' to print out each time step
202+ evalSol = 1;
203+
204+ % parameters
205+ order = 8; % order of accuracy = interpolation
206+ numResolutions = 3; % num of grid resolutions
207+ solnType = 'true'; % type of solution 'ms' or 'true'
208+
209+ tFinal = 3*pi; % final time
210+
211+ %%% 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
214+
215+ %%% sources
216+ M = 2; % number of sources
217+ maxNumNeighbors = 10; % set max number of neighbors
218+ src_dmn = [-2,2];
219+ ds = 1e-4; % set min distance between sources
220+ uniform_sgrid = 1;
221+
222+ %%% spring constants
223+ betaMax = 100; uniform_beta = 1;
224+
225+ %%% incident pulse parameter
226+ mu = 5;
227+ keepdt = 0;
228+
229+ dt0_min = 0.02; % RBC constraint
230+ dt0_min = min(0.001,dt0_min); % stability constraint
231+ case 7 % FP M = 200
232+ % testing buttons
233+ plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
234+ saveWorkspaceOpt = 1;
235+ addErrResults = 0;
236+ savePlot = 0;
237+ logScale = 0; % '1' to use log scale in heat maps
238+ printTimeStep = 0; % '1' to print out each time step
239+ evalSol = 1;
240+
241+ % parameters
242+ order = 8; % order of accuracy = interpolation
243+ numResolutions = 2; % num of grid resolutions
244+ solnType = 'true'; % type of solution 'ms' or 'true'
245+
246+ tFinal = 40*pi; % final time
247+
248+ %%% Choose number of spatial and temporal set for solution computation
249+ Nx = 800; % size of spatial grid for testing and plotting
250+ Nt_sol = 800; % size of the time grid for testing and plotting
251+
252+ %%% sources
253+ M = 200; % number of sources
254+ maxNumNeighbors = 20; % set max number of neighbors
255+ src_dmn = [-2,2];
256+ ds = 1e-4; % set min distance between sources
257+ uniform_sgrid = 1;
258+
259+ %%% spring constants
260+ betaMax = 1; uniform_beta = 1;
261+
262+ %%% incident pulse parameter
263+ mu = 5;
264+ keepdt = 0;
265+
266+ dt0_min = 0.02; % RBC constraint
267+ dt0_min = min(0.01,dt0_min); % stability constraint
268+ otherwise
269+ disp('experiment number inapplicable');
270+end
271+
272+end
273+
wfp_1Dspringscattering/get_msOrIncidentField.madded+40−0View file
@@ -0,0 +1,40 @@
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+3−6View file
@@ -1,4 +1,4 @@
1-function [h0,dataParam] = get_uIncidentInfo(mu)
1+function [dt0,dataParam] = get_uIncidentInfo(mu)
22 % GET_UINCIDENTINFO prepares information needed for an incident wave of the
33 % form uin = exp(-mu*(x - t - t0).^2);
44 %
@@ -16,11 +16,8 @@ function [h0,dataParam] = get_uIncidentInfo(mu)
1616
1717 t0 = -3;
1818
19-% Choose h0 such that the densities are resolved
20-N0 = 30; % number of grid points per standard deviation (sd)
21-bumpWidth = sqrt(log(1/eps)./mu);
22-h0_vec = bumpWidth./N0;
23-h0 = min(h0_vec);
19+% choose the initial time step
20+dt0 = pi/(4*sqrt(mu*log(1/eps)));
2421
2522 dataParam.mu = mu;
2623 dataParam.t0 = t0;
wfp_1Dspringscattering/get_xgrid.mmodified+1−0View file
@@ -17,6 +17,7 @@ h_temp = h0;
1717 xmin = (ax + 2*h_temp*W); xmax = (bx - 2*h_temp*W);
1818 x = linspace(xmin,xmax,Nx)';
1919 if(xmin>as || xmax<bs)
20+ xmin,xmax,h0
2021 error('Sources are outside allowed region; reduce the time-step or narrow down the domain where sources live');
2122 end
2223
wfp_1Dspringscattering/main.mmodified+31−45View file
@@ -2,67 +2,49 @@ startMatlabFile
22
33 %% Set file directories to save data, figures, tables and movies
44 addpath ./utils;
5-dataFile = './data';
5+dataFile = '/Users/nalhassanieh/Desktop/WFP_results/data';
66 figFile = './fig';
77
8-%% Testing Buttons
9-plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
10-saveWorkspaceOpt = 0;
11-addErrResults = 0;
12-savePlot = 0;
13-logScale = 0; % '1' to use log scale in heat maps
14-printTimeStep = 0; % '1' to print out each time step
15-evalSol = 1;
16-
178 %% Problem parameters
18-%%% order and test solution
19-order = 2; % order of accuracy = interpolation
20-numResolutions = 3; % num of grid resolutions
21-solnType = 'ms'; % type of solution 'ms' or 'true'
9+expNum = 6; % experiment number
10+
11+[plotOpt,saveWorkspaceOpt,addErrResults,savePlot,...
12+ logScale,printTimeStep,evalSol,order,numResolutions,solnType,tFinal,...
13+ Nx,Nt_sol,M,maxNumNeighbors,src_dmn,ds,...
14+ uniform_sgrid,betaMax,mu,dt0_min,uniform_beta,keepdt]...
15+ = get_experimentParameters(expNum);
2216
23-%%% domain and final time
17+%%% domain
2418 ax = -pi; bx = pi; % Domain
25-tFinal = 3*pi; % final time
2619 tsStop = -1; % stop time stepping at tsStop. ('-1' continue to tFinal)
2720
28-%%% Choose number of spatial and temporal set for solution computation
29-Nx = 10; % size of spatial grid for testing and plotting
30-Nt_sol = 10; % size of the time grid for testing and plotting
31-tSOL = tFinal/Nt_sol; % time-step for the solution evaluation
32-
33-%%% sources
34-M = 10; % number of sources
35-maxNumNeighbors = 10; % set max number of neighbors
36-src_dmn = [-1,1];
37-ds = 1e-4; % set min distance between sources
38-uniform_sgrid = 0;
21+%%% sources and spring constants
3922 s = get_sourceGrid(src_dmn,ds,M,uniform_sgrid);
40-
41-%%% spring constants
42-betaMax = 3; uniform_beta = 0;
4323 beta = get_springConstants(betaMax,M,uniform_beta);
4424
45-%%% incident pulse parameter
46-mu = 30;
25+%%% Choose number of spatial and temporal set for solution computation
26+tSOL = tFinal/Nt_sol; % time-step for the solution evaluation
4727
4828 %%% window
4929 tol = 1e-12; % error tolerance
5030 [winData,W] = setup_generalwindow(tol);
5131
52-%% Compute some values for the window function and GL
53-P = W; % GL nodes
54-gl = glwt_prep(P); % GL nodes and weights on [-1,1]
55-
5632 %% Manufactured density values for testing or incident wave for true scattering
57-if(strcmp(solnType,'ms'));[h0,dataParam] = manufacturedSolution(M,tFinal);
58-else; [h0,dataParam] = get_uIncidentInfo(mu); end
33+if(strcmp(solnType,'ms'));[dt0,dataParam] = manufacturedSolution(M,tFinal);
34+else; [dt0,dataParam] = get_uIncidentInfo(mu); end
35+dt0 = min(dt0_min,dt0); % domain RBC constraint
5936
6037 % fix initial time step
61-[h0,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,h0);
62-Nt0 = ceil(tFinal/h0); if(mod(Nt0,2) == 1); Nt0 = Nt0 + 1; end
38+[dt0,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,dt0,keepdt);
39+Nt0 = 2*ceil(0.5*tFinal/dt0);
40+dt0 = tFinal/Nt0;
41+
42+%% Compute some values for the window function and GL
43+P = W; % GL nodes
44+gl = glwt_prep(P); % GL nodes and weights on [-1,1]
6345
6446 % get the spatial grid
65-x = get_xgrid(ax,bx,src_dmn,Nx,h0,W);
47+x = get_xgrid(ax,bx,src_dmn,Nx,dt0,W);
6648
6749 %% Print title
6850 printTitle(tFinal, tol, W, P, M, typNumOfNeighbors, order, solnType);
@@ -76,14 +58,17 @@ if(evalSol == 1)
7658 t_cell = cell(1,numResolutions);
7759 x_cell = cell(1,numResolutions);
7860 sn_cell = cell(1,numResolutions);
61+ tm_cell = cell(1,numResolutions);
7962 end
8063
8164 for m = 1:numResolutions
82- N = Nt0*(2^(m-1));
8365
84- dt = 2*pi/N;
85- Nt = ceil(tFinal/dt);
86- dt = tFinal/Nt; % timestep
66+ dt = dt0/(2^(m-1));
67+ Nt = Nt0*2^(m-1); % number of time steps
68+
69+ N = 2*ceil(0.5*2*(pi/dt)); % N = number of Fourier modes case gamma = 0.5;
70+ % ensure N is even
71+
8772 h_RBC = 2*pi/N;
8873
8974 tgh = W + order;
@@ -199,6 +184,7 @@ for m = 1:numResolutions
199184 t_cell{m} = tn_sol;
200185 x_cell{m} = x';
201186 sn_cell{m} = sn_sol;
187+ tm_cell{m} = tm;
202188
203189 % Check the error at all time
204190 if(strcmp(solnType,'ms'))
wfp_1Dspringscattering/manufacturedSolution.mmodified+4−14View file
@@ -1,4 +1,4 @@
1-function [h0,dataParam] = manufacturedSolution(M,tFinal)
1+function [dt0,dataParam] = manufacturedSolution(M,tFinal)
22 % MANUFACTUREDSOLUTION analytic density functions
33 %
44 % [sig,Nx0,h0,t0,mu] = manufacturedSolution(M,tFinal,s,beta)
@@ -27,19 +27,9 @@ override = 0; % turn on '1' to avoid restrictions on chosen test function
2727 mu = linspace(40,50,M)';
2828 t0 = linspace(1,3,M)';
2929
30-% choose t0 such that sig(0) = 0.01*eps
31-% t0 = sqrt(log(amp/(0.01*eps))./mu);
32-
33-%%% Working resolution using standard deviation
34-N0 = 20; % number of grid points per standard deviation (sd)
35-bumpWidth = sqrt(log(1/eps)./mu);
36-h0_vec = bumpWidth./N0;
37-h0 = min(h0_vec);
38-
39-%%% Resolution computed from Fourier analysis
40-% omega0 = 2*sqrt(mu).*sqrt(log((1/tol)*sqrt(pi./mu)));
41-% h0 = min(pi*(1 - gam)./omega0);
42-% Nt0 = ceil(tFinal/h0); % or Nt0 = ceil(2*pi/h0);
30+% choose the initial time step
31+dt0_vec = pi./(4*sqrt(mu*log(1/eps)));
32+dt0 = min(dt0_vec);
4333
4434 % ensure that the forcing function is bandlimited within the computational time domain
4535 if(any(2*t0>tFinal)&&(override == 0))
wfp_1Dspringscattering/saveWorkspace.mmodified+2−2View file
@@ -4,7 +4,8 @@ if (uniform_sgrid == 1)
44
55 % Save only the selected variables
66 fileName = sprintf('%s/WFPworkspaceO%d_T%d_%sM%d_b%d',dataFile,order,round(tFinal),str,M,betaMax);
7-
7+ matlabDataFile =strcat(fileName,'.mat');
8+
89 if (saveWorkspaceOpt == 1)
910 % all variables in the workspace
1011 vars = who;
@@ -24,7 +25,6 @@ if (saveWorkspaceOpt == 1)
2425 fprintf('saved workspace to %s\n',matlabDataFile);
2526
2627 elseif(addErrResults == 1)
27- matlabDataFile =strcat(fileName,'.mat');
2828 if(isfile(matlabDataFile))
2929 save(fileName, 'sc_err', '-append');
3030 fprintf('added error to the data file in %s\n',matlabDataFile);