1startMatlabFile
3%% Set file directories to save data, figures, tables and movies
4addpath ./utils;
5dataFile = './data';
6figFile = './fig';
8%% Testing Buttons
9plotOpt = 0; % '1' to plot sol at final time, '2' heat maps
10saveWorkspaceOpt = 0;
11addErrResults = 0;
12savePlot = 0;
13logScale = 0; % '1' to use log scale in heat maps
14printTimeStep = 0; % '1' to print out each time step
15evalSol = 1;
17%% Problem parameters
18%%% order and test solution
19order = 2; % order of accuracy = interpolation
20numResolutions = 3; % num of grid resolutions
21solnType = 'ms'; % type of solution 'ms' or 'true'
23%%% domain and final time
24ax = -pi; bx = pi; % Domain
25tFinal = 3*pi; % final time
26tsStop = -1; % stop time stepping at tsStop. ('-1' continue to tFinal)
28%%% Choose number of spatial and temporal set for solution computation
29Nx = 10; % size of spatial grid for testing and plotting
30Nt_sol = 10; % size of the time grid for testing and plotting
31tSOL = tFinal/Nt_sol; % time-step for the solution evaluation
33%%% sources
34M = 10; % number of sources
35maxNumNeighbors = 10; % set max number of neighbors
36src_dmn = [-1,1];
37ds = 1e-4; % set min distance between sources
38uniform_sgrid = 0;
39s = get_sourceGrid(src_dmn,ds,M,uniform_sgrid);
41%%% spring constants
42betaMax = 3; uniform_beta = 0;
43beta = get_springConstants(betaMax,M,uniform_beta);
45%%% incident pulse parameter
46mu = 30;
48%%% window
49tol = 1e-12; % error tolerance
50[winData,W] = setup_generalwindow(tol);
52%% Compute some values for the window function and GL
53P = W; % GL nodes
54gl = glwt_prep(P); % GL nodes and weights on [-1,1]
56%% Manufactured density values for testing or incident wave for true scattering
57if(strcmp(solnType,'ms'));[h0,dataParam] = manufacturedSolution(M,tFinal);
58else; [h0,dataParam] = get_uIncidentInfo(mu); end
60% fix initial time step
61[h0,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,h0);
62Nt0 = ceil(tFinal/h0); if(mod(Nt0,2) == 1); Nt0 = Nt0 + 1; end
64% get the spatial grid
65x = get_xgrid(ax,bx,src_dmn,Nx,h0,W);
67%% Print title
68printTitle(tFinal, tol, W, P, M, typNumOfNeighbors, order, solnType);
69prechar = sprintf('%s%d_%d_mn%d_',solnType,betaMax,dataParam.mu(1),maxNumNeighbors); % char/string to precede title of fig or tab
71%% Grid resolution study
72if(evalSol == 1)
73 err = zeros(1,numResolutions);
74 h = zeros(1,numResolutions);
75 ug = cell(1,numResolutions);
76 t_cell = cell(1,numResolutions);
77 x_cell = cell(1,numResolutions);
78 sn_cell = cell(1,numResolutions);
79end
81for m = 1:numResolutions
82 N = Nt0*(2^(m-1));
84 dt = 2*pi/N;
85 Nt = ceil(tFinal/dt);
86 dt = tFinal/Nt; % timestep
87 h_RBC = 2*pi/N;
89 tgh = W + order;
90 tInitial = -tgh*dt; % initial time
91 Ntg = Nt + 1 + tgh; % total number of time-steps
92 tn = linspace(tInitial ,tFinal, Ntg)'; % time grid
93 it1 = tgh + 1; it2 = Ntg; % interior time indices (start from zero)
94 It = it1:it2; % time indices starting from zero
95 delta = dt*W;
97 %%% time parameter for self convergence study
98 if(m == 1); fixSolCnt = ceil((Nt + 1)./Nt_sol);
99 else; fixSolCnt = 2*fixSolCnt; end
100 Nt_sol = ceil((Ntg-tgh)/fixSolCnt);
102 % frequency paramenters
103 K0 = N/2; % max wave number
104 K = [0:(K0-1),(-K0:-1)]; % frequency set
105 zeroLoc = 1; % where k = 0 is located
107 % set up tRBC
108 tRBC = W*h_RBC; % choose tFinal to be a multiple of pi
110 %% Allocate space for time-stepping
111 % let sn hold the density function values at any time tn
112 timeLevels = tgh + 2;
113 sn = zeros(timeLevels,M);% store tgh + current time + next time
114 snIdxNow = tgh+1; % index of sn at the current time
116 % Let an hold the values of the Fourier coefficients, and bn
117 % correspond to time derivative of an, needed in the history treatment
118 an = zeros(N,1);
119 bn = an; anp1 = an; bnp1 = bn;
121 % Allocate space for computed solution for heat maps
122 u = zeros(Nx,Nt_sol); % u(x,t) to generate heat map
123 ue = zeros(Nx,Nt_sol);
124 tn_sol = zeros(1,Nt_sol);
125 sn_sol = zeros(M,Nt_sol);
127 %% Prepare for time-stepping
128 t = dt; % this represents the current time
130 tic
131 % prepare nodes and weights needed for the evaluation of self local
132 % integrals and other local integrals
133 [implicitMat,A_sp,deltaInteractions,dtInteractions] = ...
134 prepNodesAndWeights(tn(It),dt,beta,gl,s,src_dmn,snIdxNow,P,M,W,order,...
135 timeLevels,typNumOfNeighbors,winData);
137 % prepare the nodes and weights for the solution local evaluation
138 if(evalSol == 1)
139 [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = prepForLocalSolEval(x,dt,gl,s,M,P,W,order,winData);
140 end
142 % prepare integrals needed in the evaluation of Fourier coefficients
143 [p,q,p0,q0,sn_hat,phi_tn,phit_tn] = prepValuesForFourierCoefs(sn,snIdxNow,dt,K,N,W,gl,s,tol,winData);
145 tm.prep = toc;
146 %% Time Stepping
147 tRBC_now = 0; tSOL_now = tSOL; n_sol = 1; totalSteps = 1;
148 tic
149 for n = It(2:(end-1)) % note: t(n) = tInitial + (n-1)*dt = (-tgh + (n-1))*dt
151 % get the history contribution
152 [anp1,bnp1,historySum,tRBC_now] = get_historyContribution(t,an,bn,sn,zeroLoc,K,...
153 dt, tRBC, W, snIdxNow, p,q,p0,q0,sn_hat,tRBC_now,phi_tn,phit_tn,s,tol);
155 % Get the gn for the spring scattering conditions
156 gn = get_g(t,dt,s,dataParam,M,beta,solnType);
158 % Perform a sparse matvec to evaluate integrals
159 sn_v = reshape(sn,[],1);
160 RHS = gn + (beta/2).*((A_sp*sn_v)+ (1/pi)*historySum);
162 % solve the implicit system
163 sn(end,:) = implicitMat\RHS;
165 % Take the real part of the density
166 sn = real(sn);
168 % new time
169 t = tn(n+1);
171 if(evalSol == 1)
172 [u,ue,tn_sol,n_sol,tSOL_now,sn_sol] = get_solnAndErr(u,ue,tn_sol,n_sol,...
173 tSOL_now,x,t,snIdxNow,sn,anp1,interpNodesIdx_sol,...
174 interpAndGLWeights_sol,interpShift_sol,s,M,delta,tol,...
175 tSOL,solnType,dataParam,fixSolCnt,n,tgh,sn_sol);
176 end
178 % update values for the next time step
179 [an,bn,sn,sn_hat] = prepForNextStep(anp1,bnp1,sn,sn_hat,snIdxNow,tgh,s,N,tol);
181 if(printTimeStep == 1)
182 fprintf('t = %1.4e\n',t);
183 end
185 if(totalSteps == tsStop)
186 break;
187 end
188 totalSteps = totalSteps + 1;
189 end
191 timeSteppingTime = toc;
192 tm.ts = timeSteppingTime/totalSteps;
194 %% Error Analysis
195 if(evalSol == 1)
196 u = real(u);
197 ug{m} = u;
198 h(m) = dt;
199 t_cell{m} = tn_sol;
200 x_cell{m} = x';
201 sn_cell{m} = sn_sol;
203 % Check the error at all time
204 if(strcmp(solnType,'ms'))
205 err(m) = max(max(abs(u-ue))); % max-norm error
206 rate = printOutResults(m,t,Nt,K0,dt,err,h);
207 end
209 % Plot and save into a file if needed
210 if(plotOpt == 1 && strcmp(solnType,'ms'))
211 plotSolnAndError(t,x,u(:,end),ue,figFile,prechar,order,tFinal,savePlot);
212 pause(.1);
213 elseif(plotOpt == 2)
214 plotHeatMap(u,ue,x,tn_sol,s,figFile,order,tFinal,prechar,savePlot,solnType,logScale,dataFile);
215 pause(1);
216 end
217 end
218end
220checkMatchingGrids(t_cell,x_cell,numResolutions);
222if(numResolutions>2)
223 [sc_err,sc_rate] = selfConvergence(ug,h,numResolutions);
224elseif(numResolutions==2)
225 sc_err = maxNorm(ug{2} - ug{1});
226 fprintf('self convergence error = %1.2e\n',sc_err);
227end
229if(saveWorkspaceOpt == 1 || addErrResults == 1)
230 saveWorkspace;
231end