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

1#!/usr/bin/env python3 

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

3""" 

4Created on Sun Nov 12 22:14:03 2023 

5 

6Script to investigate diagonal SDC on Van der Pol with different mu parameters, 

7in particular with graphs such as : 

8 

9- error VS time-step 

10- error VS computation cost 

11 

12Note : implementation in progress ... 

13""" 

14 

15import numpy as np 

16import matplotlib.pyplot as plt 

17 

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

19 

20muVals = [0.1, 2, 10] 

21tEndVals = [6.3, 7.6, 18.9] # tEnd = 1 period for each mu 

22 

23 

24def getError(uNum, uRef): 

25 if uNum is None: 

26 return np.inf 

27 return np.linalg.norm(uRef[:, 0] - uNum[:, 0], np.inf) 

28 

29 

30def getCost(counters): 

31 nNewton, nRHS, tComp = counters 

32 return nNewton + nRHS 

33 

34 

35# Base variable parameters 

36nNodes = 4 

37quadType = 'RADAU-RIGHT' 

38nodeType = 'LEGENDRE' 

39parEfficiency = 1 / nNodes 

40 

41qDeltaList = [ 

42 'RK4', 

43 'ESDIRK43', 

44 'LU', 

45 # 'IE', 'LU', 'IEpar', 'PIC', 

46 'MIN-SR-NS', 

47 'MIN-SR-S', 

48 'MIN-SR-FLEX', 

49] 

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

51nSweepList = [1, 2, 3, 4, 5, 6] 

52 

53 

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

55 

56# qDeltaList = ['LU'] 

57nSweepList = [4] 

58 

59fig, axs = plt.subplots(2, len(muVals)) 

60 

61for j, (mu, tEnd) in enumerate(zip(muVals, tEndVals, strict=True)): 

62 print("-" * 80) 

63 print(f"mu={mu}") 

64 print("-" * 80) 

65 

66 dtVals = tEnd / nStepsList 

67 

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

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

70 uRefs = {nSteps: solutionExact(tEnd, nSteps, "VANDERPOL", mu=mu) for nSteps in nStepsList} 

71 

72 i = 0 

73 for qDelta in qDeltaList: 

74 for nSweeps in nSweepList: 

75 sym = symList[i] 

76 i += 1 

77 

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

79 try: 

80 params = getParamsRK(qDelta) 

81 name = name[:-3] 

82 except KeyError: 

83 params = getParamsSDC( 

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

85 ) 

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

87 

88 errors = [] 

89 costs = [] 

90 

91 for nSteps in nStepsList: 

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

93 

94 uRef = uRefs[nSteps] 

95 

96 uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, "VANDERPOL", mu=mu) 

97 

98 err = getError(uSDC, uRef) 

99 errors.append(err) 

100 

101 cost = getCost(counters) 

102 if parallel: 

103 cost /= nNodes * parEfficiency 

104 costs.append(cost) 

105 

106 # error VS dt 

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

108 # error VS cost 

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

110 

111 for i in range(2): 

112 if i == 0: 

113 axs[i, j].set_title(f"mu={mu}") 

114 axs[i, j].set( 

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

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

117 ylim=(1e-11, 10), 

118 ) 

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

120 axs[i, j].grid() 

121 

122fig.set_size_inches(18.2, 10.4) 

123fig.tight_layout()