75 statements

, created at 2024-09-19 09:13 +0000

1#!/usr/bin/env python3

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

3"""

4Created on Thu Jan 11 11:14:01 2024

6Figures with experiments on the Allen-Cahn problem

7"""

8import os

9import numpy as np

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

12from pySDC.helpers.testing import DataChecker

14data = DataChecker(__file__)

16PATH = '/' + os.path.join(*__file__.split('/')[:-1])

17SCRIPT = __file__.split('/')[-1].split('.')[0]

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

21# SDC parameters

22nNodes = 4

24nodeType = 'LEGENDRE'

25parEfficiency = 0.8 # 1/nNodes

26nSweeps = 4

28# Problem parameters

29pName = "ALLEN-CAHN"

30tEnd = 50

31pParams = {

32 "periodic": False,

33 "nvars": 2**11 - 1,

34 "epsilon": 0.04,

35}

37# -----------------------------------------------------------------------------

38# Trajectories (reference solution)

39# -----------------------------------------------------------------------------

40uExact = solutionExact(tEnd, 1, pName, **pParams)

41x = np.linspace(-0.5, 0.5, 2**11 + 1)[1:-1]

43figName = f"{SCRIPT}_solution"

44plt.figure(figName)

45plt.plot(x, uExact[0, :], '-', label="\$u(0)\$")

46plt.plot(x, uExact[-1, :], '--', label="\$u(T)\$")

48plt.legend()

49plt.xlabel("\$x\$")

50plt.ylabel("Solution")

51plt.gcf().set_size_inches(12, 3)

52plt.tight_layout()

53plt.savefig(f"{PATH}/{figName}.pdf")

55# -----------------------------------------------------------------------------

56# %% Convergence and error VS cost plots

57# -----------------------------------------------------------------------------

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

59dtVals = tEnd / nStepsList

62def getError(uNum, uRef):

63 if uNum is None:

64 return np.inf

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

68def getCost(counters):

69 nNewton, nRHS, tComp = counters

70 return 2 * nNewton + nRHS

73minPrec = ["MIN-SR-NS", "MIN-SR-S", "MIN-SR-FLEX"]

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

76config = [

77 (*minPrec, "VDHS", "ESDIRK43", "LU"),

78]

81i = 0

82for qDeltaList in config:

83 figNameConv = f"{SCRIPT}_conv_{i}"

84 figNameCost = f"{SCRIPT}_cost_{i}"

85 i += 1

87 for qDelta, sym in zip(qDeltaList, symList):

88 try:

89 params = getParamsRK(qDelta)

90 except KeyError:

91 params = getParamsSDC(

93 )

95 errors = []

96 costs = []

98 for nSteps in nStepsList:

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

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

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 ls = '-' if qDelta.startswith("MIN-SR-") else "--"

113 plt.figure(figNameConv)

114 plt.loglog(dtVals, errors, sym + ls, label=qDelta)

115 data.storeAndCheck(f"{figNameConv}_{qDelta}", errors, atol=1e-4, rtol=1e-4)

117 plt.figure(figNameCost)

118 plt.loglog(costs, errors, sym + ls, label=qDelta)

120 for figName in [figNameConv, figNameCost]:

121 plt.figure(figName)

122 plt.gca().set(

123 xlabel="Cost" if "cost" in figName else r"\$\Delta {t}\$",

124 ylabel=r"\$L_2\$ error at \$T\$",

125 )

126 plt.legend()

127 plt.grid(True)

128 plt.tight_layout()

129 plt.savefig(f"{PATH}/{figName}.pdf")

131data.writeToJSON()