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""" 

8 

9import numpy as np 

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

11 

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

13 

14# Script parameters 

15useRK = False 

16zoom = 20 

17reLims = -4.5 * zoom, 0.5 * zoom 

18imLims = -3.5 * zoom, 3.5 * zoom 

19nVals = 251 

20 

21 

22# RK parameters 

23rkScheme = "RK4" 

24 

25# Collocation parameters 

26nNodes = 4 

27nodeType = "LEGENDRE" 

28quadType = "RADAU-RIGHT" 

29 

30# SDC parameters 

31nSweeps = 6 

32qDeltaType = "VDHS" 

33collUpdate = False 

34 

35 

36# ----------------------------------------------------------------------------- 

37# Script execution 

38# ----------------------------------------------------------------------------- 

39 

40# Scheme instanciation 

41if useRK: # pragma: no cover 

42 params = getParamsRK(rkScheme) 

43else: 

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

45 

46# Problem instanciation 

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

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

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

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

51 

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

53stab = np.abs(uEnd) 

54 

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

56 

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

58if useRK: # pragma: no cover 

59 ax.set_title(rkScheme) 

60else: 

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

62 

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

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

65axs[1].tick_params( 

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

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

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

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

70 labelbottom=True, 

71) 

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

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

74axs[1].set_ylim(*imLims) 

75axs[1].set_aspect(0.2) 

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

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

78 

79plt.tight_layout()