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 

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

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

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

67 

68i = 0 

69for qDelta in qDeltaList: 

70 for nSweeps in nSweepList: 

71 sym = symList[i] 

72 i += 1 

73 

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

75 try: 

76 params = getParamsRK(qDelta) 

77 name = name[:-3] 

78 except KeyError: 

79 params = getParamsSDC( 

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

81 ) 

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

83 

84 errors = [] 

85 costs = [] 

86 

87 for nSteps in nStepsList: 

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

89 

90 uRef = uRefs[nSteps] 

91 

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

93 

94 err = getError(uSDC, uRef) 

95 errors.append(err) 

96 

97 cost = getCost(counters) 

98 if parallel: 

99 cost /= nNodes * parEfficiency 

100 costs.append(cost) 

101 

102 # error VS dt 

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

104 # error VS cost 

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

106 

107for i in range(2): 

108 axs[i].set( 

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

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

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

112 ) 

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

114 axs[i].grid() 

115 

116fig.set_size_inches(12, 5) 

117fig.tight_layout()