Coverage for pySDC/projects/Second_orderSDC/stability_simulation.py: 90%

1import numpy as np 

2import matplotlib.pyplot as plt 

3from pySDC.core.Errors import ProblemError 

4from pySDC.core.Step import step 

5 

6from pySDC.projects.Second_orderSDC.plot_helper import set_fixed_plot_params 

7 

8 

9class StabilityImplementation: 

10 """ 

11 This class computes and implements stability region of the harmonic oscillator problem 

12 by using different methods (SDC, Picard, RKN). 

13 

14 Parameters 

15 ----------- 

16 description: gets default paramets for the problem class 

17 kappa_max: maximum value of kappa can reach 

18 mu_max: maximum value of mu can reach 

19 Num_iter: maximum iterations for the kappa and mu on the x and y axes 

20 cwd: current working 

21 

22 """ 

23 

24 def __init__(self, description, kappa_max=20, mu_max=20, Num_iter=(400, 400), cwd=''): 

25 self.description = description 

26 self.kappa_max = kappa_max 

27 self.mu_max = mu_max 

28 self.kappa_iter = Num_iter[0] 

29 self.mu_iter = Num_iter[1] 

30 self.lambda_kappa = np.linspace(0.0, self.kappa_max, self.kappa_iter) 

31 self.lambda_mu = np.linspace(1e-10, self.mu_max, self.mu_iter) 

32 

33 self.K_iter = description['step_params']['maxiter'] 

34 self.num_nodes = description['sweeper_params']['num_nodes'] 

35 self.dt = description['level_params']['dt'] 

36 self.SDC, self.Ksdc, self.picard, self.Kpicard = self.stability_data() 

37 self.cwd = cwd 

38 

39 def stability_data(self): 

40 """ 

41 Computes stability domain matrix for the Harmonic oscillator problem 

42 Returns: 

43 numpy.ndarray: domain_SDC 

44 numpy.ndarray: domain_Ksdc 

45 numpy.ndarray: domain_picard 

46 numpy.ndarray: domain_Kpicard 

47 """ 

48 S = step(description=self.description) 

49 # Define L to get access level params and functions 

50 L = S.levels[0] 

51 # Number of nodes 

52 num_nodes = L.sweep.coll.num_nodes 

53 # Time step 

54 dt = L.params.dt 

55 

56 # Define Collocation matrix to find for the stability function 

57 Q = L.sweep.coll.Qmat[1:, 1:] 

58 QQ = np.dot(Q, Q) 

59 Q_coll = np.block([[QQ, np.zeros([num_nodes, num_nodes])], [np.zeros([num_nodes, num_nodes]), Q]]) 

60 qQ = np.dot(L.sweep.coll.weights, Q) 

61 # Matrix with all entries 1 

62 ones = np.block([[np.ones(num_nodes), np.zeros(num_nodes)], [np.zeros(num_nodes), np.ones(num_nodes)]]) 

63 # Combine all of the weights into a single matrix 

64 q_mat = np.block( 

65 [ 

66 [dt**2 * qQ, np.zeros(num_nodes)], 

67 [np.zeros(num_nodes), dt * L.sweep.coll.weights], 

68 ] 

69 ) 

70 # Zeros matrices to store the values for the stability region values 

71 domain_SDC = np.zeros((self.kappa_iter, self.mu_iter), dtype="complex") 

72 domain_picard = np.zeros((self.kappa_iter, self.mu_iter)) 

73 domain_Ksdc = np.zeros((self.kappa_iter, self.mu_iter)) 

74 domain_Kpicard = np.zeros((self.kappa_iter, self.mu_iter)) 

75 # Loop over the different values of the kappa and mu values 

76 for i in range(0, self.kappa_iter): 

77 for j in range(0, self.mu_iter): 

78 k = self.lambda_kappa[i] 

79 mu = self.lambda_mu[j] 

80 # Build right hand side matrix function for the harmonic oscillator problem 

