Coverage for pySDC / projects / parallelSDC_reloaded / scripts / fig03_lorenz.py: 100%

1#!/usr/bin/env python3 

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

3""" 

4Created on Wed Jan 10 15:34:24 2024 

5 

6Figures with experiment on the Lorenz problem 

7""" 

8 

9import os 

10import numpy as np 

11import scipy as sp 

12 

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

14from pySDC.helpers.testing import DataChecker 

15 

16data = DataChecker(__file__) 

17 

18PATH = '/' + os.path.join(*__file__.split('/')[:-1]) 

19SCRIPT = __file__.split('/')[-1].split('.')[0] 

20 

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

22 

23# SDC parameters 

24nNodes = 4 

25quadType = 'RADAU-RIGHT' 

26nodeType = 'LEGENDRE' 

27parEfficiency = 0.8 # 1/nNodes 

28 

29# ----------------------------------------------------------------------------- 

30# Trajectories (reference solution) 

31# ----------------------------------------------------------------------------- 

32tEnd = 2 

33nSteps = tEnd * 50 

34tVals = np.linspace(0, tEnd, nSteps + 1) 

35nPeriods = 2 

36 

37print(f"Computing exact solution up to t={tEnd} ...") 

38uExact = solutionExact(tEnd, nSteps, "LORENZ", u0=(5, -5, 20)) 

39 

40z = uExact[:, -1] 

41idx = sp.signal.find_peaks(z)[0][nPeriods - 1] 

42print(f'tEnd for {nPeriods} periods : {tVals[idx]}') 

43 

44figName = f"{SCRIPT}_traj" 

45plt.figure(figName) 

46me = 0.1 

47plt.plot(tVals, uExact[:, 0], 's-', label="$x(t)$", markevery=me) 

48plt.plot(tVals, uExact[:, 1], 'o-', label="$y(t)$", markevery=me) 

49plt.plot(tVals, uExact[:, 2], '^-', label="$z(t)$", markevery=me) 

50plt.vlines(tVals[idx], ymin=-20, ymax=40, linestyles="--", linewidth=1) 

51plt.legend(loc="upper right") 

52plt.xlabel("$t$") 

53plt.ylabel("Trajectory") 

54plt.gcf().set_size_inches(12, 3) 

55plt.tight_layout() 

56plt.savefig(f'{PATH}/{figName}.pdf') 

57 

58# ----------------------------------------------------------------------------- 

59# %% Convergence plots 

60# ----------------------------------------------------------------------------- 

61tEnd = 1.24 

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

63dtVals = tEnd / nStepsList 

64 

65 

66def getError(uNum, uRef): 

67 if uNum is None: # pragma: no cover 

68 return np.inf 

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

70 

71 

72config = ["PIC", "MIN-SR-NS"] 

73for qDelta, sym in zip(config, symList, strict=False): 

74 figName = f"{SCRIPT}_conv_{qDelta}" 

75 plt.figure(figName) 

76 

77 for nSweeps in [1, 2, 3, 4, 5]: 

78 params = getParamsSDC(quadType=quadType, numNodes=nNodes, nodeType=nodeType, qDeltaI=qDelta, nSweeps=nSweeps) 

79 

80 errors = [] 

81 

82 for nSteps in nStepsList: 

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

84 

85 uRef = solutionExact(tEnd, nSteps, "LORENZ", u0=(5, -5, 20)) 

86 

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

88 

89 err = getError(uSDC, uRef) 

90 errors.append(err) 

91 

92 # error VS dt 

93 label = f"$K={nSweeps}$" 

94 plt.loglog(dtVals, errors, sym + '-', label=f"$K={nSweeps}$") 

95 data.storeAndCheck(f"{figName}_{label}", errors[1:]) 

96 

97 x = dtVals[4:] 

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

99 plt.loglog(x, 1e4 * x**k, "--", color="gray", linewidth=0.8) 

100 

101 plt.gca().set( 

102 xlabel=r"$\Delta{t}$", 

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

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

105 ) 

106 plt.legend(loc="lower right") 

107 plt.grid() 

108 plt.tight_layout() 

109 plt.savefig(f"{PATH}/{figName}.pdf") 

110 

111 

112# ----------------------------------------------------------------------------- 

113# %% Error VS cost plots 

114# ----------------------------------------------------------------------------- 

115def getCost(counters): 

116 nNewton, nRHS, tComp = counters 

117 return nNewton + nRHS 

118 

119 

120minPrec = ["MIN-SR-NS", "MIN-SR-S", "MIN-SR-FLEX"] 

121 

122symList = ['^', '>', '<', 'o', 's', '*'] 

123config = [ 

124 [(*minPrec, "LU", "EE", "PIC"), 4], 

125 [(*minPrec, "VDHS", "RK4", "ESDIRK43"), 4], 

126 [(*minPrec, "PIC", "RK4", "ESDIRK43"), 5], 

127] 

128 

129 

130i = 0 

131for qDeltaList, nSweeps in config: 

132 figName = f"{SCRIPT}_cost_{i}" 

133 i += 1 

134 plt.figure(figName) 

135 

136 for qDelta, sym in zip(qDeltaList, symList, strict=False): 

137 try: 

138 params = getParamsRK(qDelta) 

139 except KeyError: 

140 params = getParamsSDC( 

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

142 ) 

143 

144 errors = [] 

145 costs = [] 

146 

147 for nSteps in nStepsList: 

148 uRef = solutionExact(tEnd, nSteps, "LORENZ", u0=(5, -5, 20)) 

149 

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

151 

152 err = getError(uSDC, uRef) 

153 errors.append(err) 

154 

155 cost = getCost(counters) 

156 if parallel: 

157 assert qDelta != "EE", "wait, whaaat ??" 

158 cost /= nNodes * parEfficiency 

159 costs.append(cost) 

160 

161 # error VS cost 

162 ls = '-' if qDelta.startswith("MIN-SR-") else "--" 

163 plt.loglog(costs, errors, sym + ls, label=qDelta) 

164 data.storeAndCheck(f"{figName}_{qDelta}", errors[2:], rtol=1e-2) 

165 

166 plt.gca().set( 

167 xlabel="Cost", 

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

169 ylim=(1e-10, 400), 

170 xlim=(30, 20000), 

171 ) 

172 plt.legend(loc="lower left") 

173 plt.grid() 

174 plt.tight_layout() 

175 plt.savefig(f"{PATH}/{figName}.pdf") 

176 

177data.writeToJSON()