## 58 statements

, created at 2024-09-19 09:13 +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, using the autonomous formulation.

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 = 2 * np.pi

21nonLinear = True

22epsilon = 0.001

23collUpdate = False

26pName = "PROTHERO-ROBINSON-A" + (nonLinear) * "-NL"

29def getError(uNum, uRef):

30 if uNum is None:

31 return np.inf

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

35def getCost(counters):

36 nNewton, nRHS, tComp = counters

37 return nNewton + nRHS

40# Base variable parameters

41nNodes = 4

43nodeType = 'LEGENDRE'

44parEfficiency = 0.8

46qDeltaList = [

47 'MIN-SR-NS',

48 'MIN-SR-S',

49 'MIN-SR-FLEX',

50 "ESDIRK43",

51 'LU',

52 'VDHS',

53 # 'IE', 'LU', 'IEpar', 'PIC',

54 # "MIN3",

55]

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

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

59# qDeltaList = ['ESDIRK43', 'ESDIRK53', 'VDHS']

60nSweepList = [4]

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

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

66dtVals = tEnd / nStepsList

68i = 0

69for qDelta in qDeltaList:

70 for nSweeps in nSweepList:

71 sym = symList[i]

72 i += 1

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

75 try:

76 params = getParamsRK(qDelta)

77 name = name.split('(')[0]

78 if nSweeps != nSweepList[0]: # pragma: no cover

79 continue

81 except KeyError:

82 params = getParamsSDC(

84 numNodes=nNodes,

85 nodeType=nodeType,

86 qDeltaI=qDelta,

87 nSweeps=nSweeps,

88 collUpdate=collUpdate,

89 initType=initSweep,

90 )

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

93 errors = []

94 costs = []

96 for nSteps in nStepsList:

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

99 uRef = solutionExact(tEnd, nSteps, pName, epsilon=epsilon)

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

103 err = getError(uSDC, uRef)

104 errors.append(err)

106 cost = getCost(counters)

107 if parallel:

108 cost /= nNodes * parEfficiency

109 costs.append(cost)

111 # error VS dt

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

113 # error VS cost

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

116for i in range(2):

117 axs[i].set(

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

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

120 ylim=(1e-12, 1e3),

121 )

122 axs[i].legend()

123 axs[i].grid()

125fig.set_size_inches(12, 5)

126fig.tight_layout()