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

1#!/usr/bin/env python3 

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

3""" 

4Created on Thu Jan 11 10:21:47 2024 

5 

6Figures with experiment on the Prothero-Robinson 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 

26 

27epsilon = 1e-3 

28 

29# ----------------------------------------------------------------------------- 

30# %% Convergence and error VS cost plots 

31# ----------------------------------------------------------------------------- 

32tEnd = 2 * np.pi 

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

34dtVals = tEnd / nStepsList 

35 

36 

37def getError(uNum, uRef): 

38 if uNum is None: 

39 return np.inf 

40 return np.linalg.norm(np.linalg.norm(uRef - uNum, np.inf, axis=-1), np.inf) 

41 

42 

43def getCost(counters): 

44 nNewton, nRHS, tComp = counters 

45 return nNewton + nRHS 

46 

47 

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

49 

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

51config = [ 

52 [(*minPrec, "VDHS", "ESDIRK43", "LU"), 4], 

53 [(*minPrec, "VDHS", "ESDIRK43", "LU"), 6], 

54] 

55 

56 

57i = 0 

58for qDeltaList, nSweeps in config: 

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

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

61 i += 1 

62 

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

64 try: 

65 params = getParamsRK(qDelta) 

66 except KeyError: 

67 params = getParamsSDC( 

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

69 ) 

70 

71 errors = [] 

72 costs = [] 

73 

74 for nSteps in nStepsList: 

75 uRef = solutionExact(tEnd, nSteps, "PROTHERO-ROBINSON", epsilon=epsilon) 

76 

77 uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, "PROTHERO-ROBINSON", epsilon=epsilon) 

78 

79 err = getError(uSDC, uRef) 

80 errors.append(err) 

81 

82 cost = getCost(counters) 

83 if parallel: 

84 cost /= nNodes * parEfficiency 

85 costs.append(cost) 

86 

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

88 

89 plt.figure(figNameConv) 

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

91 data.storeAndCheck(f"{figNameConv}_{qDelta}", errors) 

92 

93 plt.figure(figNameCost) 

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

95 

96 for figName in [figNameConv, figNameCost]: 

97 plt.figure(figName) 

98 plt.gca().set( 

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

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

101 ) 

102 plt.legend() 

103 plt.grid(True) 

104 plt.tight_layout() 

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

106 

107data.writeToJSON()