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""" 

14 

15import numpy as np 

16import matplotlib.pyplot as plt 

17 

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

19 

20# Problem parameters 

21tEnd = 2 * np.pi 

22nonLinear = False 

23epsilon = 1e-3 

24collUpdate = False 

25initSweep = "copy" 

26 

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

28 

29 

30def getError(uNum, uRef): 

31 if uNum is None: 

32 return np.inf 

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

34 

35 

36def getCost(counters): 

37 nNewton, nRHS, tComp = counters 

38 return nNewton + nRHS 

39 

40 

41# Base variable parameters 

42nNodes = 4 

43quadType = 'RADAU-RIGHT' 

44nodeType = 'LEGENDRE' 

45parEfficiency = 0.8 

46 

47qDeltaList = [ 

48 'MIN-SR-NS', 

49 'MIN-SR-S', 

50 'MIN-SR-FLEX', 

51 "ESDIRK43", 

52 'LU', 

53 'VDHS', 

54 # 'IE', 'LU', 'IEpar', 'PIC', 

55 # "MIN3", 

56] 

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

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

59 

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

61nSweepList = [6] 

62 

63 

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

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

66 

67dtVals = tEnd / nStepsList 

68 

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

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

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

72 

73i = 0 

74for qDelta in qDeltaList: 

75 for nSweeps in nSweepList: 

76 sym = symList[i] 

77 i += 1 

78 

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

80 try: 

81 params = getParamsRK(qDelta) 

82 name = name.split('(')[0] 

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

84 continue 

85 

86 except KeyError: 

87 params = getParamsSDC( 

88 quadType=quadType, 

89 numNodes=nNodes, 

90 nodeType=nodeType, 

91 qDeltaI=qDelta, 

92 nSweeps=nSweeps, 

93 collUpdate=collUpdate, 

94 initType=initSweep, 

95 ) 

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

97 

98 errors = [] 

99 costs = [] 

100 

101 for nSteps in nStepsList: 

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

103 

104 uRef = uRefs[nSteps] 

105 

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

107 

108 err = getError(uSDC, uRef) 

109 errors.append(err) 

110 

111 cost = getCost(counters) 

112 if parallel: 

113 cost /= nNodes * parEfficiency 

114 costs.append(cost) 

115 

116 # error VS dt 

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

118 # error VS cost 

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

120 

121for i in range(2): 

122 axs[i].set( 

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

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

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

126 ) 

127 axs[i].legend() 

128 axs[i].grid() 

129 

130fig.set_size_inches(12, 5) 

131fig.tight_layout()