Coverage for pySDC/projects/parallelSDC_reloaded/allenCahn_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 Allen-Cahn problem : 

7 

8- error VS time-step 

9- error VS computation cost 

10 

11Note : implementation in progress ... 

12""" 

13import numpy as np 

14import matplotlib.pyplot as plt 

15 

16from pySDC.projects.parallelSDC_reloaded.utils import getParamsSDC, getParamsRK, solutionSDC, solutionExact 

17 

18# Problem parameters 

19tEnd = 50 

20pName = "ALLEN-CAHN" 

21periodic = False 

22pParams = { 

23 "periodic": periodic, 

24 "nvars": 2**11 - (not periodic), 

25 "epsilon": 0.04, 

26} 

27 

28 

29def getError(uNum, uRef): 

30 if uNum is None: 

31 return np.inf 

32 return np.linalg.norm(uRef[-1, :] - uNum[-1, :], ord=2) 

33 

34 

35def getCost(counters): 

36 nNewton, nRHS, tComp = counters 

37 return 2 * nNewton + nRHS 

38 

39 

40# Base variable parameters 

41nNodes = 4 

42quadType = 'RADAU-RIGHT' 

43nodeType = 'LEGENDRE' 

44parEfficiency = 1 / nNodes 

45 

46qDeltaList = [ 

47 'RK4', 

48 'ESDIRK53', 

49 'VDHS', 

50 'MIN', 

51 # 'IE', 'LU', 'IEpar', 'PIC', 

52 'MIN-SR-NS', 

53 'MIN-SR-S', 

54 'MIN-SR-FLEX', 

55 "PIC", 

56 # "MIN3", 

57] 

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

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

60 

61qDeltaList = ['ESDIRK43', 'MIN-SR-FLEX'] 

62nSweepList = [4] 

63 

64 

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

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

67 

68dtVals = tEnd / nStepsList 

69 

70i = 0 

71for qDelta in qDeltaList: 

72 for nSweeps in nSweepList: 

73 sym = symList[i] 

74 i += 1 

75 

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

77 try: 

78 params = getParamsRK(qDelta) 

79 name = name[:-3] 

80 except KeyError: 

81 params = getParamsSDC( 

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

83 ) 

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

85 

86 errors = [] 

87 costs = [] 

88 

89 for nSteps in nStepsList: 

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

91 

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

93 

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

95 

96 err = getError(uSDC, uRef) 

97 errors.append(err) 

98 

99 cost = getCost(counters) 

100 if parallel: 

101 cost /= nNodes * parEfficiency 

102 costs.append(cost) 

103 

104 # error VS dt 

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

106 # error VS cost 

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

108 

109for i in range(2): 

110 axs[i].set( 

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

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

113 ylim=(1e-5, 1e1), 

114 ) 

115 axs[i].legend() 

116 axs[i].grid() 

117 

118fig.set_size_inches(12, 5) 

119fig.tight_layout()