Coverage for bim2sim/plugins/PluginOpenFOAM/bim2sim_openfoam/utils/postProcessingPlots.py: 0%
358 statements
« prev ^ index » next coverage.py v7.10.7, created at 2025-10-01 10:24 +0000
1import os
2from pathlib import Path
3import re
4import pandas as pd
5import numpy as np
6import matplotlib.pyplot as plt
8"""
9A collection of post-processing python scripts for plotting different
10parameters from OpenFOAM simulations.
11"""
12from itertools import cycle
15# Function to get the next unused color
16def get_next_unused_color(ax):
17 # Get current color cycle
18 color_cycle = cycle(plt.rcParams['axes.prop_cycle'].by_key()['color'])
20 # Get colors already used in the Axes
21 used_colors = [line.get_color() for line in ax.get_lines()]
23 # Find the next unused color
24 for color in color_cycle:
25 if color not in used_colors:
26 return color
28 # If all colors are used, return the first color in the cycle (fallback)
29 return next(color_cycle)
31def convergencePlot(of_directory: str):
32 """
33 Creates a plot of velocity (U) and temperature (T, directly proportional
34 to h) residuals over the simulated iteration steps.
35 """
36 Ux_Init = []
37 Ux_Final = []
38 Uy_Init = []
39 Uy_Final = []
40 Uz_Init = []
41 Uz_Final = []
42 h_Init = []
43 h_Final = []
44 values = [Ux_Init, Ux_Final, Uy_Init, Uy_Final, Uz_Init, Uz_Final,
45 h_Init, h_Final]
46 # Collect directories from all simulation phases and sort them in
47 # ascending time steps
48 time_dirs = os.listdir(of_directory / 'postProcessing/solverInfo/')
49 time_dirs.sort()
50 time_dirs.sort(key=len) # yes, these are both necessary
51 for timestep in time_dirs:
52 with open(of_directory / 'postProcessing/solverInfo/' /
53 timestep / 'solverInfo.dat', 'r') as f:
54 lines = f.readlines()
55 track_U = track_h = False
56 indices = lines[1].split('\t')
57 if 'U_solver' in lines[1]:
58 track_U = True
59 ind_u = indices.index('U_solver ')
60 if 'h_solver' in lines[1]:
61 track_h = True
62 ind_h = indices.index('h_solver ')
63 for line in lines[2:]:
64 line = line.split('\t')
65 if track_U and len(line) > 5:
66 [uxinit, uxfinal] = line[ind_u+1:ind_u+3]
67 [uyinit, uyfinal] = line[ind_u+4:ind_u+6]
68 [uzinit, uzfinal] = line[ind_u+7:ind_u+9]
69 else:
70 uxinit = uxfinal = uyinit = uyfinal = uzinit = uzfinal = 0
71 Ux_Init.append(float(uxinit))
72 Ux_Final.append(float(uxfinal))
73 Uy_Init.append(float(uyinit))
74 Uy_Final.append(float(uyfinal))
75 Uz_Init.append(float(uzinit))
76 Uz_Final.append(float(uzfinal))
77 if track_h and len(line) > 5:
78 [hinit, hfinal] = line[ind_h+1:ind_h+3]
79 else:
80 hinit = hfinal = 0
81 h_Init.append(float(hinit))
82 h_Final.append(float(hfinal))
84 legend = []
85 if track_U:
86 legend.extend(['Ux_Init', 'Ux_Final', 'Uy_Init', 'Uy_Final', 'Uz_Init',
87 'Uz_Final'])
88 if track_h:
89 legend.extend(['h_Init', 'h_Final'])
90 for value in values:
91 plt.plot(value)
92 plt.ylabel('Residual')
93 plt.xlabel('Iteration')
94 plt.legend(legend)
95 plt.yscale("log")
96 plt.title('Residuals')
97 plt.savefig(of_directory / 'Residuals.png')
98 # plt.show()
101def convergencePlot2(of_directory: str):
102 """
103 Creates a plot of velocity (U) and temperature (T, directly proportional
104 to h) residuals over the simulated iteration steps and detects the
105 stabilization point where residuals converge with minimal variations.
106 """
107 Ux_Init = []
108 Ux_Final = []
109 Uy_Init = []
110 Uy_Final = []
111 Uz_Init = []
112 Uz_Final = []
113 h_Init = []
114 h_Final = []
115 values = [Ux_Init, Ux_Final, Uy_Init, Uy_Final, Uz_Init, Uz_Final,
116 h_Init, h_Final]
118 # Collect directories from all simulation phases and sort them
119 time_dirs = os.listdir(of_directory / 'postProcessing/solverInfo/')
120 time_dirs.sort()
121 time_dirs.sort(key=len) # both sorts needed for correct order
123 # Parse residual data
124 for timestep in time_dirs:
125 with open(of_directory / 'postProcessing/solverInfo/' /
126 timestep / 'solverInfo.dat', 'r') as f:
127 lines = f.readlines()
128 track_U = track_h = False
129 indices = lines[1].split('\t')
130 if 'U_solver' in lines[1]:
131 track_U = True
132 ind_u = indices.index('U_solver ')
133 if 'h_solver' in lines[1]:
134 track_h = True
135 ind_h = indices.index('h_solver ')
136 for line in lines[2:]:
137 line = line.split('\t')
138 if track_U and len(line) > 5:
139 [uxinit, uxfinal] = line[ind_u + 1:ind_u + 3]
140 [uyinit, uyfinal] = line[ind_u + 4:ind_u + 6]
141 [uzinit, uzfinal] = line[ind_u + 7:ind_u + 9]
142 else:
143 uxinit = uxfinal = uyinit = uyfinal = uzinit = uzfinal = 0
144 Ux_Init.append(float(uxinit))
145 Ux_Final.append(float(uxfinal))
146 Uy_Init.append(float(uyinit))
147 Uy_Final.append(float(uyfinal))
148 Uz_Init.append(float(uzinit))
149 Uz_Final.append(float(uzfinal))
150 if track_h and len(line) > 5:
151 [hinit, hfinal] = line[ind_h + 1:ind_h + 3]
152 else:
153 hinit = hfinal = 0
154 h_Init.append(float(hinit))
155 h_Final.append(float(hfinal))
157 # Helper functions
158 def smooth_data(data, window_size=500):
159 """Apply a moving average to smooth out small peaks in the data."""
160 return pd.DataFrame(data).rolling(window=window_size,
161 center=True).median().values
163 def find_stabilization_iteration(smoothed_residuals, min_tolerance=1e-6,
164 stop_tolerance=1e-8,
165 consecutive=1000):
166 """Find the first iteration where residuals stabilize within a given tolerance."""
168 curr_abs = 1
169 final_iter = None
170 for start_iter in range(len(smoothed_residuals) - consecutive):
171 if start_iter < 1500:
172 continue
173 window = smoothed_residuals[start_iter:start_iter + consecutive]
174 this_abs = np.abs(np.diff(window, axis=0))
175 this_abs = pd.DataFrame(this_abs).median().values
176 if this_abs[0] > min_tolerance:
177 continue
178 if this_abs[0] < curr_abs and this_abs[0] > 0:
179 curr_abs = this_abs[0]
180 final_iter = start_iter
181 if this_abs[0] < stop_tolerance:
182 return final_iter
183 return final_iter # If stabilization is not found
185 # Apply smoothing and find stabilization points for each residual
186 stabilization_points = []
187 smoothed_values = [smooth_data(value) for value in values]
188 global_max = 0
189 for smoothed_value in smoothed_values:
190 stabilization_iteration = find_stabilization_iteration(smoothed_value)
191 if stabilization_iteration and stabilization_iteration > global_max:
192 global_max = stabilization_iteration
193 stabilization_points.append(stabilization_iteration)
195 # Prepare the legend
196 legend = []
197 if 'U_solver ' in indices: # Check if 'U_solver' is in the data
198 legend.extend(['Ux_Init', 'Ux_Final', 'Uy_Init', 'Uy_Final', 'Uz_Init',
199 'Uz_Final'])
200 if 'h_solver ' in indices: # Check if 'h_solver' is in the data
201 legend.extend(['h_Init', 'h_Final'])
203 # Check if the number of legend labels matches the number of values
204 if len(legend) != len(values):
205 print(
206 f"Warning: The number of legend labels ({len(legend)}) doesn't match the number of values ({len(values)}).")
207 print("Legend and values may not correspond correctly.")
208 # Adjust the legend size or values accordingly to match
210 # Plotting
211 plt.figure(figsize=(10, 6))
212 for idx, value in enumerate(values):
213 plt.plot(value, label=legend[idx] if idx < len(
214 legend) else f"Parameter {idx + 1}")
215 if stabilization_points[idx] == global_max:
216 plt.axvline(stabilization_points[idx], color='r', linestyle='--')
217 if stabilization_points[idx] is not None:
218 plt.plot(stabilization_points[idx], value[stabilization_points[
219 idx]], 'ro')
220 plt.text(stabilization_points[idx], value[stabilization_points[
221 idx]], f'(it={stabilization_points[idx]:d},'
222 f' v={value[stabilization_points[idx]]:.2e})',
223 ha='right', va='bottom', color='red', fontsize=10)
224 print('global max: ', global_max)
226 plt.ylabel('Residual')
227 plt.xlabel('Iteration')
228 plt.legend(legend)
229 plt.yscale("log")
230 plt.title('Residuals with Stabilization Points')
231 plt.savefig(of_directory / 'Residuals.png')
232 plt.show()
233 return global_max
235def MinMaxPlot(of_directory: str, min_max_dir='MinMax',
236 vol_field_dir='volFieldValue', max_dir='', min_dir='', name=''):
237 """
238 Creates a plot of minimal and maximal values of the velocity's magnitude
239 (U) and temperature (T) over the simulated iteration steps.
240 """
241 T_min = []
242 T_max = []
243 T_av = []
244 AoA_min = []
245 AoA_max = []
246 U_mag_min = []
247 U_mag_max = []
248 max_dirs = None
249 min_dirs = None
250 time_dirs = None
251 if not (max_dir and min_dir):
252 time_dirs = os.listdir(of_directory / f'postProcessing/{min_max_dir}/')
253 time_dirs.sort()
254 time_dirs.sort(key=len) # yes, these are both necessary
255 else:
256 max_dirs = os.listdir(of_directory / f'postProcessing/{max_dir}/')
257 min_dirs = os.listdir(of_directory / f'postProcessing/{min_dir}/')
258 min_dirs.sort()
259 min_dirs.sort(key=len) # yes, these are both necessary
260 max_dirs.sort()
261 max_dirs.sort(key=len) # yes, these are both necessary
262 mean_dirs = os.listdir(of_directory / f'postProcessing/{vol_field_dir}')
263 mean_dirs.sort()
264 mean_dirs.sort(key=len)
265 if time_dirs:
266 for timestep in time_dirs:
267 with open(of_directory / f'postProcessing/{min_max_dir}' /
268 timestep / 'fieldMinMax.dat', 'r') as f:
269 lines = f.readlines()
270 for i in range(11, len(lines)):
271 line = lines[i].split('\t')
272 if len(line) > 5 and 'T' in line[1]:
273 tmin = line[2]
274 tmax = line[5]
275 T_min.append(float(tmin))
276 T_max.append(float(tmax))
277 elif len(line) > 5 and 'mag(U)' in line[1]:
278 umin = line[2]
279 umax = line[5]
280 U_mag_min.append(float(umin))
281 U_mag_max.append(float(umax))
282 elif len(line) > 5 and 'AoA' in line[1]:
283 AoA_min.append(float(line[2]))
284 AoA_max.append(float(line[5]))
285 if min_dirs and max_dirs:
286 for timestep in min_dirs:
287 with open(of_directory / f'postProcessing/{min_dir}' /
288 timestep / 'volFieldValue.dat', 'r') as f3:
289 lines = f3.readlines()
290 for i in range(4, len(lines)):
291 line = lines[i].split('\t')
292 tmin = line[1]
293 T_min.append(float(tmin))
294 for timestep in max_dirs:
295 with open(of_directory / f'postProcessing/{max_dir}' /
296 timestep / 'volFieldValue.dat', 'r') as f4:
297 lines = f4.readlines()
298 for i in range(4, len(lines)):
299 line = lines[i].split('\t')
300 tmax = line[1]
301 T_max.append(float(tmax))
302 for timestep in mean_dirs:
303 with open(of_directory/ f'postProcessing/{vol_field_dir}' /
304 timestep / 'volFieldValue.dat', 'r') as f2:
305 lines = f2.readlines()
306 for i in range(4, len(lines)):
307 line = lines[i].split('\t')
308 avT = line[1]
309 T_av.append(float(avT))
310 plt.plot(T_min, linewidth=1)
311 plt.plot(T_max, linewidth=1)
312 plt.plot(T_av, linewidth=1)
313 plt.text(len(T_av), T_av[-1], f'T_mean_final: {T_av[-1]:.2f} K',
314 ha='right', va='bottom', fontsize=12)
315 plt.ylabel('T')
316 plt.xlabel('Iteration')
317 plt.legend(['T_min', 'T_max', 'T_av'])
318 plt.title(f'{name} Temperature min/max/average')
319 plt.savefig(of_directory / f'{name}minmaxavT.png')
320 plt.show()
321 plt.close()
322 plt.plot(AoA_min, linewidth=1)
323 plt.plot(AoA_max, linewidth=1)
324 plt.ylabel('AoA')
325 plt.xlabel('Iteration')
326 plt.legend(['AoA_min', 'AoA_max'])
327 plt.title(f'{name} AoA min/max')
328 plt.savefig(of_directory / f'{name}_AoA.png')
329 plt.show()
330 plt.close()
331 if time_dirs:
332 plt.plot(U_mag_min)
333 plt.plot(U_mag_max)
334 plt.ylabel('|U|')
335 plt.xlabel('Iteration')
336 plt.legend(['|U|_min', '|U|_max'])
337 plt.title(f'{name} mag(U) min/max')
338 plt.savefig(of_directory / f'{name}minmax_U_.png')
339 # plt.show()
342def analyze_execution_times(of_directory, target_iterations=[1000, 'final']):
343 """
344 Extracts execution times from an OpenFOAM log file, prints times for specific iterations,
345 and plots all execution times over the course of the iterations with target times annotated.
347 Parameters:
348 - of_directory (str or Path): Directory containing the OpenFOAM log file.
349 - target_iterations (list): List of iteration numbers to print and annotate (e.g., [1000, 'final']).
350 """
351 # Initialize lists to store all iterations and their execution times
352 all_iterations = []
353 all_execution_times = []
355 # Dictionary to store the target iteration times for printing
356 target_times = {}
357 last_time = None
359 # Regular expressions to match lines
360 iteration_re = re.compile(r"Time = (\d+)")
361 execution_re = re.compile(r"ExecutionTime = ([\d.]+)")
363 # Parse the log file
364 with open(of_directory / 'logSimulation.compress', 'r') as file:
365 current_iteration = None
366 iteration_done = True
367 for line in file:
368 # Match iteration number
369 iteration_match = iteration_re.search(line)
370 if iteration_match and iteration_done:
371 current_iteration = int(iteration_match.group(1))
372 iteration_done = False
374 # Match execution time if we have an iteration
375 execution_match = execution_re.search(line)
376 if execution_match and current_iteration is not None:
377 execution_time = float(execution_match.group(1))
378 last_time = (current_iteration, execution_time)
380 # Store all iterations and execution times for plotting
381 all_iterations.append(current_iteration)
382 all_execution_times.append(execution_time)
384 # Store target times for specified iterations
385 if current_iteration in target_iterations:
386 target_times[current_iteration] = execution_time
387 iteration_done = True
389 # Store the last recorded execution time if 'final' is specified in target_iterations
390 if 'final' in target_iterations and last_time:
391 target_times[last_time[0]] = last_time[1] # Use the final iteration number as the key
393 # Print the execution times for the specified target iterations
394 print("Execution times for specified iterations:")
395 for iteration, exec_time in target_times.items():
396 print(f"Iteration {iteration}: {exec_time} seconds")
398 # Plot all execution times over the iterations
399 plt.figure(figsize=(10, 6))
400 plt.plot(all_iterations, all_execution_times, 'b-', label="Execution Time")
401 plt.xlabel("Iteration")
402 plt.ylabel("Execution Time (s)")
403 plt.title("Execution Time Over Iterations")
404 plt.legend()
405 plt.grid(True)
407 # Annotate target times on the plot
408 for iteration, exec_time in target_times.items():
409 plt.text(iteration, exec_time, f'{exec_time:.2f}s',
410 ha='right', va='bottom', color='red', fontsize=10)
412 # Save the plot
413 plt.savefig(of_directory / 'iteration_time.png')
415 # plt.show()
416 plt.close()
417 return plt
419def add_simulation_times(fig, of_directory, name='', number=0):
420 all_iterations = []
421 all_execution_times = []
423 # Regular expressions to match lines
424 iteration_re = re.compile(r"Time = (\d+)")
425 execution_re = re.compile(r"ExecutionTime = ([\d.]+)")
427 # Parse the log file
428 with open(of_directory / 'logSimulation.compress', 'r') as file:
429 current_iteration = None
430 iteration_done = True
431 for line in file:
432 # Match iteration number
433 iteration_match = iteration_re.search(line)
434 if iteration_match and iteration_done:
435 current_iteration = int(iteration_match.group(1))
436 iteration_done = False
438 # Match execution time if we have an iteration
439 execution_match = execution_re.search(line)
440 if execution_match and current_iteration is not None:
441 execution_time = float(execution_match.group(1))
442 last_time = (current_iteration, execution_time)
444 # Store all iterations and execution times for plotting
445 all_iterations.append(current_iteration)
446 all_execution_times.append(execution_time)
448 # Store target times for specified iterations
449 iteration_done = True
451 # Store the last recorded execution time if 'final' is specified in target_iterations
453 # Plot all execution times over the iterations
454 if not fig:
455 fig, ax = plt.subplots(figsize=(9, 6))
456 ax.set_xlabel("Iteration")
457 ax.set_ylabel("Execution Time (s)")
458 ax.set_title("Execution Time Over Iterations")
459 ax.grid(True)
460 else:
461 ax = fig.axes[0]
462 this_color = get_next_unused_color(ax)
463 ax.plot(all_iterations, all_execution_times, color=this_color,
464 label=f"Execution Time {name}")
465 # line = ax.lines[0] # Access the first plotted line
466 # line.set_label(f"Execution Time {name}")
467 ax.legend(loc='upper right')
469 # Annotate target times on the plot
470 ax.text(all_iterations[-1], all_execution_times[-1], f'{all_execution_times[-1]:.2f}s',
471 ha='right', va='bottom', fontsize=18, color=this_color)
472 for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
473 ax.get_xticklabels() + ax.get_yticklabels() + ax.get_legend().get_texts()):
474 item.set_fontsize(14)
475 # Save the plot
476 fig.savefig(of_directory / f'iteration_time_V{number}.png')
478 #fig.show()
479 return fig
481if __name__ == '__main__':
482 # this_path = Path(r'C:\Users\richter\sciebo\03-Paperdrafts\00-Promotion\05'
483 # r'-SIM-Data\02-CFD\diss_fv_100\OpenFOAM')
484 # global_conv_iter = convergencePlot2(this_path)
485 # MinMaxPlot(this_path)
486 # if global_conv_iter:
487 # analyze_execution_times(this_path, target_iterations=[global_conv_iter, 'final'])
488 # else:
489 # analyze_execution_times(this_path, target_iterations=[1000, 'final'])
490 #
491 directory = Path(r'C:\Users\richter\sciebo\03-Paperdrafts\BS2025\SimData\seminar_1OG_south')
492 # Iterate through directories that start with 'diss_'
494 fig_temp = None
495 counter=0
496 for diss_dir in directory.glob('*'):
497 # Check if "OpenFOAM" subdirectory exists within the current directory
498 # openfoam_dir = diss_dir / 'OpenFOAM'
499 openfoam_dir = diss_dir /'export' / 'OpenFOAM'
500 if openfoam_dir.is_dir():
501 print(openfoam_dir)
502 MinMaxPlot(openfoam_dir,
503 name=str(diss_dir.name + '_MinMax'))
504 try:
505 global_conv_iter = convergencePlot2(openfoam_dir)
506 MinMaxPlot(openfoam_dir, min_max_dir=None,
507 vol_field_dir='evaluate_air_volume_average',
508 max_dir='evaluate_air_volume_max',
509 min_dir='evaluate_air_volume_min',
510 name=str(diss_dir.name+'_volume'))
511 try:
512 MinMaxPlot(openfoam_dir, min_max_dir=None,
513 vol_field_dir='person_Person0_average',
514 max_dir='person_Person0_max',
515 min_dir='person_Person0_min',
516 name=str(diss_dir.name+'_Person0')
517 )
518 except:
519 print('No person found!')
520 if global_conv_iter:
521 analyze_execution_times(openfoam_dir,
522 target_iterations=[global_conv_iter,
523 'final'])
524 else:
525 analyze_execution_times(openfoam_dir,
526 target_iterations=[1000, 'final'])
527 fig_temp = add_simulation_times(fig_temp, openfoam_dir,
528 name=diss_dir.name.replace(
529 'diss_', ''),
530 number=counter)
531 plt.close('all')
532 counter+=1
533 except:
534 print(f"failed plot for {diss_dir}")
535 fig_temp.show()