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 

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

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

60uRefs = {nSteps: solutionExact(tEnd, nSteps, "PROTHERO-ROBINSON", epsilon=epsilon) for nSteps in nStepsList} 

61 

62i = 0 

63for qDeltaList, nSweeps in config: 

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

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

66 i += 1 

67 

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

69 try: 

70 params = getParamsRK(qDelta) 

71 except KeyError: 

72 params = getParamsSDC( 

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

74 ) 

75 

76 errors = [] 

77 costs = [] 

78 

79 for nSteps in nStepsList: 

80 uRef = uRefs[nSteps] 

81 

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

83 

84 err = getError(uSDC, uRef) 

85 errors.append(err) 

86 

87 cost = getCost(counters) 

88 if parallel: 

89 cost /= nNodes * parEfficiency 

90 costs.append(cost) 

91 

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

93 

94 plt.figure(figNameConv) 

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

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

97 

98 plt.figure(figNameCost) 

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

100 

101 for figName in [figNameConv, figNameCost]: 

102 plt.figure(figName) 

103 plt.gca().set( 

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

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

106 ) 

107 plt.legend() 

108 plt.grid(True) 

109 plt.tight_layout() 

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

111 

112data.writeToJSON()