56 statements

, created at 2024-09-20 16:55 +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 = 1

21pName = "KAPS"

22pParams = {

23 "epsilon": 1e-6,

24}

27def getError(uNum, uRef):

28 if uNum is None:

29 return np.inf

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

33def getCost(counters):

34 nNewton, nRHS, tComp = counters

35 return nNewton + nRHS

38# Base variable parameters

39nNodes = 4

41nodeType = 'LEGENDRE'

42parEfficiency = 0.8 # 1/nNodes

44qDeltaList = [

45 'RK4',

46 'ESDIRK53',

47 'DIRK43',

48 # 'IE', 'LU', 'IEpar', 'PIC',

49 'MIN-SR-NS',

50 'MIN-SR-S',

51 'MIN-SR-FLEX',

52 # "MIN3",

53]

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

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

57# qDeltaList = ['MIN-SR-S']

58nSweepList = [4]

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

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

64dtVals = tEnd / nStepsList

66i = 0

67for qDelta in qDeltaList:

68 for nSweeps in nSweepList:

69 sym = symList[i]

70 i += 1

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

73 try:

74 params = getParamsRK(qDelta)

75 name = name[:-3]

76 except KeyError:

77 params = getParamsSDC(

79 )

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

82 errors = []

83 costs = []

85 for nSteps in nStepsList:

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

88 uRef = solutionExact(tEnd, nSteps, pName, **pParams)

90 uSDC, counters, parallel = solutionSDC(tEnd, nSteps, params, pName, **pParams)

92 err = getError(uSDC, uRef)

93 errors.append(err)

95 cost = getCost(counters)

96 if parallel:

97 cost /= nNodes * parEfficiency

98 costs.append(cost)

100 # error VS dt

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

102 # error VS cost

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

105for i in range(2):

106 axs[i].set(

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

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

109 ylim=(1e-17, 1e0),

110 )

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

112 axs[i].grid()

114fig.set_size_inches(12, 5)

115fig.tight_layout()