81 F = np.block( 

82 [ 

83 [-k * np.eye(num_nodes), -mu * np.eye(num_nodes)], 

84 [-k * np.eye(num_nodes), -mu * np.eye(num_nodes)], 

85 ] 

86 ) 

87 

88 if self.K_iter != 0: 

89 # num iteration is not equal to zero then do SDC and Picard iteration 

90 lambdas = [k, mu] 

91 SDC_mat_sweep, Ksdc_eigval = L.sweep.get_scalar_problems_manysweep_mats( 

92 nsweeps=self.K_iter, lambdas=lambdas 

93 ) 

94 # If picard_mats_sweep=True then do also Picard iteration 

95 if L.sweep.params.picard_mats_sweep: 

96 ( 

97 Picard_mat_sweep, 

98 Kpicard_eigval, 

99 ) = L.sweep.get_scalar_problems_picardsweep_mats(nsweeps=self.K_iter, lambdas=lambdas) 

100 else: 

101 ProblemError("Picard iteration is not enabled. Set 'picard_mats_sweep' to True to enable.") 

102 domain_Ksdc[i, j] = Ksdc_eigval 

103 if L.sweep.params.picard_mats_sweep: 

104 domain_Kpicard[i, j] = Kpicard_eigval 

105 

106 else: 

107 # Otherwise Collocation problem 

108 SDC_mat_sweep = np.linalg.inv(np.eye(2 * num_nodes) - dt * np.dot(Q_coll, F)) 

109 # Collation update for both Picard and SDC iterations 

110 if L.sweep.params.do_coll_update: 

111 FSDC = np.dot(F, SDC_mat_sweep) 

112 Rsdc_mat = np.array([[1.0, dt], [0, 1.0]]) + np.dot(q_mat, FSDC) @ ones.T 

113 stab_func, v = np.linalg.eig(Rsdc_mat) 

114 

115 if L.sweep.params.picard_mats_sweep: 

116 FPicard = np.dot(F, Picard_mat_sweep) 

117 Rpicard_mat = np.array([[1.0, dt], [0, 1.0]]) + np.dot(q_mat, FPicard) @ ones.T 

118 stab_func_picard, v = np.linalg.eig(Rpicard_mat) 

119 else: 

120 raise ProblemError("Collocation update step works only when 'do_coll_update' is set to True.") 

121 # Find and store spectral radius 

122 domain_SDC[i, j] = np.max(np.abs(stab_func)) 

123 if L.sweep.params.picard_mats_sweep: 

124 domain_picard[i, j] = np.max(np.abs(stab_func_picard)) 

125 

126 return ( 

127 dt * domain_SDC.real, 

128 dt * domain_Ksdc.real, 

129 dt * domain_picard.real, 

130 dt * domain_Kpicard.real, 

131 ) 

132 

133 def stability_function_RKN(self, k, mu, dt): 

134 """ 

135 Stability function of RKN method 

136 

137 Returns: 

138 float: maximum absolute values of eigvales 

139 """ 

140 A = np.array([[0, 0, 0, 0], [0.5, 0, 0, 0], [0, 0.5, 0, 0], [0, 0, 1, 0]]) 

141 B = np.array([[0, 0, 0, 0], [0.125, 0, 0, 0], [0.125, 0, 0, 0], [0, 0, 0.5, 0]]) 

142 c = np.array([0, 0.5, 0.5, 1]) 

143 b = np.array([1 / 6, 2 / 6, 2 / 6, 1 / 6]) 

144 bA = np.array([1 / 6, 1 / 6, 1 / 6, 0]) 

145 L = np.eye(4) + k * (dt**2) * B + mu * dt * A 

146 R = np.block([[-k * np.ones(4)], [-(k * c + mu * np.ones(4))]]) 

147 

148 K = np.linalg.inv(L) @ R.T 

149 C = np.block([[dt**2 * bA], [dt * b]]) 

150 Y = np.array([[1, dt], [0, 1]]) + C @ K 

