Coverage for pySDC / projects / parallelSDC_reloaded / scripts / fig05_allenCahn.py: 100%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Created on Thu Jan 11 11:14:01 2024 

5 

6Figures with experiments on the Allen-Cahn problem 

7""" 

8 

9import os 

10import numpy as np 

11 

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

13from pySDC.helpers.testing import DataChecker 

14 

15data = DataChecker(__file__) 

16 

17PATH = '/' + os.path.join(*__file__.split('/')[:-1]) 

18SCRIPT = __file__.split('/')[-1].split('.')[0] 

19 

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

21 

22# SDC parameters 

23nNodes = 4 

24quadType = 'RADAU-RIGHT' 

25nodeType = 'LEGENDRE' 

26parEfficiency = 0.8 # 1/nNodes 

27nSweeps = 4 

28 

29# Problem parameters 

30pName = "ALLEN-CAHN" 

31tEnd = 50 

32pParams = { 

33 "periodic": False, 

34 "nvars": 2**11 - 1, 

35 "epsilon": 0.04, 

36} 

37 

38# ----------------------------------------------------------------------------- 

39# Trajectories (reference solution) 

40# ----------------------------------------------------------------------------- 

41uExact = solutionExact(tEnd, 1, pName, **pParams) 

42x = np.linspace(-0.5, 0.5, 2**11 + 1)[1:-1] 

43 

44figName = f"{SCRIPT}_solution" 

45plt.figure(figName) 

46plt.plot(x, uExact[0, :], '-', label="$u(0)$") 

47plt.plot(x, uExact[-1, :], '--', label="$u(T)$") 

48 

49plt.legend() 

50plt.xlabel("$x$") 

51plt.ylabel("Solution") 

52plt.gcf().set_size_inches(12, 3) 

53plt.tight_layout() 

54plt.savefig(f"{PATH}/{figName}.pdf") 

55 

56# ----------------------------------------------------------------------------- 

57# %% Convergence and error VS cost plots 

58# ----------------------------------------------------------------------------- 

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

60dtVals = tEnd / nStepsList 

61 

62 

63def getError(uNum, uRef): 

64 if uNum is None: 

65 return np.inf 

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

67 

68 

69def getCost(counters): 

70 nNewton, nRHS, tComp = counters 

71 return 2 * nNewton + nRHS 

72 

73 

74minPrec = ["MIN-SR-NS", "MIN-SR-S", "MIN-SR-FLEX"] 

75 

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

77config = [ 

78 (*minPrec, "VDHS", "ESDIRK43", "LU"), 

79] 

80 

81 

82i = 0 

83for qDeltaList in config: 

84 figNameConv = f"{SCRIPT}_conv_{i}" 

85 figNameCost = f"{SCRIPT}_cost_{i}" 

86 i += 1 

87 

88 for qDelta, sym in zip(qDeltaList, symList, strict=False): 

89 try: 

90 params = getParamsRK(qDelta) 

91 except KeyError: 

92 params = getParamsSDC( 

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

94 ) 

95 

96 errors = [] 

97 costs = [] 

98 

99 for nSteps in nStepsList: 

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

101 

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

103 

104 err = getError(uSDC, uRef) 

105 errors.append(err) 

106 

107 cost = getCost(counters) 

108 if parallel: 

109 cost /= nNodes * parEfficiency 

110 costs.append(cost) 

111 

112 ls = '-' if qDelta.startswith("MIN-SR-") else "--" 

113 

114 plt.figure(figNameConv) 

115 plt.loglog(dtVals, errors, sym + ls, label=qDelta) 

116 data.storeAndCheck(f"{figNameConv}_{qDelta}", errors, atol=1e-4, rtol=1e-4) 

117 

118 plt.figure(figNameCost) 

119 plt.loglog(costs, errors, sym + ls, label=qDelta) 

120 

121 for figName in [figNameConv, figNameCost]: 

122 plt.figure(figName) 

123 plt.gca().set( 

124 xlabel="Cost" if "cost" in figName else r"$\Delta {t}$", 

125 ylabel=r"$L_2$ error at $T$", 

126 ) 

127 plt.legend() 

128 plt.grid(True) 

129 plt.tight_layout() 

130 plt.savefig(f"{PATH}/{figName}.pdf") 

131 

132data.writeToJSON()