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

8import os 

9import numpy as np 

10 

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

12from pySDC.helpers.testing import DataChecker 

13 

14data = DataChecker(__file__) 

15 

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

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

18 

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

20 

21# SDC parameters 

22nNodes = 4 

23quadType = 'RADAU-RIGHT' 

24nodeType = 'LEGENDRE' 

25parEfficiency = 0.8 # 1/nNodes 

26nSweeps = 4 

27 

28# Problem parameters 

29pName = "ALLEN-CAHN" 

30tEnd = 50 

31pParams = { 

32 "periodic": False, 

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

34 "epsilon": 0.04, 

35} 

36 

37# ----------------------------------------------------------------------------- 

38# Trajectories (reference solution) 

39# ----------------------------------------------------------------------------- 

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

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

42 

43figName = f"{SCRIPT}_solution" 

44plt.figure(figName) 

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

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

47 

48plt.legend() 

49plt.xlabel("$x$") 

50plt.ylabel("Solution") 

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

52plt.tight_layout() 

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

54 

55# ----------------------------------------------------------------------------- 

56# %% Convergence and error VS cost plots 

57# ----------------------------------------------------------------------------- 

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

59dtVals = tEnd / nStepsList 

60 

61 

62def getError(uNum, uRef): 

63 if uNum is None: 

64 return np.inf 

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

66 

67 

68def getCost(counters): 

69 nNewton, nRHS, tComp = counters 

70 return 2 * nNewton + nRHS 

71 

72 

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

74 

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

76config = [ 

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

78] 

79 

80 

81i = 0 

82for qDeltaList in config: 

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

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

85 i += 1 

86 

87 for qDelta, sym in zip(qDeltaList, symList): 

88 try: 

89 params = getParamsRK(qDelta) 

90 except KeyError: 

91 params = getParamsSDC( 

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

93 ) 

94 

95 errors = [] 

96 costs = [] 

97 

98 for nSteps in nStepsList: 

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

100 

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

102 

103 err = getError(uSDC, uRef) 

104 errors.append(err) 

105 

106 cost = getCost(counters) 

107 if parallel: 

108 cost /= nNodes * parEfficiency 

109 costs.append(cost) 

110 

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

112 

113 plt.figure(figNameConv) 

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

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

116 

117 plt.figure(figNameCost) 

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

119 

120 for figName in [figNameConv, figNameCost]: 

121 plt.figure(figName) 

122 plt.gca().set( 

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

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

125 ) 

126 plt.legend() 

127 plt.grid(True) 

128 plt.tight_layout() 

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

130 

131data.writeToJSON()