151 eigval = np.linalg.eigvals(Y) 

152 

153 return np.max(np.abs(eigval)) 

154 

155 def stability_data_RKN(self): 

156 """ 

157 Compute and store values into a matrix 

158 

159 Returns: 

160 numpy.ndarray: stab_RKN 

161 """ 

162 stab_RKN = np.zeros([self.kappa_iter, self.mu_iter]) 

163 for ii, kk in enumerate(self.lambda_kappa): 

164 for jj, mm in enumerate(self.lambda_mu): 

165 stab_RKN[jj, ii] = self.stability_function_RKN(kk, mm, self.dt) 

166 

167 return stab_RKN 

168 

169 def plot_stability(self, region, title=""): # pragma: no cover 

170 """ 

171 Plotting runtine for moduli 

172 

173 Args: 

174 stabval (numpy.ndarray): moduli 

175 title: title for the plot 

176 """ 

177 set_fixed_plot_params() 

178 lam_k_max = np.amax(self.lambda_kappa) 

179 lam_mu_max = np.amax(self.lambda_mu) 

180 

181 plt.figure() 

182 levels = np.array([0.25, 0.5, 0.75, 0.9, 1.0, 1.1]) 

183 

184 CS1 = plt.contour(self.lambda_kappa, self.lambda_mu, region.T, levels, colors='k', linestyles="dashed") 

185 # CS2 = plt.contour(self.lambda_k, self.lambda_mu, np.absolute(region.T), [1.0], colors='r') 

186 

187 plt.clabel(CS1, inline=True, fmt="%3.2f") 

188 

189 plt.gca().set_xticks(np.arange(0, int(lam_k_max) + 3, 3)) 

190 plt.gca().set_yticks(np.arange(0, int(lam_mu_max) + 3, 3)) 

191 plt.gca().tick_params(axis="both", which="both") 

192 plt.xlim([0.0, lam_k_max]) 

193 plt.ylim([0.0, lam_mu_max]) 

194 

195 plt.xlabel(r"$\Delta t\cdot \kappa$", labelpad=0.0) 

196 plt.ylabel(r"$\Delta t\cdot \mu$", labelpad=0.0) 

197 if self.RKN: 

198 plt.title(f"{title}") 

199 if self.radius: 

200 plt.title("{} $M={}$".format(title, self.num_nodes)) 

201 else: 

202 plt.title(r"{} $M={},\ K={}$".format(title, self.num_nodes, self.K_iter)) 

203 plt.tight_layout() 

204 plt.savefig(self.cwd + "data/M={}_K={}_redion_{}.pdf".format(self.num_nodes, self.K_iter, title)) 

205 

206 def run_SDC_stability(self): # pragma: no cover 

207 self.RKN = False 

208 self.radius = False 

209 self.plot_stability(self.SDC, title="SDC stability region") 

210 

211 def run_Picard_stability(self): # pragma: no cover 

212 self.RKN = False 

213 self.radius = False 

214 self.plot_stability(self.picard, title="Picard stability region") 

215 

216 def run_Ksdc(self): # pragma: no cover 

217 self.radius = True 

218 self.plot_stability(self.Ksdc, title="$K_{sdc}$ spectral radius") 

219 

220 def run_Kpicard(self): # pragma: no cover 

221 self.radius = True 

222 self.plot_stability(self.Kpicard, title="$K_{picard}$ spectral radius") 

223 

224 def run_RKN_stability(self): # pragma: no cover 

225 self.RKN = True 

226 self.radius = False 

227 region_RKN = self.stability_data_RKN() 

228 self.plot_stability(region_RKN.T, title='RKN-4 stability region') 

229 

230 

231def check_points_and_interval( 

232 description, helper_params, point, compute_interval=False, check_stability_point=False, Picard=False 

233): 

234 # Storage for stability interval and stability check 

235 interval_data = [] 

236 points_data = [] 

237 

