Coverage for pySDC / projects / RayleighBenard / analysis_scripts / process_RBC3D_data.py: 95%

1from pySDC.projects.RayleighBenard.RBC3D_configs import get_config 

2from pySDC.projects.RayleighBenard.run_experiment import parse_args 

3from pySDC.helpers.fieldsIO import FieldsIO 

4import matplotlib.pyplot as plt 

5from tqdm import tqdm 

6from mpi4py import MPI 

7import numpy as np 

8import pickle 

9import os 

10 

11 

12def process_RBC3D_data(base_path='./data/processed', plot=True, args=None, config=None): 

13 # prepare problem instance 

14 args = args if args else parse_args() 

15 comm = MPI.COMM_WORLD 

16 args['procs'] = [1, 1, comm.size] 

17 config = config if config else get_config(args) 

18 desc = config.get_description(**args) 

19 P = desc['problem_class']( 

20 **{ 

21 **desc['problem_params'], 

22 'spectral_space': False, 

23 'comm': comm, 

24 'Dirichlet_recombination': False, 

25 'left_preconditioner': False, 

26 } 

27 ) 

28 P.setUpFieldsIO() 

29 zInt = P.axes[-1].get_integration_weights() 

30 xp = P.xp 

31 

32 # prepare paths 

33 os.makedirs(base_path, exist_ok=True) 

34 fname = config.get_file_name() 

35 fname_trim = fname[fname.index('RBC3D') : fname.index('.pySDC')] 

36 path = f'{base_path}/{fname_trim}.pickle' 

37 

38 # open simulation data 

39 data = FieldsIO.fromFile(fname) 

40 

41 # prepare arrays to store data in 

42 Nu = { 

43 'V': [], 

44 'b': [], 

45 't': [], 

46 'thermal': [], 

47 'kinetic': [], 

48 } 

49 t = [] 

50 profiles = {key: [] for key in ['T', 'u', 'v', 'w']} 

51 rms_profiles = {key: [] for key in profiles.keys()} 

52 spectrum = [] 

53 spectrum_all = [] 

54 

55 # try to load time averaged values 

56 u_mean_profile = P.u_exact() 

57 if os.path.isfile(path): 

58 with open(path, 'rb') as file: 

59 avg_data = pickle.load(file) 

60 if comm.rank == 0: 

61 print(f'Read data from file {path!r}') 

62 for key in profiles.keys(): 

63 if f'profile_{key}' in avg_data.keys(): 

64 u_mean_profile[P.index(key)] = avg_data[f'profile_{key}'][P.local_slice(False)[-1]] 

65 elif comm.rank == 0: 

66 print('No mean profiles available yet. Please rerun script after completion to get correct RMS profiles') 

67 

68 # prepare progress bar 

69 indeces = range(args['restart_idx'], data.nFields) 

70 if P.comm.rank == 0: 

71 indeces = tqdm(indeces) 

72 

73 # loop through all data points and compute stuff 

74 for i in indeces: 

75 _t, u = data.readField(i) 

76 

77 # Nusselt numbers 

78 _Nu = P.compute_Nusselt_numbers(u) 

79 if any(me > 1e3 for me in _Nu.values()): 

80 continue 

81 

82 for key in Nu.keys(): 

83 Nu[key].append(_Nu[key]) 

84 

85 t.append(_t) 

86 

87 # profiles 

88 _profiles = P.get_vertical_profiles(u, list(profiles.keys())) 

89 _rms_profiles = P.get_vertical_profiles((u - u_mean_profile) ** 2, list(profiles.keys())) 

90 for key in profiles.keys(): 

91 profiles[key].append(_profiles[key]) 

92 rms_profiles[key].append(_rms_profiles[key]) 

93 

94 # spectrum 

95 k, s = P.get_frequency_spectrum(u) 

96 s_mean = zInt @ P.axes[-1].transform(s[0], axes=(0,)) 

97 spectrum.append(s_mean) 

98 spectrum_all.append(s) 

99 

100 # make a plot of the results 

101 t = xp.array(t) 

102 z = P.axes[-1].get_1dgrid() 

103 

104 converged = config.converged 

105 if config.converged == 0: 

106 print('Warning: no convergence has been set for this configuration!') 

107 if xp.max(t) < converged: 

108 converged = 0 

109 print(f'Warning: Convergence time {config.converged} has not been reached! Simulation only goes to {xp.max(t)}') 

110 

111 fig, axs = plt.subplots(1, 4, figsize=(18, 4)) 

112 for key in Nu.keys(): 

113 axs[0].plot(t, Nu[key], label=f'$Nu_{{{key}}}$') 

114 if converged > 0: 

115 axs[0].axvline(converged, color='black') 

116 axs[0].set_ylabel('$Nu$') 

