580df75mip load finufft in main.mJeremy Magland 1mip load --install flatironinstitute/flatironinstitute/finufft
5%% Set file directories to save data, figures, tables and movies
6addpath ./utils;
10%% Problem parameters
ae869c2updated codenalhassanieh 12
13[plotOpt,saveWorkspaceOpt,addErrResults,savePlot,...
14 logScale,printTimeStep,evalSol,order,numResolutions,solnType,tFinal,...
15 Nx,Nt_sol,M,maxNumNeighbors,src_dmn,ds,...
520bc24updated readme file and cleaned up codenalhassanieh 16 uniform_sgrid,betaMax,dataParam,dt0_min,uniform_beta,keepdt]...
8748c66added codenalhassanieh 18
21tsStop = -1; % stop time stepping at tsStop. ('-1' continue to tFinal)
25beta = get_springConstants(betaMax,M,uniform_beta);
ae869c2updated codenalhassanieh 27%%% Choose number of spatial and temporal set for solution computation
28tSOL = tFinal/Nt_sol; % time-step for the solution evaluation
8748c66added codenalhassanieh 29
30%%% window
31tol = 1e-12; % error tolerance
32[winData,W] = setup_generalwindow(tol);
34%% Manufactured density values for testing or incident wave for true scattering
ae869c2updated codenalhassanieh 35if(strcmp(solnType,'ms'));[dt0,dataParam] = manufacturedSolution(M,tFinal);
520bc24updated readme file and cleaned up codenalhassanieh 36else; [dt0,dataParam] = get_uIncidentInfo(dataParam); end
8748c66added codenalhassanieh 38
39% fix initial time step
ae869c2updated codenalhassanieh 40[dt0,typNumOfNeighbors] = fixDtBasedOnTypNeighborsNum(maxNumNeighbors,src_dmn,W,M,dt0,keepdt);
41Nt0 = 2*ceil(0.5*tFinal/dt0);
42dt0 = tFinal/Nt0;
520bc24updated readme file and cleaned up codenalhassanieh 44%% Compute some values for the window function and GL*
46gl = glwt_prep(P); % GL nodes and weights on [-1,1]
8748c66added codenalhassanieh 47
48% get the spatial grid
8748c66added codenalhassanieh 50
51%% Print title
52printTitle(tFinal, tol, W, P, M, typNumOfNeighbors, order, solnType);
53prechar = sprintf('%s%d_%d_mn%d_',solnType,betaMax,dataParam.mu(1),maxNumNeighbors); % char/string to precede title of fig or tab
55%% Grid resolution study
56if(evalSol == 1)
57 err = zeros(1,numResolutions);
58 h = zeros(1,numResolutions);
59 ug = cell(1,numResolutions);
60 t_cell = cell(1,numResolutions);
61 x_cell = cell(1,numResolutions);
62 sn_cell = cell(1,numResolutions);
66for m = 1:numResolutions
69 Nt = Nt0*2^(m-1); % number of time steps
71 N = 2*ceil(0.5*2*(pi/dt)); % N = number of Fourier modes case gamma = 0.5;
72 % ensure N is even
76 tgh = W + order;
77 tInitial = -tgh*dt; % initial time
78 Ntg = Nt + 1 + tgh; % total number of time-steps
79 tn = linspace(tInitial ,tFinal, Ntg)'; % time grid
80 it1 = tgh + 1; it2 = Ntg; % interior time indices (start from zero)
81 It = it1:it2; % time indices starting from zero
82 delta = dt*W;
84 %%% time parameter for self convergence study
85 if(m == 1); fixSolCnt = ceil((Nt + 1)./Nt_sol);
86 else; fixSolCnt = 2*fixSolCnt; end
87 Nt_sol = ceil((Ntg-tgh)/fixSolCnt);
89 % frequency paramenters
90 K0 = N/2; % max wave number
91 K = [0:(K0-1),(-K0:-1)]; % frequency set
92 zeroLoc = 1; % where k = 0 is located
94 % set up tRBC
95 tRBC = W*h_RBC; % choose tFinal to be a multiple of pi
97 %% Allocate space for time-stepping
98 % let sn hold the density function values at any time tn
99 timeLevels = tgh + 2;
100 sn = zeros(timeLevels,M);% store tgh + current time + next time
101 snIdxNow = tgh+1; % index of sn at the current time
103 % Let an hold the values of the Fourier coefficients, and bn
104 % correspond to time derivative of an, needed in the history treatment
105 an = zeros(N,1);
106 bn = an; anp1 = an; bnp1 = bn;
108 % Allocate space for computed solution for heat maps
109 u = zeros(Nx,Nt_sol); % u(x,t) to generate heat map
110 ue = zeros(Nx,Nt_sol);
111 tn_sol = zeros(1,Nt_sol);
112 sn_sol = zeros(M,Nt_sol);
114 %% Prepare for time-stepping
115 t = dt; % this represents the current time
117 tic
118 % prepare nodes and weights needed for the evaluation of self local
119 % integrals and other local integrals
120 [implicitMat,A_sp,deltaInteractions,dtInteractions] = ...
121 prepNodesAndWeights(tn(It),dt,beta,gl,s,src_dmn,snIdxNow,P,M,W,order,...
122 timeLevels,typNumOfNeighbors,winData);
124 % prepare the nodes and weights for the solution local evaluation
125 if(evalSol == 1)
126 [interpNodesIdx_sol,interpAndGLWeights_sol,interpShift_sol] = prepForLocalSolEval(x,dt,gl,s,M,P,W,order,winData);
127 end
129 % prepare integrals needed in the evaluation of Fourier coefficients
130 [p,q,p0,q0,sn_hat,phi_tn,phit_tn] = prepValuesForFourierCoefs(sn,snIdxNow,dt,K,N,W,gl,s,tol,winData);
132 tm.prep = toc;
133 %% Time Stepping
134 tRBC_now = 0; tSOL_now = tSOL; n_sol = 1; totalSteps = 1;
135 tic
136 for n = It(2:(end-1)) % note: t(n) = tInitial + (n-1)*dt = (-tgh + (n-1))*dt
138 % get the history contribution
139 [anp1,bnp1,historySum,tRBC_now] = get_historyContribution(t,an,bn,sn,zeroLoc,K,...
140 dt, tRBC, W, snIdxNow, p,q,p0,q0,sn_hat,tRBC_now,phi_tn,phit_tn,s,tol);
142 % Get the gn for the spring scattering conditions
143 gn = get_g(t,dt,s,dataParam,M,beta,solnType);
145 % Perform a sparse matvec to evaluate integrals
146 sn_v = reshape(sn,[],1);
147 RHS = gn + (beta/2).*((A_sp*sn_v)+ (1/pi)*historySum);
149 % solve the implicit system
150 sn(end,:) = implicitMat\RHS;
152 % Take the real part of the density
153 sn = real(sn);
155 % new time
156 t = tn(n+1);
158 if(evalSol == 1)
159 [u,ue,tn_sol,n_sol,tSOL_now,sn_sol] = get_solnAndErr(u,ue,tn_sol,n_sol,...
160 tSOL_now,x,t,snIdxNow,sn,anp1,interpNodesIdx_sol,...
161 interpAndGLWeights_sol,interpShift_sol,s,M,delta,tol,...
162 tSOL,solnType,dataParam,fixSolCnt,n,tgh,sn_sol);
163 end
165 % update values for the next time step
166 [an,bn,sn,sn_hat] = prepForNextStep(anp1,bnp1,sn,sn_hat,snIdxNow,tgh,s,N,tol);
520bc24updated readme file and cleaned up codenalhassanieh 168 if(printTimeStep == 1 && mod(totalSteps,100) == 0)
170 end
172 if(totalSteps == tsStop)
173 break;
174 end
175 totalSteps = totalSteps + 1;
176 end
178 timeSteppingTime = toc;
179 tm.ts = timeSteppingTime/totalSteps;
181 %% Error Analysis
182 if(evalSol == 1)
183 u = real(u);
184 ug{m} = u;
185 h(m) = dt;
186 t_cell{m} = tn_sol;
187 x_cell{m} = x';
188 sn_cell{m} = sn_sol;
8748c66added codenalhassanieh 190
191 % Check the error at all time
192 if(strcmp(solnType,'ms'))
193 err(m) = max(max(abs(u-ue))); % max-norm error
194 rate = printOutResults(m,t,Nt,K0,dt,err,h);
195 end
197 % Plot and save into a file if needed
198 if(plotOpt == 1 && strcmp(solnType,'ms'))
199 plotSolnAndError(t,x,u(:,end),ue,figFile,prechar,order,tFinal,savePlot);
200 pause(.1);
201 elseif(plotOpt == 2)
202 plotHeatMap(u,ue,x,tn_sol,s,figFile,order,tFinal,prechar,savePlot,solnType,logScale,dataFile);
203 pause(1);
204 end
205 end
206end
208checkMatchingGrids(t_cell,x_cell,numResolutions);
210if(numResolutions>2)
211 [sc_err,sc_rate] = selfConvergence(ug,h,numResolutions);
212elseif(numResolutions==2)
213 sc_err = maxNorm(ug{2} - ug{1});
214 fprintf('self convergence error = %1.2e\n',sc_err);
215end
217if(saveWorkspaceOpt == 1 || addErrResults == 1)
218 saveWorkspace;