Coverage for pySDC/projects/parallelSDC_reloaded/kaps_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 = 1 

21pName = "KAPS" 

22pParams = { 

23 "epsilon": 1e-6, 

24} 

25 

26 

27def getError(uNum, uRef): 

28 if uNum is None: 

29 return np.inf 

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

31 

32 

33def getCost(counters): 

34 nNewton, nRHS, tComp = counters 

35 return nNewton + nRHS 

36 

37 

38# Base variable parameters 

39nNodes = 4 

40quadType = 'RADAU-RIGHT' 

41nodeType = 'LEGENDRE' 

42parEfficiency = 0.8 # 1/nNodes 

43 

44qDeltaList = [ 

45 'RK4', 

46 'ESDIRK53', 

47 'DIRK43', 

48 # 'IE', 'LU', 'IEpar', 'PIC', 

49 'MIN-SR-NS', 

50 'MIN-SR-S', 

51 'MIN-SR-FLEX', 

52 # "MIN3", 

53] 

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

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

56 

57# qDeltaList = ['MIN-SR-S'] 

58nSweepList = [4] 

59 

60 

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

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

63 

64dtVals = tEnd / nStepsList 

65 

66i = 0 

67for qDelta in qDeltaList: 

68 for nSweeps in nSweepList: 

69 sym = symList[i] 

70 i += 1 

71 

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

73 try: 

74 params = getParamsRK(qDelta) 

75 name = name[:-3] 

76 except KeyError: 

77 params = getParamsSDC( 

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

79 ) 

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

81 

82 errors = [] 

83 costs = [] 

84 

85 for nSteps in nStepsList: 

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

87 

88 uRef = solutionExact(tEnd, nSteps, pName, **pParams) 

89 

90 uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, pName, **pParams) 

91 

92 err = getError(uSDC, uRef) 

93 errors.append(err) 

94 

95 cost = getCost(counters) 

96 if parallel: 

97 cost /= nNodes * parEfficiency 

98 costs.append(cost) 

99 

100 # error VS dt 

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

102 # error VS cost 

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

104 

105for i in range(2): 

106 axs[i].set( 

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

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

109 ylim=(1e-17, 1e0), 

110 ) 

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

112 axs[i].grid() 

113 

114fig.set_size_inches(12, 5) 

115fig.tight_layout()