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

13 

14import numpy as np 

15import matplotlib.pyplot as plt 

16 

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

18 

19# Problem parameters 

20tEnd = 50 

21pName = "ALLEN-CAHN" 

22periodic = False 

23pParams = { 

24 "periodic": periodic, 

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

26 "epsilon": 0.04, 

27} 

28 

29 

30def getError(uNum, uRef): 

31 if uNum is None: 

32 return np.inf 

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

34 

35 

36def getCost(counters): 

37 nNewton, nRHS, tComp = counters 

38 return 2 * nNewton + nRHS 

39 

40 

41# Base variable parameters 

42nNodes = 4 

43quadType = 'RADAU-RIGHT' 

44nodeType = 'LEGENDRE' 

45parEfficiency = 1 / nNodes 

46 

47qDeltaList = [ 

48 'RK4', 

49 'ESDIRK53', 

50 'VDHS', 

51 'MIN', 

52 # 'IE', 'LU', 'IEpar', 'PIC', 

53 'MIN-SR-NS', 

54 'MIN-SR-S', 

55 'MIN-SR-FLEX', 

56 "PIC", 

57 # "MIN3", 

58] 

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

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

61 

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

63nSweepList = [4] 

64 

65 

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

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

68 

69dtVals = tEnd / nStepsList 

70 

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

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

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

74 

75i = 0 

76for qDelta in qDeltaList: 

77 for nSweeps in nSweepList: 

78 sym = symList[i] 

79 i += 1 

80 

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

82 try: 

83 params = getParamsRK(qDelta) 

84 name = name[:-3] 

85 except KeyError: 

86 params = getParamsSDC( 

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

88 ) 

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

90 

91 errors = [] 

92 costs = [] 

93 

94 for nSteps in nStepsList: 

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

96 

97 uRef = uRefs[nSteps] 

98 

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

100 

101 err = getError(uSDC, uRef) 

102 errors.append(err) 

103 

104 cost = getCost(counters) 

105 if parallel: 

106 cost /= nNodes * parEfficiency 

107 costs.append(cost) 

108 

109 # error VS dt 

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

111 # error VS cost 

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

113 

114for i in range(2): 

115 axs[i].set( 

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

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

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

119 ) 

120 axs[i].legend() 

121 axs[i].grid() 

122 

123fig.set_size_inches(12, 5) 

124fig.tight_layout()