Coverage for pySDC/projects/Second_orderSDC/stability_simulation.py: 90%
107 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +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}")
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))
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")
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")
216 def run_Ksdc(self): # pragma: no cover
217 self.radius = True
218 self.plot_stability(self.Ksdc, title="$K_{sdc}$ spectral radius")
220 def run_Kpicard(self): # pragma: no cover
221 self.radius = True
222 self.plot_stability(self.Kpicard, title="$K_{picard}$ spectral radius")
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')
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
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
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
258 for ii in range(len(mask)):
259 if mask[ii]:
260 kappa_max_interval = stab_model.lambda_kappa[ii]
261 else:
262 break
264 # Add row to the interval data
265 interval_data.append([quad_type, num_nodes, max_iter, kappa_max_interval])
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',
294 quadrature_list=('GAUSS', 'LOBATTO'),
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
309 interval_headers = ["Quad Type", "Num Nodes", "Max Iter", 'kappa_max']
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}")
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}")
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"))
333 if check_stability_point:
334 print("\nStability Results Table:")
335 print(tabulate(points_data, headers=points_headers, tablefmt="grid"))