1startMatlabFile
3%% Set file directories to save data, figures, tables and movies
4addpath ./utils;
5dataFile = './data';
6figFile = './fig';
8%% Problem parameters
9expNum = 1; % experiment number
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,dataParam,dt0_min,uniform_beta,keepdt]...
15 = get_experimentParameters(expNum);
17%%% domain
18ax = -pi; bx = pi; % Domain
19tsStop = -1; % stop time stepping at tsStop. ('-1' continue to tFinal)
21%%% sources and spring constants
22s = get_sourceGrid(src_dmn,ds,M,uniform_sgrid);
23beta = get_springConstants(betaMax,M,uniform_beta);
25%%% Choose number of spatial and temporal set for solution computation
26tSOL = tFinal/Nt_sol; % time-step for the solution evaluation
28%%% window
29tol = 1e-12; % error tolerance
30[winData,W] = setup_generalwindow(tol);
32%% Manufactured density values for testing or incident wave for true scattering
33if(strcmp(solnType,'ms'));[dt0,dataParam] = manufacturedSolution(M,tFinal);
34else; [dt0,dataParam] = get_uIncidentInfo(dataParam); end
35dt0 = min(dt0_min,dt0); % domain RBC constraint
37% fix initial time step
38[dt0,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,dt0,keepdt);
39Nt0 = 2*ceil(0.5*tFinal/dt0);
40dt0 = tFinal/Nt0;
42%% Compute some values for the window function and GL*
43P = W; % GL nodes
44gl = glwt_prep(P); % GL nodes and weights on [-1,1]
46% get the spatial grid
47x = get_xgrid(ax,bx,src_dmn,Nx,dt0,W);
49%% Print title
50printTitle(tFinal, tol, W, P, M, typNumOfNeighbors, order, solnType);
51prechar = sprintf('%s%d_%d_mn%d_',solnType,betaMax,dataParam.mu(1),maxNumNeighbors); % char/string to precede title of fig or tab
53%% Grid resolution study
54if(evalSol == 1)
55 err = zeros(1,numResolutions);
56 h = zeros(1,numResolutions);
57 ug = cell(1,numResolutions);
58 t_cell = cell(1,numResolutions);
59 x_cell = cell(1,numResolutions);
60 sn_cell = cell(1,numResolutions);
61 tm_cell = cell(1,numResolutions);
62end
64for m = 1:numResolutions
66 dt = dt0/(2^(m-1));
67 Nt = Nt0*2^(m-1); % number of time steps
69 N = 2*ceil(0.5*2*(pi/dt)); % N = number of Fourier modes case gamma = 0.5;
70 % ensure N is even
72 h_RBC = 2*pi/N;
74 tgh = W + order;
75 tInitial = -tgh*dt; % initial time
76 Ntg = Nt + 1 + tgh; % total number of time-steps
77 tn = linspace(tInitial ,tFinal, Ntg)'; % time grid
78 it1 = tgh + 1; it2 = Ntg; % interior time indices (start from zero)
79 It = it1:it2; % time indices starting from zero
80 delta = dt*W;
82 %%% time parameter for self convergence study
83 if(m == 1); fixSolCnt = ceil((Nt + 1)./Nt_sol);
84 else; fixSolCnt = 2*fixSolCnt; end
85 Nt_sol = ceil((Ntg-tgh)/fixSolCnt);
87 % frequency paramenters
88 K0 = N/2; % max wave number
89 K = [0:(K0-1),(-K0:-1)]; % frequency set
90 zeroLoc = 1; % where k = 0 is located
92 % set up tRBC
93 tRBC = W*h_RBC; % choose tFinal to be a multiple of pi
95 %% Allocate space for time-stepping
96 % let sn hold the density function values at any time tn
97 timeLevels = tgh + 2;
98 sn = zeros(timeLevels,M);% store tgh + current time + next time
99 snIdxNow = tgh+1; % index of sn at the current time
101 % Let an hold the values of the Fourier coefficients, and bn
102 % correspond to time derivative of an, needed in the history treatment
103 an = zeros(N,1);
104 bn = an; anp1 = an; bnp1 = bn;
106 % Allocate space for computed solution for heat maps
107 u = zeros(Nx,Nt_sol); % u(x,t) to generate heat map
108 ue = zeros(Nx,Nt_sol);
109 tn_sol = zeros(1,Nt_sol);
110 sn_sol = zeros(M,Nt_sol);
112 %% Prepare for time-stepping
113 t = dt; % this represents the current time
115 tic
116 % prepare nodes and weights needed for the evaluation of self local
117 % integrals and other local integrals
118 [implicitMat,A_sp,deltaInteractions,dtInteractions] = ...
119 prepNodesAndWeights(tn(It),dt,beta,gl,s,src_dmn,snIdxNow,P,M,W,order,...
120 timeLevels,typNumOfNeighbors,winData);
122 % prepare the nodes and weights for the solution local evaluation
123 if(evalSol == 1)
124 [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = prepForLocalSolEval(x,dt,gl,s,M,P,W,order,winData);
125 end
127 % prepare integrals needed in the evaluation of Fourier coefficients
128 [p,q,p0,q0,sn_hat,phi_tn,phit_tn] = prepValuesForFourierCoefs(sn,snIdxNow,dt,K,N,W,gl,s,tol,winData);
130 tm.prep = toc;
131 %% Time Stepping
132 tRBC_now = 0; tSOL_now = tSOL; n_sol = 1; totalSteps = 1;
133 tic
134 for n = It(2:(end-1)) % note: t(n) = tInitial + (n-1)*dt = (-tgh + (n-1))*dt
136 % get the history contribution
137 [anp1,bnp1,historySum,tRBC_now] = get_historyContribution(t,an,bn,sn,zeroLoc,K,...
138 dt, tRBC, W, snIdxNow, p,q,p0,q0,sn_hat,tRBC_now,phi_tn,phit_tn,s,tol);
140 % Get the gn for the spring scattering conditions
141 gn = get_g(t,dt,s,dataParam,M,beta,solnType);
143 % Perform a sparse matvec to evaluate integrals
144 sn_v = reshape(sn,[],1);
145 RHS = gn + (beta/2).*((A_sp*sn_v)+ (1/pi)*historySum);
147 % solve the implicit system
148 sn(end,:) = implicitMat\RHS;
150 % Take the real part of the density
151 sn = real(sn);
153 % new time
154 t = tn(n+1);
156 if(evalSol == 1)
157 [u,ue,tn_sol,n_sol,tSOL_now,sn_sol] = get_solnAndErr(u,ue,tn_sol,n_sol,...
158 tSOL_now,x,t,snIdxNow,sn,anp1,interpNodesIdx_sol,...
159 interpAndGLWeights_sol,interpShift_sol,s,M,delta,tol,...
160 tSOL,solnType,dataParam,fixSolCnt,n,tgh,sn_sol);
161 end
163 % update values for the next time step
164 [an,bn,sn,sn_hat] = prepForNextStep(anp1,bnp1,sn,sn_hat,snIdxNow,tgh,s,N,tol);
166 if(printTimeStep == 1 && mod(totalSteps,100) == 0)
167 fprintf('t = %1.4e\n',t);
168 end
170 if(totalSteps == tsStop)
171 break;
172 end
173 totalSteps = totalSteps + 1;
174 end
176 timeSteppingTime = toc;
177 tm.ts = timeSteppingTime/totalSteps;
179 %% Error Analysis
180 if(evalSol == 1)
181 u = real(u);
182 ug{m} = u;
183 h(m) = dt;
184 t_cell{m} = tn_sol;
185 x_cell{m} = x';
186 sn_cell{m} = sn_sol;
187 tm_cell{m} = tm;
189 % Check the error at all time
190 if(strcmp(solnType,'ms'))
191 err(m) = max(max(abs(u-ue))); % max-norm error
192 rate = printOutResults(m,t,Nt,K0,dt,err,h);
193 end
195 % Plot and save into a file if needed
196 if(plotOpt == 1 && strcmp(solnType,'ms'))
197 plotSolnAndError(t,x,u(:,end),ue,figFile,prechar,order,tFinal,savePlot);
198 pause(.1);
199 elseif(plotOpt == 2)
200 plotHeatMap(u,ue,x,tn_sol,s,figFile,order,tFinal,prechar,savePlot,solnType,logScale,dataFile);
201 pause(1);
202 end
203 end
204end
206checkMatchingGrids(t_cell,x_cell,numResolutions);
208if(numResolutions>2)
209 [sc_err,sc_rate] = selfConvergence(ug,h,numResolutions);
210elseif(numResolutions==2)
211 sc_err = maxNorm(ug{2} - ug{1});
212 fprintf('self convergence error = %1.2e\n',sc_err);
213end
215if(saveWorkspaceOpt == 1 || addErrResults == 1)
216 saveWorkspace;
217end