## 52 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 JacobianElliptic problem

8- error VS time-step

9- error VS computation cost

11Note : implementation in progress ...

12"""

13import numpy as np

14import matplotlib.pyplot as plt

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

18# Problem parameters

19tEnd = 10

20pName = "JACELL"

23def getError(uNum, uRef):

24 if uNum is None: # pragma: no cover

25 return np.inf

26 return np.linalg.norm(uRef[-1] - uNum[-1], np.inf)

29def getCost(counters):

30 nNewton, nRHS, tComp = counters

31 return nNewton + nRHS

34# Base variable parameters

35nNodes = 4

37nodeType = 'LEGENDRE'

38parEfficiency = 1 / nNodes

40qDeltaList = [

41 'RK4',

42 'ESDIRK53',

43 'ESDIRK43',

44 'PIC',

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

46 'MIN-SR-NS',

47 'MIN-SR-S',

48 'MIN-SR-FLEX',

49 # "MIN3",

50]

51nStepsList = np.array([10, 20, 50, 100, 200])

52# nSweepList = [1, 2, 3, 4]

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

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 if nSweeps != nSweepList[0]: # pragma: no cover

74 continue

75 except KeyError:

76 params = getParamsSDC(

78 )

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

81 errors = []

82 costs = []

84 for nSteps in nStepsList:

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

87 uRef = solutionExact(tEnd, nSteps, pName)

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

91 err = getError(uSDC, uRef)

92 errors.append(err)

94 cost = getCost(counters)

95 if parallel:

96 cost /= nNodes * parEfficiency

97 costs.append(cost)

99 # error VS dt

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

101 # error VS cost

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

104for i in range(2):

105 axs[i].set(

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

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

108 # ylim=(1e-9, 1e0),

109 )

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

111 axs[i].grid()

113fig.set_size_inches(12, 5)

114fig.tight_layout()