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

## 107 statements

, created at 2024-09-20 17:10 +0000

1import numpy as np

2import matplotlib.pyplot as plt

3from pySDC.core.errors import ProblemError

4from pySDC.core.step import Step

6from pySDC.projects.Second_orderSDC.plot_helper import set_fixed_plot_params

9class StabilityImplementation:

10 """

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

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

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

22 """

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)

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

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

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 )

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

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)

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))

126 return (

127 dt * domain_SDC.real,

128 dt * domain_Ksdc.real,

129 dt * domain_picard.real,

130 dt * domain_Kpicard.real,

131 )

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

134 """

135 Stability function of RKN method

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))]])

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)

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

155 def stability_data_RKN(self):

156 """

157 Compute and store values into a matrix

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)

167 return stab_RKN

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

170 """

171 Plotting runtine for moduli

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)

181 plt.figure()

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

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')

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

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])

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}")

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))

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

207 self.RKN = False

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

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

212 self.RKN = False

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

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

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

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

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

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

225 self.RKN = True

227 region_RKN = self.stability_data_RKN()

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

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 = []

238 # Loop through different numbers of nodes and maximum iterations

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

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

247 # Create Stability_implementation instance for stability check

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

260 kappa_max_interval = stab_model.lambda_kappa[ii]

261 else:

262 break

264 # Add row to the interval data

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"

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

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',

295 Picard=False,

296): # pragma: no cover

297 from tabulate import tabulate

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 )

308 # Define column names for interval data

311 # Define column names for results data

313 # Print or save the tables using tabulate

314 if save_interval_file and compute_interval:

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

317 file.write(interval_table_str)

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

320 if save_points_table and check_stability_point:

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}")

326 if compute_interval:

327 if Picard:

328 print("Picard stability Interval Table:")

329 else:

330 print("SDC stability Interval Table:")