## 54 statements

, created at 2024-09-20 17:10 +0000

1#!/usr/bin/env python3

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

3"""

4Created on Tue Dec 5 11:02:39 2023

6Script to investigate diagonal SDC on the ProtheroRobinson

7(linear and non-linear) problem :

9- error VS time-step

10- error VS computation cost

12Note : implementation in progress ...

13"""

14import numpy as np

15import matplotlib.pyplot as plt

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

19# Problem parameters

20tEnd = 300

21pName = "CHEMREC"

24def getError(uNum, uRef):

25 if uNum is None: # pragma: no cover

26 return np.inf

27 return max(np.linalg.norm(uRef[:, 0] - uNum[:, 0], np.inf), np.linalg.norm(uRef[:, 1] - uNum[:, 1], np.inf))

30def getCost(counters):

31 nNewton, nRHS, tComp = counters

32 return nNewton + nRHS

35# Base variable parameters

36nNodes = 4

38nodeType = 'LEGENDRE'

39parEfficiency = 1 / nNodes

41qDeltaList = [

42 'RK4',

43 'ESDIRK53',

44 'ESDIRK43',

45 # 'IE', 'LU', 'IEpar', 'PIC',

46 'MIN-SR-NS',

47 'MIN-SR-S',

48 'MIN-SR-FLEX',

49 # "MIN3",

50]

51nStepsList = np.array([2, 5, 10, 20])

52nSweepList = [1, 2, 3, 4, 5, 6]

54qDeltaList = ['ESDIRK43', 'MIN-SR-S', 'MIN-SR-FLEX']

55nSweepList = [4]

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

59fig, axs = plt.subplots(1, 2)

61dtVals = tEnd / nStepsList

63i = 0

64for qDelta in qDeltaList:

65 for nSweeps in nSweepList:

66 sym = symList[i]

67 i += 1

69 name = f"{qDelta}({nSweeps})"

70 try:

71 params = getParamsRK(qDelta)

72 name = name[:-3]

73 except KeyError:

74 params = getParamsSDC(

76 )

77 print(f'computing for {name} ...')

79 errors = []

80 costs = []

82 for nSteps in nStepsList:

83 print(f' -- nSteps={nSteps} ...')

85 uRef = solutionExact(tEnd, nSteps, pName)

87 uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, pName)

89 err = getError(uSDC, uRef)

90 errors.append(err)

92 cost = getCost(counters)

93 if parallel:

94 cost /= nNodes * parEfficiency

95 costs.append(cost)

97 # error VS dt

98 axs[0].loglog(dtVals, errors, sym + '-', label=name)

99 # error VS cost

100 axs[1].loglog(costs, errors, sym + '-', label=name)

102for i in range(2):

103 axs[i].set(

104 xlabel=r"$\Delta{t}$" if i == 0 else "cost",

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

106 ylim=(1e-9, 1e0),

107 )

108 axs[i].legend(loc="lower right" if i == 0 else "lower left")

109 axs[i].grid()

111fig.set_size_inches(12, 5)

112fig.tight_layout()