117 axs[0].set_xlabel('$t$') 

118 axs[0].legend(frameon=False) 

119 

120 # compute differences in Nusselt numbers 

121 avg_Nu = {} 

122 std_Nu = {} 

123 for key in Nu.keys(): 

124 _Nu = [Nu[key][i] for i in range(len(Nu[key])) if t[i] > converged] 

125 avg_Nu[key] = xp.mean(_Nu) 

126 std_Nu[key] = xp.std(_Nu) 

127 

128 rel_error = { 

129 key: abs(avg_Nu[key] - avg_Nu['V']) / avg_Nu['V'] 

130 for key in [ 

131 't', 

132 'b', 

133 'thermal', 

134 'kinetic', 

135 ] 

136 } 

137 if comm.rank == 0: 

138 print( 

139 f'With Ra={P.Rayleigh:.0e} got Nu={avg_Nu["V"]:.2f}+-{std_Nu["V"]:.2f} with errors: Top {rel_error["t"]:.2e}, bottom: {rel_error["b"]:.2e}, thermal: {rel_error["thermal"]:.2e}, kinetic: {rel_error["kinetic"]:.2e}' 

140 ) 

141 

142 # compute average profiles 

143 avg_profiles = {} 

144 for key, values in profiles.items(): 

145 values_from_convergence = [values[i] for i in range(len(values)) if t[i] >= converged] 

146 

147 avg_profiles[key] = xp.mean(values_from_convergence, axis=0) 

148 

149 avg_rms_profiles = {} 

150 for key, values in rms_profiles.items(): 

151 values_from_convergence = [values[i] for i in range(len(values)) if t[i] >= converged] 

152 avg_rms_profiles[key] = xp.sqrt(xp.mean(values_from_convergence, axis=0)) 

153 

154 # average T 

155 avg_T = avg_profiles['T'] 

156 axs[1].axvline(0.5, color='black') 

157 axs[1].plot(avg_T, z) 

158 axs[1].set_xlabel('$T$') 

159 axs[1].set_ylabel('$z$') 

160 

161 # rms profiles 

162 avg_T = avg_rms_profiles['T'] 

163 max_idx = xp.argmax(avg_T) 

164 res_in_boundary_layer = max_idx if max_idx < len(z) / 2 else len(z) - max_idx 

165 boundary_layer = z[max_idx] if max_idx > len(z) / 2 else P.axes[-1].L - z[max_idx] 

166 if comm.rank == 0: 

167 print( 

168 f'Thermal boundary layer of thickness {boundary_layer:.2f} is resolved with {res_in_boundary_layer} points' 

169 ) 

170 axs[2].axhline(z[max_idx], color='black') 

171 axs[2].plot(avg_T, z) 

172 axs[2].scatter(avg_T, z) 

173 axs[2].set_xlabel(r'$T_\text{rms}$') 

174 axs[2].set_ylabel('$z$') 

175 

176 # spectrum 

177 _s = xp.array(spectrum) 

178 avg_spectrum = xp.mean(_s[t >= converged], axis=0) 

179 axs[3].loglog(k[avg_spectrum > 1e-15], avg_spectrum[avg_spectrum > 1e-15]) 

180 axs[3].set_xlabel('$k$') 

181 axs[3].set_ylabel(r'$\|\hat{u}_x\|$') 

182 

183 if P.comm.rank == 0: 

184 write_data = { 

185 't': t, 

186 'Nu': Nu, 

187 'avg_Nu': avg_Nu, 

188 'std_Nu': std_Nu, 

189 'z': P.axes[-1].get_1dgrid(), 

190 'k': k, 

191 'spectrum': spectrum_all, 

192 'avg_spectrum': avg_spectrum, 

193 'boundary_layer_thickness': boundary_layer, 

194 'res_in_boundary_layer': res_in_boundary_layer, 

195 } 

196 for key, values in avg_profiles.items(): 

197 write_data[f'profile_{key}'] = values 

198 for key, values in avg_rms_profiles.items(): 

199 write_data[f'rms_profile_{key}'] = values 

200 

201 with open(path, 'wb') as file: 

202 pickle.dump(write_data, file) 

203 print(f'Wrote data to file {path!r}') 

204 

205 if plot: 

206 fig.tight_layout() 

207 fig.savefig(f'{base_path}/{fname_trim}.pdf') 

208 plt.show() 

209 return path 

210 

211 

212def get_pySDC_data(res=-1, dt=-1, config_name='RBC3DG4', base_path='data/processed'): 

213 path = f'{base_path}/{config_name}-res{res}-dt{dt:.0e}.pickle' 

214 with open(path, 'rb') as file: 

215 data = pickle.load(file) 

216 

217 return data 

218 

219 

220if __name__ == '__main__': 

221 process_RBC3D_data()