Coverage for pySDC/projects/compression/order.py: 77%

1import numpy as np 

2import matplotlib.pyplot as plt 

3 

4from pySDC.projects.Resilience.advection import run_advection 

5 

6from pySDC.helpers.stats_helper import get_sorted 

7from pySDC.helpers.plot_helper import figsize_by_journal 

8from pySDC.implementations.hooks.log_errors import LogGlobalErrorPostRun 

9 

10from pySDC.projects.compression.compression_convergence_controller import Compression 

11 

12MACHINEPRECISION = ( 

13 1e-8 # generous tolerance below which we ascribe errors to floating point rounding errors rather than compression 

14) 

15LOGGER_LEVEL = 30 

16 

17 

18def single_run(problem, description=None, thresh=1e-10, Tend=2e-1, useMPI=False, num_procs=1): 

19 description = {} if description is None else description 

20 

21 compressor_args = {} 

22 compressor_args['compressor_config'] = {'pressio:abs': thresh} 

23 

24 if thresh > 0: 

25 description['convergence_controllers'] = {Compression: {'compressor_args': compressor_args}} 

26 

27 controller_params = { 

28 'mssdc_jac': False, 

29 'logger_level': LOGGER_LEVEL, 

30 } 

31 

32 stats, _, _ = problem( 

33 custom_description=description, 

34 hook_class=LogGlobalErrorPostRun, 

35 Tend=Tend, 

36 use_MPI=useMPI, 

37 num_procs=num_procs, 

38 custom_controller_params=controller_params, 

39 ) 

40 

41 if useMPI: 

42 from mpi4py import MPI 

43 

44 comm = MPI.COMM_WORLD 

45 else: 

46 comm = None 

47 e = min([me[1] for me in get_sorted(stats, type='e_global_post_run', comm=comm)]) 

48 return e 

49 

50 

51def multiple_runs(problem, values, expected_order, mode='dt', thresh=1e-10, useMPI=False, num_procs=1, **kwargs): 

52 errors = np.zeros_like(values) 

53 

54 description = { 

55 'level_params': {}, 

56 'problam_params': {}, 

57 'step_params': {}, 

58 } 

59 if mode == 'dt': 

60 description['step_params'] = {'maxiter': expected_order} 

61 elif mode == 'nvars': 

62 description['problem_params'] = {'order': expected_order} 

63 

64 for i in range(len(values)): 

65 if mode == 'dt': 

66 description['level_params']['dt'] = values[i] 

67 Tend = values[i] * (5 if num_procs == 1 else 2 * num_procs) 

68 elif mode == 'nvars': 

69 description['problem_params']['nvars'] = values[i] 

70 Tend = 2e-1 

71 

72 errors[i] = single_run(problem, description, thresh=thresh, Tend=Tend, useMPI=useMPI, num_procs=num_procs) 

73 return values, errors 

74 

75 

76def get_order(values, errors, thresh=1e-16, expected_order=None): 

77 values = np.array(values) 

78 idx = np.argsort(values) 

79 local_orders = np.log(errors[idx][1:] / errors[idx][:-1]) / np.log(values[idx][1:] / values[idx][:-1]) 

80 order = np.mean(local_orders[errors[idx][1:] > max([thresh, MACHINEPRECISION])]) 

81 if expected_order is not None: 

82 assert np.isclose(order, expected_order, atol=0.5), f"Expected order {expected_order}, but got {order:.2f}!" 

83 return order 

84 

85 

86def plot_order(values, errors, ax, thresh=1e-16, color='black', expected_order=None, **kwargs): 

87 values = np.array(values) 

88 order = get_order(values, errors, thresh=thresh, expected_order=expected_order) 

89 ax.scatter(values, errors, color=color, **kwargs) 

90 ax.loglog(values, errors[0] * (values / values[0]) ** order, color=color, label=f'p={order:.2f}', **kwargs) 

91 

92 

93def plot_order_in_time(ax, thresh, useMPI=False, num_procs=1): 

94 problem = run_advection 

95 

96 base_configs_dt = { 

97 'values': np.array([2.0 ** (-i) for i in [2, 3, 4, 5, 6, 7, 8, 9]]), 

98 'mode': 'dt', 

99 'ax': ax, 

100 'thresh': thresh, 

101 } 

102 

103 configs_dt = {} 

104 configs_dt[2] = {**base_configs_dt, 'color': 'black'} 

105 configs_dt[3] = {**base_configs_dt, 'color': 'magenta'} 

106 configs_dt[4] = {**base_configs_dt, 'color': 'teal'} 

107 configs_dt[5] = {**base_configs_dt, 'color': 'orange'} 

108 # configs_dt[6] = {**base_configs_dt, 'color': 'blue'} 

109 

110 for key in configs_dt.keys(): 

111 values, errors = multiple_runs( 

112 problem, expected_order=key, useMPI=useMPI, **configs_dt[key], num_procs=num_procs 

113 ) 

114 plot_order( 

115 values, 

116 errors, 

117 ax=configs_dt[key]['ax'], 

118 thresh=configs_dt[key]['thresh'] * 1e2, 

119 color=configs_dt[key]['color'], 

120 expected_order=key + 1, 

121 ) 

122 base_configs_dt['ax'].set_xlabel(r'$\Delta t$') 

123 base_configs_dt['ax'].set_ylabel('local error') 

124 base_configs_dt['ax'].axhline( 

125 base_configs_dt['thresh'], color='grey', ls='--', label=rf'$\|\delta\|={ {thresh:.0e}} $' 

126 ) 

127 base_configs_dt['ax'].legend(frameon=False) 

128 

129 

130def order_in_time_different_error_bounds(): 

131 fig, axs = plt.subplots( 

132 2, 2, figsize=figsize_by_journal('Springer_Numerical_Algorithms', 1.0, 1.0), sharex=True, sharey=True 

133 ) 

134 threshs = [1e-6, 1e-8, 1e-10, 1e-12] 

135 

136 for i in range(len(threshs)): 

137 ax = axs.flatten()[i] 

138 plot_order_in_time(ax, threshs[i]) 

139 if i != 2: 

140 ax.set_ylabel('') 

141 ax.set_xlabel('') 

142 fig.suptitle('Order in time for advection problem') 

143 fig.tight_layout() 

144 fig.savefig('compression-order-time.pdf') 

145 

146 

147if __name__ == '__main__': 

148 order_in_time_different_error_bounds() 

149 

150 # base_configs_nvars = { 

151 # 'values': [128, 256, 512, 1024], 

152 # # 'values': np.array([2**(i) for i in [4, 5, 6, 7, 8, 9]]), 

153 # 'mode': 'nvars', 

154 # } 

155 

156 # configs_nvars = {} 

157 # configs_nvars[2] = {**base_configs_nvars, 'color': 'black'} 

158 # configs_nvars[4] = {**base_configs_nvars, 'color': 'magenta'} 

159 

160 # for key in configs_nvars.keys(): 

161 # values, errors = multiple_runs(problem, expected_order=key, **configs_nvars[key]) 

162 # plot_order(values, errors, axs[1], color=configs_nvars[key]['color']) 

163 

164 plt.show()