Coverage for pySDC/projects/parallelSDC_reloaded/jacobiElliptic_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 the JacobianElliptic problem 

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 

19# Problem parameters 

20tEnd = 10 

21pName = "JACELL" 

22 

23 

24def getError(uNum, uRef): 

25 if uNum is None: # pragma: no cover 

26 return np.inf 

27 return np.linalg.norm(uRef[-1] - uNum[-1], 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 'ESDIRK53', 

44 'ESDIRK43', 

45 'PIC', 

46 # 'IE', 'LU', 'IEpar', 'PIC', 

47 'MIN-SR-NS', 

48 'MIN-SR-S', 

49 'MIN-SR-FLEX', 

50 # "MIN3", 

51] 

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

53# nSweepList = [1, 2, 3, 4] 

54 

55# qDeltaList = ['RK4', 'ESDIRK43', 'MIN-SR-S'] 

56nSweepList = [4] 

57 

58 

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

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

61 

62dtVals = tEnd / nStepsList 

63 

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

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

66uRefs = {nSteps: solutionExact(tEnd, nSteps, pName) for nSteps in nStepsList} 

67 

68i = 0 

69for qDelta in qDeltaList: 

70 for nSweeps in nSweepList: 

71 sym = symList[i] 

72 i += 1 

73 

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

75 try: 

76 params = getParamsRK(qDelta) 

77 name = name[:-3] 

78 if nSweeps != nSweepList[0]: # pragma: no cover 

79 continue 

80 except KeyError: 

81 params = getParamsSDC( 

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

83 ) 

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

85 

86 errors = [] 

87 costs = [] 

88 

89 for nSteps in nStepsList: 

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

91 

92 uRef = uRefs[nSteps] 

93 

94 uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, pName) 

95 

96 err = getError(uSDC, uRef) 

97 errors.append(err) 

98 

99 cost = getCost(counters) 

100 if parallel: 

101 cost /= nNodes * parEfficiency 

102 costs.append(cost) 

103 

104 # error VS dt 

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

106 # error VS cost 

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

108 

109for i in range(2): 

110 axs[i].set( 

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

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

113 # ylim=(1e-9, 1e0), 

114 ) 

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

116 axs[i].grid() 

117 

118fig.set_size_inches(12, 5) 

119fig.tight_layout()