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

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 = 1 

22pName = "KAPS" 

23pParams = { 

24 "epsilon": 1e-6, 

25} 

26 

27 

28def getError(uNum, uRef): 

29 if uNum is None: 

30 return np.inf 

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

32 

33 

34def getCost(counters): 

35 nNewton, nRHS, tComp = counters 

36 return nNewton + nRHS 

37 

38 

39# Base variable parameters 

40nNodes = 4 

41quadType = 'RADAU-RIGHT' 

42nodeType = 'LEGENDRE' 

43parEfficiency = 0.8 # 1/nNodes 

44 

45qDeltaList = [ 

46 'RK4', 

47 'ESDIRK53', 

48 'DIRK43', 

49 # 'IE', 'LU', 'IEpar', 'PIC', 

50 'MIN-SR-NS', 

51 'MIN-SR-S', 

52 'MIN-SR-FLEX', 

53 # "MIN3", 

54] 

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

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

57 

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

59nSweepList = [4] 

60 

61 

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

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

64 

65dtVals = tEnd / nStepsList 

66 

67i = 0 

68for qDelta in qDeltaList: 

69 for nSweeps in nSweepList: 

70 sym = symList[i] 

71 i += 1 

72 

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

74 try: 

75 params = getParamsRK(qDelta) 

76 name = name[:-3] 

77 except KeyError: 

78 params = getParamsSDC( 

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

80 ) 

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

82 

83 errors = [] 

84 costs = [] 

85 

86 for nSteps in nStepsList: 

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

88 

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

90 

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

92 

93 err = getError(uSDC, uRef) 

94 errors.append(err) 

95 

96 cost = getCost(counters) 

97 if parallel: 

98 cost /= nNodes * parEfficiency 

99 costs.append(cost) 

100 

101 # error VS dt 

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

103 # error VS cost 

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

105 

106for i in range(2): 

107 axs[i].set( 

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

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

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

111 ) 

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

113 axs[i].grid() 

114 

115fig.set_size_inches(12, 5) 

116fig.tight_layout()