Coverage for pySDC / projects / StroemungsRaum / run_accuracy / run_accuracy_heat_equation_FEniCS.py: 80%

1from pathlib import Path 

2import numpy as np 

3from pySDC.projects.StroemungsRaum.run_heat_equation_FEniCS import setup, run_simulation 

4import matplotlib.pyplot as plt 

5 

6 

7def run_accuracy(c_nvars=64, num_nodes=2): 

8 """ 

9 Routine to run the accuracy test for the heat equation with FEniCS. It runs the simulation 

10 for different dt values, collects the errors, computes, and plots the observed order of accuracy. 

11 

12 Args: 

13 ----- 

14 c_nvars: int, 

15 spatial resolution (number of degrees of freedom in space) 

16 num_nodes: int, 

17 number of collocation nodes in time 

18 

19 Returns 

20 ------- 

21 results: dict 

22 Dictionary containing the errors for each dt, as well as the list of dts for easier access. 

23 order: float 

24 The observed order of accuracy, obtained by fitting a line in the log-log space to the errors and dts. 

25 """ 

26 # set parameters 

27 Tend = 1.0 

28 t0 = 0.0 

29 

30 description, controller_params = setup(t0=t0) 

31 description['problem_params']['c_nvars'] = c_nvars 

32 description['sweeper_params']['num_nodes'] = num_nodes 

33 

34 # assemble list of dt 

35 dt_list = [0.2 / 2**p for p in range(0, 4)] 

36 

37 results = {} 

38 # loop over all dt values 

39 for dt in dt_list: 

40 

41 # update dt in description for this run 

42 description['level_params']['dt'] = dt 

43 

44 # run the simulation and get the error 

45 _, _, err = run_simulation(description, controller_params, Tend) 

46 

47 # store the error in the results dictionary 

48 results[dt] = err 

49 

50 # errors in the correct order 

51 err_list = [results[dt] for dt in results.keys()] 

52 

53 # global slope (fit in log-log) 

54 order, _ = np.polyfit(np.log(dt_list), np.log(err_list), 1) 

55 

56 return results, order 

57 

58 

59def plot_accuracy(results, p=3): # pragma: no cover 

60 """ 

61 Routine to visualize the errors as well as the expected errors 

62 

63 Args: 

64 ----- 

65 results: dict 

66 The dictionary containing the errors for each dt, as well as the list of dts for easier access. 

67 p: int 

68 The expected order of accuracy 

69 """ 

70 

71 # get suffix for order (e.g. 1st, 2nd, 3rd, 4th, etc.) 

72 suffix = "th" if 10 <= p % 100 <= 13 else {1: "st", 2: "nd", 3: "rd"}.get(p % 10, "th") 

73 

74 # get the list of dt and errors from the results dictionary 

75 dt_list = sorted(results.keys()) 

76 err_list = [results[dt] for dt in dt_list] 

77 

78 # Set up plotting parameters 

79 params = { 

80 'legend.fontsize': 20, 

81 'figure.figsize': (12, 8), 

82 'axes.labelsize': 20, 

83 'axes.titlesize': 20, 

84 'xtick.labelsize': 16, 

85 'ytick.labelsize': 16, 

86 'lines.linewidth': 3, 

87 } 

88 plt.rcParams.update(params) 

89 

90 # create new figure 

91 plt.figure() 

92 # take x-axis limits from dt_list + some spacning left and right 

93 plt.xlim([min(dt_list) / 1.5, max(dt_list) * 1.5]) 

94 plt.xlabel('dt') 

95 plt.ylabel('abs. error') 

96 plt.grid() 

97 

98 # assemble optimal errors for 3rd order method and plot 

99 order_guide = [err_list[np.argmax(dt_list)] / (dt_list[np.argmax(dt_list)] / dt) ** p for dt in dt_list] 

100 

101 plt.loglog(dt_list, order_guide, color='k', ls='--', label=f"{p}{suffix} order") 

102 plt.loglog(dt_list, err_list, ls=' ', marker='o', markersize=10, label='experiment') 

103 

104 # adjust y-axis limits, add legend 

105 plt.ylim([min(order_guide) / 2, max(order_guide) * 2]) 

106 plt.legend(loc=2, ncol=1, numpoints=1) 

107 

108 plt.grid(True, which="minor", ls="--", color='0.8') 

109 plt.grid(True, which="major", ls="-", color='0.001') 

110 

111 # get the data directory 

112 import os 

113 

114 path = f"{os.path.dirname(__file__)}/../data/heat_equation/" 

115 

116 # if it does not exist, create the 'data' directory at the specified path, including any necessary parent directories 

117 Path(path).mkdir(parents=True, exist_ok=True) 

118 fname = path + f'heat_equation_{p}{suffix}_order_time_FEniCS.png' 

119 plt.savefig(fname, bbox_inches='tight') 

120 

121 return None 

122 

123 

124def main(): 

125 

126 # parameters for 3rd order accuracy test 

127 c_nvars = 64 # spatial resolution 

128 num_nodes = 2 # number of collocation nodes in time 

129 p = 2 * num_nodes - 1 # expected order of accuracy 

130 

131 # run the simulation and get the results 

132 results, _ = run_accuracy(c_nvars, num_nodes) 

133 plot_accuracy(results, p) 

134 

135 

136if __name__ == "__main__": 

137 main()