238 # Loop through different numbers of nodes and maximum iterations 

239 for quad_type in helper_params['quad_type_list']: 

240 for num_nodes in helper_params['num_nodes_list']: 

241 for max_iter in helper_params['max_iter_list']: 

242 # Update simulation parameters 

243 description['sweeper_params']['num_nodes'] = num_nodes 

244 description['sweeper_params']['quad_type'] = quad_type 

245 description['step_params']['maxiter'] = max_iter 

246 

247 # Create Stability_implementation instance for stability check 

248 

249 stab_model = StabilityImplementation( 

250 description, kappa_max=point[0], mu_max=point[1], Num_iter=helper_params['Num_iter'] 

251 ) 

252 if compute_interval: 

253 # Extract the values where SDC is less than or equal to 1 

254 if Picard: 

255 mask = stab_model.picard <= 1 + 1e-14 

256 else: 

257 mask = stab_model.SDC <= 1.0 

258 for ii in range(len(mask)): 

259 if mask[ii]: 

260 kappa_max_interval = stab_model.lambda_kappa[ii] 

261 else: 

262 break 

263 

264 # Add row to the interval data 

265 interval_data.append([quad_type, num_nodes, max_iter, kappa_max_interval]) 

266 

267 if check_stability_point: 

268 # Check stability and print results 

269 if stab_model.SDC[-1, -1] <= 1: 

270 stability_result = "Stable" 

271 else: 

272 stability_result = "Unstable. Increase M or K" 

273 

274 # Add row to the results data 

275 points_data.append( 

276 [quad_type, num_nodes, max_iter, point, stability_result, stab_model.SDC[-1, -1]] 

277 ) 

278 if compute_interval: 

279 return interval_data 

280 else: 

281 return points_data 

282 

283 

284def compute_and_generate_table( 

285 description, 

286 helper_params, 

287 point, 

288 compute_interval=False, 

289 save_interval_file=False, 

290 interval_filename='./data/stab_interval.txt', 

291 check_stability_point=False, 

292 save_points_table=False, 

293 points_table_filename='./data/point_table.txt', 

294 quadrature_list=('GAUSS', 'LOBATTO'), 

295 Picard=False, 

296): # pragma: no cover 

297 from tabulate import tabulate 

298 

299 if compute_interval: 

300 interval_data = check_points_and_interval( 

301 description, helper_params, point, compute_interval=compute_interval, Picard=Picard 

302 ) 

303 else: 

304 points_data = check_points_and_interval( 

305 description, helper_params, point, check_stability_point=check_stability_point 

306 ) 

307 

308 # Define column names for interval data 

309 interval_headers = ["Quad Type", "Num Nodes", "Max Iter", 'kappa_max'] 

310 

311 # Define column names for results data 

312 points_headers = ["Quad Type", "Num Nodes", "Max Iter", "(kappa, mu)", "Stability", "Spectral Radius"] 

313 # Print or save the tables using tabulate 

314 if save_interval_file and compute_interval: 

315 interval_table_str = tabulate(interval_data, headers=interval_headers, tablefmt="grid") 

316 with open(interval_filename, 'w') as file: 

317 file.write(interval_table_str) 

318 print(f"Stability Interval Table saved to {interval_filename}") 

319 

320 if save_points_table and check_stability_point: 

321 points_table_str = tabulate(points_data, headers=points_headers, tablefmt="grid") 

322 with open(points_table_filename, 'w') as file: 

323 file.write(points_table_str) 

324 print(f"Stability Results Table saved to {points_table_filename}") 

325 

326 if compute_interval: 

327 if Picard: 

328 print("Picard stability Interval Table:") 

329 else: 

330 print("SDC stability Interval Table:") 

331 print(tabulate(interval_data, headers=interval_headers, tablefmt="grid")) 

332 

333 if check_stability_point: 

334 print("\nStability Results Table:") 

335 print(tabulate(points_data, headers=points_headers, tablefmt="grid"))