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): # pragma: no cover 

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

76 if Ra == 1e5: 

77 names = ['RK', 'Euler', 'SDC22', 'SDC23', 'SDC34', 'SDC44'][::-1] 

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

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

80 

81 else: 

82 raise NotImplementedError 

83 

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

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

86 errors = pickle.load(file) 

87 

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

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

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

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

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

93 

94 for _dt in dt: 

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

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

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

98 

99 ax.legend(frameon=False) 

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

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

102 savefig(fig, 'RBC3D_order_Ra1e5') 

103 

104 

105def run(args, dt, Tend): 

106 from pySDC.projects.RayleighBenard.run_experiment import run_experiment 

107 from pySDC.core.errors import ConvergenceError 

108 

109 args['mode'] = 'run' 

110 args['dt'] = dt 

111 

112 config = get_config(args) 

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

114 

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

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

117 

118 ic_config_name = type(config).__name__ 

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

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

121 

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

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

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

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

126 

127 config.get_LogToFile = no_logging_hook 

128 config.Tend = Tend 

129 

130 u_hat = run_experiment(args, config) 

131 u = prob.itransform(u_hat) 

132 

133 return u 

134 

135 

136if __name__ == '__main__': 

137 from pySDC.projects.RayleighBenard.run_experiment import parse_args 

138 

139 args = parse_args() 

140 

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

142 config = get_config(args) 

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

144 compute_errors(args, dts, max(dts)) 

145 

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

147 plot_error_all_components(args) 

148 

149 compare_order(1e5) 

150 if MPI.COMM_WORLD.rank == 0: 

151 plt.show()