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

1import os 

2import pickle 

3import numpy as np 

4from pySDC.helpers.fieldsIO import FieldsIO 

5from pySDC.projects.RayleighBenard.RBC3D_configs import get_config 

6from pySDC.implementations.problem_classes.RayleighBenard3D import RayleighBenard3D 

7from mpi4py import MPI 

8import matplotlib.pyplot as plt 

9from pySDC.helpers.plot_helper import figsize_by_journal 

10from pySDC.projects.RayleighBenard.analysis_scripts.plotting_utils import get_plotting_style, savefig 

11 

12step_sizes = { 

13 'RBC3DG4R4SDC22Ra1e5': [5e-3 * 2**i for i in range(8)], 

14 'RBC3DG4R4SDC23Ra1e5': [5e-3 * 2**i for i in range(8)], 

15 'RBC3DG4R4SDC34Ra1e5': [5e-3 * 2**i for i in range(8)], 

16 'RBC3DG4R4SDC44Ra1e5': [5e-3 * 2**i for i in range(8)], 

17 'RBC3DG4R4RKRa1e5': [5e-3 * 2**i for i in range(8)], 

18 'RBC3DG4R4EulerRa1e5': [5e-3 * 2**i for i in range(8)], 

19} 

20n_freefall_times = {} 

21 

22 

23def no_logging_hook(*args, **kwargs): 

24 return None 

25 

26 

27def get_path(args): 

28 config = get_config(args) 

29 fname = config.get_file_name() 

30 return f'{fname[:fname.index('dt')]}order.pickle' 

31 

32 

33def compute_errors(args, dts, Tend): 

34 errors = {'u': [], 'v': [], 'w': [], 'T': [], 'p': [], 'dt': []} 

35 prob = RayleighBenard3D(nx=4, ny=4, nz=4, comm=MPI.COMM_SELF) 

36 

37 dts = np.sort(dts)[::-1] 

38 ref = run(args, dts[-1], Tend) 

39 for dt in dts[:-1]: 

40 u = run(args, dt, Tend) 

41 e = u - ref 

42 for comp in ['u', 'v', 'w', 'T', 'p']: 

43 i = prob.index(comp) 

44 e_comp = np.max(np.abs(e[i])) / np.max(np.abs(ref[i])) 

45 e_comp = MPI.COMM_WORLD.allreduce(e_comp, op=MPI.MAX) 

46 errors[comp].append(float(e_comp)) 

47 errors['dt'].append(dt) 

48 

49 path = get_path(args) 

50 if MPI.COMM_WORLD.rank == 0: 

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

52 pickle.dump(errors, file) 

53 print(f'Saved errors to {path}', flush=True) 

54 

55 

56def plot_error_all_components(args): # pragma: no cover 

57 fig, ax = plt.subplots() 

58 with open(get_path(args), 'rb') as file: 

59 errors = pickle.load(file) 

60 

61 for comp in ['u', 'v', 'w', 'T', 'p']: 

62 e = np.array(errors[comp]) 

63 dt = np.array(errors['dt']) 

64 order = np.log(e[1:] / e[:-1]) / np.log(dt[1:] / dt[:-1]) 

65 ax.loglog(errors['dt'], errors[comp], label=f'{comp} order {np.mean(order):.1f}') 

66 

67 ax.loglog(errors['dt'], np.array(errors['dt']) ** 4, label='Order 4', ls='--') 

68 ax.loglog(errors['dt'], np.array(errors['dt']) ** 2, label='Order 2', ls='--') 

69 ax.legend() 

70 ax.set_xlabel(r'$\Delta t$') 

71 ax.set_ylabel(r'$e$') 

72 

73 

74def compare_order(Ra, ax=None): # pragma: no cover 

75 if ax is None: 

76 fig, ax = plt.subplots(figsize=figsize_by_journal('Nature_CS', 1, 0.6)) 

77 else: 

78 fig = None 

79 

80 if Ra == 1e5: 

81 names = ['RK', 'Euler', 'SDC23', 'SDC44'][::-1] 

82 configs = [f'RBC3DG4R4{me}Ra1e5' for me in names] 

83 paths = [f'./data/RBC3DG4R4{me}Ra1e5-res-1-order.pickle' for me in names] 

84 

85 else: 

86 raise NotImplementedError 

87 

88 for path, config in zip(paths, configs, strict=True): 

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

90 errors = pickle.load(file) 

91 

92 e = np.array(errors['T']) 

93 dt = np.array(errors['dt']) 

94 order = np.log(e[1:] / e[:-1]) / np.log(dt[1:] / dt[:-1]) 

95 print(f'{config}: order: mean={np.mean(order):.1f} / median={np.median(order):.1f}') 

96 ax.loglog(dt, e, **get_plotting_style(config)) 

97 

98 for _dt in dt: 

99 for i in [1, 2, 3, 4]: 

100 ax.text(_dt, _dt**i, i, fontweight='bold', fontsize=14, ha='center', va='center') 

101 ax.loglog(dt, dt**i, ls=':', color='black') 

102 

103 ax.legend(frameon=False) 

104 ax.set_xlabel(r'$\Delta t$') 

105 ax.set_ylabel(r'$\frac{\|T-T_\mathrm{ref}\|_\infty}{\|T_\mathrm{ref}\|_\infty}$') 

106 if fig is not None: 

107 savefig(fig, 'RBC3D_order_Ra1e5') 

108 

109 

110def run(args, dt, Tend): 

111 from pySDC.projects.RayleighBenard.run_experiment import run_experiment 

112 from pySDC.core.errors import ConvergenceError 

113 

114 args['mode'] = 'run' 

115 args['dt'] = dt 

116 

117 config = get_config(args) 

118 config.Tend = n_freefall_times.get(type(config).__name__, 3) 

119 

120 desc = config.get_description(res=args['res']) 

121 prob = desc['problem_class'](**desc['problem_params']) 

122 

123 ic_config_name = type(config).__name__ 

124 for name in ['RK', 'Euler', 'O3', 'O4', 'SDC23', 'SDC34', 'SDC44', 'SDC22']: 

125 ic_config_name = ic_config_name.replace(name, 'SDC34') 

126 

127 ic_config = get_config({**args, 'config': ic_config_name}) 

128 config.ic_config['config'] = type(ic_config) 

129 config.ic_config['res'] = ic_config.res 

130 config.ic_config['dt'] = ic_config.dt 

131 

132 config.get_LogToFile = no_logging_hook 

133 config.Tend = Tend 

134 

135 u_hat = run_experiment(args, config) 

136 u = prob.itransform(u_hat) 

137 

138 return u 

139 

140 

141if __name__ == '__main__': 

142 from pySDC.projects.RayleighBenard.run_experiment import parse_args 

143 

144 args = parse_args() 

145 

146 if args['mode'] == 'run': 

147 config = get_config(args) 

148 dts = step_sizes[type(config).__name__] 

149 compute_errors(args, dts, max(dts)) 

150 

151 if args['config'] is not None: 

152 plot_error_all_components(args) 

153 

154 compare_order(1e5) 

155 if MPI.COMM_WORLD.rank == 0: 

156 plt.show()