Coverage for pySDC / projects / RayleighBenard / analysis_scripts / plot_Nu.py: 100%

1import pickle 

2import matplotlib.pyplot as plt 

3import numpy as np 

4from scipy import integrate 

5from pySDC.projects.RayleighBenard.analysis_scripts.process_RBC3D_data import get_pySDC_data 

6from pySDC.projects.RayleighBenard.analysis_scripts.plotting_utils import figsize, savefig 

7 

8 

9def interpolate_NuV_to_reference_times(data, reference_data, order=12): 

10 from qmat.lagrange import getSparseInterpolationMatrix 

11 

12 t_in = np.array(data['t']) 

13 t_out = np.array([me for me in reference_data['t'] if me <= max(t_in)]) 

14 

15 order = min([order, len(t_in), len(t_out)]) 

16 interpolation_matrix = getSparseInterpolationMatrix(t_in, t_out, order=order) 

17 return interpolation_matrix @ t_in, interpolation_matrix @ data['Nu']['V'] 

18 

19 

20def plot_Nu(res, dts, config_name, ref, ax, title, converged_from=0, **plotting_args): # pragma: no cover 

21 ax.plot(ref['t'], ref['Nu']['V'], color='black', ls='--') 

22 ax.set_title(title) 

23 Nu_ref = np.array(ref['Nu']['V']) 

24 

25 for dt in dts: 

26 data = get_pySDC_data(res=res, dt=dt, config_name=config_name) 

27 t_i, Nu_i = interpolate_NuV_to_reference_times(data, ref) 

28 

29 ax.plot(data['t'], data['Nu']['V'], **{'label': rf'$\Delta t={{{dt}}}$', **plotting_args}) 

30 

31 error = np.abs(Nu_ref[: Nu_i.shape[0]] - Nu_i) / np.abs(Nu_ref[: Nu_i.shape[0]]) 

32 

33 # compute mean Nu 

34 mask = np.logical_and(t_i >= converged_from, t_i <= np.inf) 

35 Nu_mean = np.mean(Nu_i[mask]) 

36 Nu_std = np.std(Nu_i[mask]) 

37 

38 last_line = ax.get_lines()[-1] 

39 if any(error > 1e-2): 

40 deviates = min(t_i[error > 1e-2]) 

41 ax.axvline(deviates, color=last_line.get_color(), ls=':') 

42 print(f'{title} dt={dt:.4f} Nu={Nu_mean:.3f}+={Nu_std:.3f}, deviates more than 1% from t={deviates:.2f}') 

43 else: 

44 print(f'{title} dt={dt:.4f} Nu={Nu_mean:.3f}+={Nu_std:.3f}') 

45 ax.legend(frameon=True, loc='upper left') 

46 

47 

48def plot_Nu_over_time_Ra1e5(): # pragma: no cover 

49 Nu_fig, Nu_axs = plt.subplots(4, 1, sharex=True, sharey=True, figsize=figsize(scale=1, ratio=1.4)) 

50 

51 res = 32 

52 converged_from = 100 

53 

54 ref_data = get_pySDC_data(res=res, dt=0.01, config_name='RBC3DG4R4SDC44Ra1e5') 

55 

56 plot_Nu(res, [0.06, 0.04, 0.02], 'RBC3DG4R4SDC44Ra1e5', ref_data, Nu_axs[0], 'SDC44', converged_from) 

57 plot_Nu(res, [0.06, 0.05, 0.02, 0.01], 'RBC3DG4R4SDC23Ra1e5', ref_data, Nu_axs[1], 'SDC23', converged_from) 

58 plot_Nu(res, [0.05, 0.04, 0.02, 0.01, 0.005], 'RBC3DG4R4RKRa1e5', ref_data, Nu_axs[2], 'RK443', converged_from) 

59 plot_Nu(res, [0.02, 0.01, 0.005], 'RBC3DG4R4EulerRa1e5', ref_data, Nu_axs[3], 'RK111', converged_from) 

60 

61 Nu_axs[-1].set_xlabel('$t$') 

62 Nu_axs[-1].set_ylabel('$Nu$') 

63 

64 Nu_fig.tight_layout() 

65 Nu_fig.savefig('./plots/Nu_over_time_Ra1e5.pdf', bbox_inches='tight') 

66 

67 

68def plot_Nu_over_time_Ra1e6(): # pragma: no cover 

69 Nu_fig, Nu_axs = plt.subplots(4, 1, sharex=True, sharey=True, figsize=figsize(scale=1, ratio=1.4)) 

70 

71 res = 64 

72 converged_from = 25 

73 

74 ref_data = get_pySDC_data(res=res, dt=0.002, config_name='RBC3DG4R4SDC34Ra1e6') 

75 

76 plot_Nu(res, [0.02, 0.01], 'RBC3DG4R4SDC44Ra1e6', ref_data, Nu_axs[0], 'SDC44', converged_from) 

77 plot_Nu(res, [0.01, 0.005, 0.002], 'RBC3DG4R4SDC23Ra1e6', ref_data, Nu_axs[1], 'SDC23', converged_from) 

78 plot_Nu(res, [0.01, 0.005, 0.002], 'RBC3DG4R4RKRa1e6', ref_data, Nu_axs[2], 'RK443', converged_from) 

79 plot_Nu(res, [0.005, 0.002], 'RBC3DG4R4EulerRa1e6', ref_data, Nu_axs[3], 'RK111', converged_from) 

80 

81 Nu_axs[-1].set_xlabel('$t$') 

82 Nu_axs[-1].set_ylabel('$Nu$') 

83 

84 Nu_fig.tight_layout() 

85 Nu_fig.savefig('./plots/Nu_over_time_Ra1e6.pdf', bbox_inches='tight') 

86 

87 

88def plot_Nusselt_Ra1e5_same_dt(): # pragma: no cover 

89 fig, ax = plt.subplots(figsize=figsize(scale=0.9, ratio=0.45)) 

90 

91 res = 32 

92 converged_from = 100 

93 

94 ref_data = get_pySDC_data(res=res, dt=0.01, config_name='RBC3DG4R4SDC44Ra1e5') 

95 

96 plot_Nu(res, [0.02], 'RBC3DG4R4SDC44Ra1e5', ref_data, ax, 'SDC44', converged_from) 

97 plot_Nu(res, [0.02], 'RBC3DG4R4SDC23Ra1e5', ref_data, ax, 'SDC23', converged_from) 

98 plot_Nu(res, [0.02], 'RBC3DG4R4RKRa1e5', ref_data, ax, 'RK443', converged_from) 

99 plot_Nu(res, [0.02], 'RBC3DG4R4EulerRa1e5', ref_data, ax, 'RK111', converged_from) 

100 

101 ax.set_xlabel('$t$') 

102 ax.set_ylabel('$Nu$') 

103 ax.set_title(r'$Ra=10^5 \quad \Delta t=0.02$') 

104 

105 lines = ax.get_lines() 

106 for line, label in zip(lines[1::3], ['SDC44', 'SDC23', 'RK443', 'RK111'], strict=True): 

107 line.set_label(label) 

108 ax.legend(frameon=False) 

109 

110 fig.tight_layout() 

111 savefig(fig, 'NuRa1e5SameDt', pad_inches=0.1) 

112 

113 

114if __name__ == '__main__': 

115 

116 # plot_Nu_over_time_Ra1e5() 

117 # plot_Nu_over_time_Ra1e6() 

118 plot_Nusselt_Ra1e5_same_dt() 

119 

120 plt.show()