Coverage for pySDC/projects/parallelSDC_reloaded/stability.py: 100%

1#!/usr/bin/env python3 

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

3""" 

4Created on Tue Jan 9 10:12:09 2024 

5 

6Compute stability regions for SDC wit given parameters 

7""" 

8import numpy as np 

9from pySDC.projects.parallelSDC_reloaded.utils import getParamsRK, getParamsSDC, solutionSDC, plotStabContour, plt 

10 

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

12 

13# Script parameters 

14useRK = False 

15zoom = 20 

16reLims = -4.5 * zoom, 0.5 * zoom 

17imLims = -3.5 * zoom, 3.5 * zoom 

18nVals = 251 

19 

20 

21# RK parameters 

22rkScheme = "RK4" 

23 

24# Collocation parameters 

25nNodes = 4 

26nodeType = "LEGENDRE" 

27quadType = "RADAU-RIGHT" 

28 

29# SDC parameters 

30nSweeps = 6 

31qDeltaType = "VDHS" 

32collUpdate = False 

33 

34 

35# ----------------------------------------------------------------------------- 

36# Script execution 

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

38 

39# Scheme instanciation 

40if useRK: # pragma: no cover 

41 params = getParamsRK(rkScheme) 

42else: 

43 params = getParamsSDC(quadType, nNodes, qDeltaType, nSweeps, nodeType, collUpdate) 

44 

45# Problem instanciation 

46reVals = np.linspace(*reLims, num=nVals) 

47imVals = np.linspace(*imLims, num=nVals) 

48lambdas = reVals[None, :] + 1j * imVals[:, None] 

49uNum, counters, parallel = solutionSDC(1, 1, params, 'DAHLQUIST', lambdas=lambdas.ravel()) 

50 

51uEnd = uNum[-1, :].reshape(lambdas.shape) 

52stab = np.abs(uEnd) 

53 

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

55 

56ax = plotStabContour(reVals, imVals, stab, ax=axs[0]) 

57if useRK: # pragma: no cover 

58 ax.set_title(rkScheme) 

59else: 

60 ax.set_title(f"{qDeltaType}, K={nSweeps}") 

61 

62imStab = stab[:, np.argwhere(reVals == 0)].ravel() 

63axs[1].semilogx(imStab, imVals) 

64axs[1].tick_params( 

65 axis='x', # changes apply to the x-axis 

66 which='both', # both major and minor ticks are affected 

67 bottom=True, # ticks along the bottom edge are off 

68 top=False, # ticks along the top edge are off 

69 labelbottom=True, 

70) 

71axs[1].vlines(1, *imLims, linestyles='--', colors='black', linewidth=1) 

72axs[1].set_xlim([0.1, 10]) 

73axs[1].set_ylim(*imLims) 

74axs[1].set_aspect(0.2) 

75axs[1].set_xticks([0.1, 1, 10]) 

76axs[1].set_title("Imaginary axis") 

77 

78plt.tight_layout()