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

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 

31 # prepare paths 

32 os.makedirs(base_path, exist_ok=True) 

33 fname = config.get_file_name() 

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

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

36 

37 # open simulation data 

38 data = FieldsIO.fromFile(fname) 

39 

40 # prepare arrays to store data in 

41 Nu = { 

42 'V': [], 

43 'b': [], 

44 't': [], 

45 'thermal': [], 

46 'kinetic': [], 

47 } 

48 t = [] 

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

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

51 spectrum = [] 

52 spectrum_all = [] 

53 

54 # try to load time averaged values 

55 u_mean_profile = P.u_exact() 

56 if os.path.isfile(path): 

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

58 avg_data = pickle.load(file) 

59 if comm.rank == 0: 

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

61 for key in profiles.keys(): 

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

63 if P.useGPU: 

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

65 else: 

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

67 elif comm.rank == 0: 

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

69 

70 # prepare progress bar 

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

72 if P.comm.rank == 0: 

73 indeces = tqdm(indeces) 

74 

75 # loop through all data points and compute stuff 

76 for i in indeces: 

77 _t, u = data.readField(i) 

78 

79 if P.useGPU: 

80 u = P.xp.array(u) 

81 

82 # Nusselt numbers 

83 _Nu = P.compute_Nusselt_numbers(u) 

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

85 continue 

86 

87 for key in Nu.keys(): 

88 if P.useGPU: 

89 Nu[key].append(_Nu[key].get()) 

90 else: 

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

92 

93 t.append(_t) 

94 

95 # profiles 

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

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

98 for key in profiles.keys(): 

99 if P.useGPU: 

100 profiles[key].append(_profiles[key].get()) 

101 rms_profiles[key].append(_rms_profiles[key].get()) 

102 else: 

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

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

105 

106 # spectrum 

107 k, s = P.get_frequency_spectrum(u) 

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

109 if P.useGPU: 

110 spectrum.append(s_mean.get()) 

111 spectrum_all.append(s.get()) 

112 else: 

113 spectrum.append(s_mean) 

114 spectrum_all.append(s) 

115 

116 # make a plot of the results 

117 t = np.array(t) 

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

119 if P.useGPU: 

120 z = z.get() 

121 k = k.get() 

122 

123 converged = config.converged 

124 if config.converged == 0: 

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

126 if np.max(t) < converged: 

127 converged = 0 

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

129 

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

131 for key in Nu.keys(): 

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

133 if converged > 0: 

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

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

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

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

138 

139 # compute differences in Nusselt numbers 

140 avg_Nu = {} 

141 std_Nu = {} 

142 for key in Nu.keys(): 

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

144 avg_Nu[key] = np.mean(_Nu) 

145 std_Nu[key] = np.std(_Nu) 

146 

147 rel_error = { 

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

149 for key in [ 

150 't', 

151 'b', 

152 'thermal', 

153 'kinetic', 

154 ] 

155 } 

156 if comm.rank == 0: 

157 print( 

158 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}' 

159 ) 

160 

161 # compute average profiles 

162 avg_profiles = {} 

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

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

165 

166 avg_profiles[key] = np.mean(values_from_convergence, axis=0) 

167 

168 avg_rms_profiles = {} 

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

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

171 avg_rms_profiles[key] = np.sqrt(np.mean(values_from_convergence, axis=0)) 

172 

173 # average T 

174 avg_T = avg_profiles['T'] 

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

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

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

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

179 

180 # rms profiles 

181 avg_T = avg_rms_profiles['T'] 

182 max_idx = np.argmax(avg_T) 

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

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

185 if comm.rank == 0: 

186 print( 

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

188 ) 

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

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

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

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

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

194 

195 # spectrum 

196 _s = np.array(spectrum) 

197 avg_spectrum = np.mean(_s[t >= converged], axis=0) 

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

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

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

201 

202 if P.comm.rank == 0: 

203 write_data = { 

204 't': t, 

205 'Nu': Nu, 

206 'avg_Nu': avg_Nu, 

207 'std_Nu': std_Nu, 

208 'z': z, 

209 'k': k, 

210 'spectrum': spectrum_all, 

211 'avg_spectrum': avg_spectrum, 

212 'boundary_layer_thickness': int(boundary_layer), 

213 'res_in_boundary_layer': int(res_in_boundary_layer), 

214 } 

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

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

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

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

219 

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

221 pickle.dump(write_data, file) 

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

223 

224 if plot: 

225 fig.tight_layout() 

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

227 plt.show() 

228 return path 

229 

230 

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

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

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

234 data = pickle.load(file) 

235 

236 return data 

237 

238 

239if __name__ == '__main__': 

240 process_RBC3D_data()