30 statements

, created at 2024-09-09 14:59 +0000

1#!/usr/bin/env python3

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

3"""

4Created on Thu Dec 7 13:12:40 2023

6Setup script for the ProtheroRobinson (linear and non-linear) problem,

7using the autonomous formulation

8"""

9import numpy as np

10import matplotlib.pyplot as plt

12from pySDC.projects.parallelSDC_reloaded.utils import solutionExact, getParamsRK, solutionSDC

14script = __file__.split('/')[-1].split('.')[0]

16tEnd = 2 * np.pi

17nSteps = 2

18epsilon = 1e-3

19rkScheme = "ESDIRK43"

21tVals = np.linspace(0, tEnd, nSteps + 1)

22tValsFine = np.linspace(0, tEnd, nSteps * 10 + 1)

24for pType in ["linear", "nonlinear"]:

25 probName = "PROTHERO-ROBINSON-A"

26 if pType == "nonlinear":

27 probName += "-NL"

29 print(f"Computing {pType} ODE solution")

30 uExact = solutionExact(tEnd, nSteps, probName, epsilon=epsilon)

32 params = getParamsRK(rkScheme)

33 uNum, counters, parallel = solutionSDC(tEnd, nSteps, params, probName, epsilon=epsilon)

34 uNumFine, counters, parallel = solutionSDC(tEnd, nSteps * 10, params, probName, epsilon=epsilon)

36 figName = f"{script}_{pType}"

37 plt.figure(figName)

38 plt.plot(tVals, uExact[:, 0], '-', label="Exact")

39 plt.plot(tVals, uNum[:, 0], '--', label="Numerical")

40 plt.plot(tValsFine, uNumFine[:, 0], '-.', label="Numerical (fine)")

43for figName in [f"{script}_linear", f"{script}_nonlinear"]:

44 plt.figure(figName)

45 plt.legend()

46 plt.xlabel("time")

47 plt.ylabel("solution")

48 plt.tight_layout()