## 58 statements

, created at 2024-09-09 14:59 +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 Allen-Cahn 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 = 50

20pName = "ALLEN-CAHN"

21periodic = False

22pParams = {

23 "periodic": periodic,

24 "nvars": 2**11 - (not periodic),

25 "epsilon": 0.04,

26}

29def getError(uNum, uRef):

30 if uNum is None:

31 return np.inf

32 return np.linalg.norm(uRef[-1, :] - uNum[-1, :], ord=2)

35def getCost(counters):

36 nNewton, nRHS, tComp = counters

37 return 2 * nNewton + nRHS

40# Base variable parameters

41nNodes = 4

43nodeType = 'LEGENDRE'

44parEfficiency = 1 / nNodes

46qDeltaList = [

47 'RK4',

48 'ESDIRK53',

49 'VDHS',

50 'MIN',

51 # 'IE', 'LU', 'IEpar', 'PIC',

52 'MIN-SR-NS',

53 'MIN-SR-S',

54 'MIN-SR-FLEX',

55 "PIC",

56 # "MIN3",

57]

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

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

61qDeltaList = ['ESDIRK43', 'MIN-SR-FLEX']

62nSweepList = [4]

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

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

68dtVals = tEnd / nStepsList

70i = 0

71for qDelta in qDeltaList:

72 for nSweeps in nSweepList:

73 sym = symList[i]

74 i += 1

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

77 try:

78 params = getParamsRK(qDelta)

79 name = name[:-3]

80 except KeyError:

81 params = getParamsSDC(

83 )

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

86 errors = []

87 costs = []

89 for nSteps in nStepsList:

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

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

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

96 err = getError(uSDC, uRef)

97 errors.append(err)

99 cost = getCost(counters)

100 if parallel:

101 cost /= nNodes * parEfficiency

102 costs.append(cost)

104 # error VS dt

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

106 # error VS cost

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

109for i in range(2):

110 axs[i].set(

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

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

113 ylim=(1e-5, 1e1),

114 )

115 axs[i].legend()

116 axs[i].grid()

118fig.set_size_inches(12, 5)

119fig.tight_layout()