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 

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

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

69uRefs = {nSteps: solutionExact(tEnd, nSteps, pName, **pParams) for nSteps in nStepsList} 

70 

71i = 0 

72for qDelta in qDeltaList: 

73 for nSweeps in nSweepList: 

74 sym = symList[i] 

75 i += 1 

76 

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

78 try: 

79 params = getParamsRK(qDelta) 

80 name = name[:-3] 

81 except KeyError: 

82 params = getParamsSDC( 

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

84 ) 

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

86 

87 errors = [] 

88 costs = [] 

89 

90 for nSteps in nStepsList: 

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

92 

93 uRef = uRefs[nSteps] 

94 

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

96 

97 err = getError(uSDC, uRef) 

98 errors.append(err) 

99 

100 cost = getCost(counters) 

101 if parallel: 

102 cost /= nNodes * parEfficiency 

103 costs.append(cost) 

104 

105 # error VS dt 

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

107 # error VS cost 

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

109 

110for i in range(2): 

111 axs[i].set( 

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

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

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

115 ) 

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

117 axs[i].grid() 

118 

119fig.set_size_inches(12, 5) 

120fig.tight_layout()