Coverage for pySDC/projects/parallelSDC_reloaded/lorenz_accuracy.py: 100%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Created on Tue Dec 5 11:02:39 2023 

5 

6Script to investigate diagonal SDC on Lorenz system 

7 

8- error VS time-step 

9- error VS computation cost 

10 

11Note : implementation in progress ... 

12""" 

13 

14import numpy as np 

15import matplotlib.pyplot as plt 

16 

17from pySDC.projects.parallelSDC_reloaded.utils import getParamsSDC, getParamsRK, solutionSDC, solutionExact 

18 

19tEnd = 1.24 

20 

21 

22def getError(uNum, uRef): 

23 if uNum is None: # pragma: no cover 

24 return np.inf 

25 return np.linalg.norm(np.linalg.norm(uRef - uNum, np.inf, axis=-1), np.inf) 

26 

27 

28def getCost(counters): 

29 nNewton, nRHS, tComp = counters 

30 return nNewton + nRHS 

31 

32 

33# Base variable parameters 

34nNodes = 4 

35quadType = 'RADAU-RIGHT' 

36nodeType = 'LEGENDRE' 

37parEfficiency = 0.8 # 1/nNodes 

38 

39qDeltaList = [ 

40 'RK4', 

41 'ESDIRK53', 

42 'VDHS', 

43 'MIN', 

44 # 'IE', 'LU', 'IEpar', 'PIC', 

45 'MIN-SR-NS', 

46 'MIN-SR-S', 

47 'MIN-SR-FLEX', 

48 "PIC", 

49 # "MIN3", 

50] 

51nStepsList = np.array([2, 5, 10, 20, 50, 100, 200]) 

52nSweepList = [1, 2, 3, 4] 

53 

54 

55symList = ['o', '^', 's', '>', '*', '<', 'p', '>'] * 10 

56 

57qDeltaList = ['MIN-SR-S', 'RK4'] 

58# nSweepList = [4] 

59 

60fig, axs = plt.subplots(1, 2) 

61 

62 

63dtVals = tEnd / nStepsList 

64 

65# The reference solution depends only on nSteps, so compute it once per nSteps 

66# instead of once per (qDelta, nSweeps, nSteps). 

67uRefs = {nSteps: solutionExact(tEnd, nSteps, "LORENZ", u0=(5, -5, 20)) for nSteps in nStepsList} 

68 

69i = 0 

70for qDelta in qDeltaList: 

71 for nSweeps in nSweepList: 

72 sym = symList[i] 

73 i += 1 

74 

75 name = f"{qDelta}({nSweeps})" 

76 try: 

77 params = getParamsRK(qDelta) 

78 name = name[:-3] 

79 if nSweeps != nSweepList[0]: 

80 continue 

81 except KeyError: 

82 params = getParamsSDC( 

83 quadType=quadType, numNodes=nNodes, nodeType=nodeType, qDeltaI=qDelta, nSweeps=nSweeps 

84 ) 

85 print(f'computing for {name} ...') 

86 

87 errors = [] 

88 costs = [] 

89 

90 for nSteps in nStepsList: 

91 print(f' -- nSteps={nSteps} ...') 

92 

93 uRef = uRefs[nSteps] 

94 

95 uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, "LORENZ", u0=(5, -5, 20)) 

96 

97 err = getError(uSDC, uRef) 

98 errors.append(err) 

99 

100 cost = getCost(counters) 

101 if parallel: 

102 cost /= nNodes * parEfficiency 

103 costs.append(cost) 

104 

105 # error VS dt 

106 axs[0].loglog(dtVals, errors, sym + '-', label=name) 

107 # error VS cost 

108 axs[1].loglog(costs, errors, sym + '-', label=name) 

109 

110x = dtVals[4:] 

111for k in [1, 2, 3, 4, 5]: 

112 axs[0].loglog(x, 1e4 * x**k, "--", color="gray", linewidth=0.8) 

113 

114for i in range(2): 

115 axs[i].set( 

116 xlabel=r"$\Delta{t}$" if i == 0 else "cost", 

117 ylabel=r"$L_\infty$ error", 

118 ylim=(8.530627786509715e-12, 372.2781393394293), 

119 ) 

120 axs[i].legend(loc="lower right" if i == 0 else "lower left") 

121 axs[i].grid() 

122 

123fig.set_size_inches(12, 5) 

124fig.tight_layout()