Coverage for pySDC/projects/parallelSDC_reloaded/protheroRobinson_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 ProtheroRobinson 

7(linear and non-linear) problem : 

8 

9- error VS time-step 

10- error VS computation cost 

11 

12Note : implementation in progress ... 

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 = 2 * np.pi 

21nonLinear = False 

22epsilon = 1e-3 

23collUpdate = False 

24initSweep = "copy" 

25 

26pName = "PROTHERO-ROBINSON" + (nonLinear) * "-NL" 

27 

28 

29def getError(uNum, uRef): 

30 if uNum is None: 

31 return np.inf 

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

33 

34 

35def getCost(counters): 

36 nNewton, nRHS, tComp = counters 

37 return nNewton + nRHS 

38 

39 

40# Base variable parameters 

41nNodes = 4 

42quadType = 'RADAU-RIGHT' 

43nodeType = 'LEGENDRE' 

44parEfficiency = 0.8 

45 

46qDeltaList = [ 

47 'MIN-SR-NS', 

48 'MIN-SR-S', 

49 'MIN-SR-FLEX', 

50 "ESDIRK43", 

51 'LU', 

52 'VDHS', 

53 # 'IE', 'LU', 'IEpar', 'PIC', 

54 # "MIN3", 

55] 

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

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

58 

59# qDeltaList = ['ESDIRK43', 'ESDIRK53', 'VDHS'] 

60nSweepList = [6] 

61 

62 

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

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

65 

66dtVals = tEnd / 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.split('(')[0] 

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

79 continue 

80 

81 except KeyError: 

82 params = getParamsSDC( 

83 quadType=quadType, 

84 numNodes=nNodes, 

85 nodeType=nodeType, 

86 qDeltaI=qDelta, 

87 nSweeps=nSweeps, 

88 collUpdate=collUpdate, 

89 initType=initSweep, 

90 ) 

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

92 

93 errors = [] 

94 costs = [] 

95 

96 for nSteps in nStepsList: 

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

98 

99 uRef = solutionExact(tEnd, nSteps, pName, epsilon=epsilon) 

100 

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

102 

103 err = getError(uSDC, uRef) 

104 errors.append(err) 

105 

106 cost = getCost(counters) 

107 if parallel: 

108 cost /= nNodes * parEfficiency 

109 costs.append(cost) 

110 

111 # error VS dt 

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

113 # error VS cost 

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

115 

116for i in range(2): 

117 axs[i].set( 

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

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

120 ylim=(1e-12, 1e3), 

121 ) 

122 axs[i].legend() 

123 axs[i].grid() 

124 

125fig.set_size_inches(12, 5) 

126fig.tight_layout()