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

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 if config.converged == 0: 

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

106 

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

108 for key in Nu.keys(): 

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

110 if config.converged > 0: 

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

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

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

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

115 

116 # compute differences in Nusselt numbers 

117 avg_Nu = {} 

118 std_Nu = {} 

119 for key in Nu.keys(): 

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

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

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

123 

124 rel_error = { 

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

126 for key in [ 

127 't', 

128 'b', 

129 'thermal', 

130 'kinetic', 

131 ] 

132 } 

133 if comm.rank == 0: 

134 print( 

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

136 ) 

137 

138 # compute average profiles 

139 avg_profiles = {} 

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

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

142 

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

144 

145 avg_rms_profiles = {} 

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

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

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

149 

150 # average T 

151 avg_T = avg_profiles['T'] 

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

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

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

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

156 

157 # rms profiles 

158 avg_T = avg_rms_profiles['T'] 

159 max_idx = xp.argmax(avg_T) 

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

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

162 if comm.rank == 0: 

163 print( 

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

165 ) 

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

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

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

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

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

171 

172 # spectrum 

173 _s = xp.array(spectrum) 

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

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

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

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

178 

179 if P.comm.rank == 0: 

180 write_data = { 

181 't': t, 

182 'Nu': Nu, 

183 'avg_Nu': avg_Nu, 

184 'std_Nu': std_Nu, 

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

186 'k': k, 

187 'spectrum': spectrum_all, 

188 'avg_spectrum': avg_spectrum, 

189 'boundary_layer_thickness': boundary_layer, 

190 'res_in_boundary_layer': res_in_boundary_layer, 

191 } 

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

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

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

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

196 

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

198 pickle.dump(write_data, file) 

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

200 

201 if plot: 

202 fig.tight_layout() 

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

204 plt.show() 

205 return path 

206 

207 

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

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

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

211 data = pickle.load(file) 

212 

213 return data 

214 

215 

216if __name__ == '__main__': 

217 process_RBC3D_data()