Coverage for pySDC / projects / parallelSDC_reloaded / chemicalReaction_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 = 300 

22pName = "CHEMREC" 

23 

24 

25def getError(uNum, uRef): 

26 if uNum is None: # pragma: no cover 

27 return np.inf 

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

29 

30 

31def getCost(counters): 

32 nNewton, nRHS, tComp = counters 

33 return nNewton + nRHS 

34 

35 

36# Base variable parameters 

37nNodes = 4 

38quadType = 'RADAU-RIGHT' 

39nodeType = 'LEGENDRE' 

40parEfficiency = 1 / nNodes 

41 

42qDeltaList = [ 

43 'RK4', 

44 'ESDIRK53', 

45 'ESDIRK43', 

46 # 'IE', 'LU', 'IEpar', 'PIC', 

47 'MIN-SR-NS', 

48 'MIN-SR-S', 

49 'MIN-SR-FLEX', 

50 # "MIN3", 

51] 

52nStepsList = np.array([2, 5, 10, 20]) 

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

54 

55qDeltaList = ['ESDIRK43', 'MIN-SR-S', 'MIN-SR-FLEX'] 

56nSweepList = [4] 

57 

58 

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

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

61 

62dtVals = tEnd / nStepsList 

63 

64i = 0 

65for qDelta in qDeltaList: 

66 for nSweeps in nSweepList: 

67 sym = symList[i] 

68 i += 1 

69 

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

71 try: 

72 params = getParamsRK(qDelta) 

73 name = name[:-3] 

74 except KeyError: 

75 params = getParamsSDC( 

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

77 ) 

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

79 

80 errors = [] 

81 costs = [] 

82 

83 for nSteps in nStepsList: 

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

85 

86 uRef = solutionExact(tEnd, nSteps, pName) 

87 

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

89 

90 err = getError(uSDC, uRef) 

91 errors.append(err) 

92 

93 cost = getCost(counters) 

94 if parallel: 

95 cost /= nNodes * parEfficiency 

96 costs.append(cost) 

97 

98 # error VS dt 

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

100 # error VS cost 

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

102 

103for i in range(2): 

104 axs[i].set( 

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

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

107 ylim=(1e-9, 1e0), 

108 ) 

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

110 axs[i].grid() 

111 

112fig.set_size_inches(12, 5) 

113fig.tight_layout()