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

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 

27 

28epsilon = 1e-3 

29 

30# ----------------------------------------------------------------------------- 

31# %% Convergence and error VS cost plots 

32# ----------------------------------------------------------------------------- 

33tEnd = 2 * np.pi 

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

35dtVals = tEnd / nStepsList 

36 

37 

38def getError(uNum, uRef): 

39 if uNum is None: 

40 return np.inf 

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

42 

43 

44def getCost(counters): 

45 nNewton, nRHS, tComp = counters 

46 return nNewton + nRHS 

47 

48 

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

50 

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

52config = [ 

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

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

55] 

56 

57 

58i = 0 

59for qDeltaList, nSweeps in config: 

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

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

62 i += 1 

63 

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

65 try: 

66 params = getParamsRK(qDelta) 

67 except KeyError: 

68 params = getParamsSDC( 

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

70 ) 

71 

72 errors = [] 

73 costs = [] 

74 

75 for nSteps in nStepsList: 

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

77 

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

79 

80 err = getError(uSDC, uRef) 

81 errors.append(err) 

82 

83 cost = getCost(counters) 

84 if parallel: 

85 cost /= nNodes * parEfficiency 

86 costs.append(cost) 

87 

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

89 

90 plt.figure(figNameConv) 

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

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

93 

94 plt.figure(figNameCost) 

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

96 

97 for figName in [figNameConv, figNameCost]: 

98 plt.figure(figName) 

99 plt.gca().set( 

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

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

102 ) 

103 plt.legend() 

104 plt.grid(True) 

105 plt.tight_layout() 

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

107 

108data.writeToJSON()