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 

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

83# instead of once per (qDelta, nSteps). 

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

85 

86i = 0 

87for qDeltaList in config: 

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

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

90 i += 1 

91 

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

93 try: 

94 params = getParamsRK(qDelta) 

95 except KeyError: 

96 params = getParamsSDC( 

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

98 ) 

99 

100 errors = [] 

101 costs = [] 

102 

103 for nSteps in nStepsList: 

104 uRef = uRefs[nSteps] 

105 

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

107 

108 err = getError(uSDC, uRef) 

109 errors.append(err) 

110 

111 cost = getCost(counters) 

112 if parallel: 

113 cost /= nNodes * parEfficiency 

114 costs.append(cost) 

115 

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

117 

118 plt.figure(figNameConv) 

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

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

121 

122 plt.figure(figNameCost) 

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

124 

125 for figName in [figNameConv, figNameCost]: 

126 plt.figure(figName) 

127 plt.gca().set( 

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

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

130 ) 

131 plt.legend() 

132 plt.grid(True) 

133 plt.tight_layout() 

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

135 

136data.writeToJSON()