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

1#!/usr/bin/env python3 

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

3""" 

4Created on Fri Dec 22 17:17:14 2023 

5 

6Evaluate the nilpotency of diagonal preconditioners MIN-SR-S and MIN-SR-NS 

7with increasing number of nodes. 

8""" 

9 

10import numpy as np 

11import matplotlib.pyplot as plt 

12 

13from pySDC.core.sweeper import Sweeper 

14 

15quadType = "LOBATTO" 

16nodeType = "LEGENDRE" 

17 

18 

19def nilpotencyS(d, Q): 

20 if quadType in ['LOBATTO', 'RADAU-LEFT']: 

21 d = d[1:] 

22 Q = Q[1:, 1:] 

23 M = d.size 

24 D = np.diag(1 / d) 

25 K = np.eye(M) - D @ Q 

26 return np.linalg.norm(np.linalg.matrix_power(K, M), ord=np.inf) 

27 

28 

29def nilpotencyNS(d, Q): 

30 M = d.size 

31 D = np.diag(d) 

32 K = D - Q 

33 return np.linalg.norm(np.linalg.matrix_power(K, M), ord=np.inf) 

34 

35 

36nil_MIN_SR_S = [] 

37nil_MIN_SR_NS = [] 

38nNodes = range(2, 20) 

39for m in nNodes: 

40 s = Sweeper({"num_nodes": m, "quad_type": quadType, "node_type": nodeType}, None) 

41 Q = s.coll.Qmat[1:, 1:] 

42 nodes = s.coll.nodes 

43 

44 qDelta = s.get_Qdelta_implicit(qd_type="MIN-SR-S") 

45 d = np.diag(qDelta)[1:] 

46 nil_MIN_SR_S.append([nilpotencyS(d, Q), nilpotencyNS(d, Q)]) 

47 

48 qDelta = s.get_Qdelta_implicit(qd_type="MIN-SR-NS") 

49 d = np.diag(qDelta)[1:] 

50 nil_MIN_SR_NS.append([nilpotencyS(d, Q), nilpotencyNS(d, Q)]) 

51 

52nil_MIN_SR_NS = np.array(nil_MIN_SR_NS).T 

53nil_MIN_SR_S = np.array(nil_MIN_SR_S).T 

54 

55fig = plt.figure("nilpotency") 

56nNodes = np.array(nNodes, np.float128) 

57plt.semilogy(nNodes, nil_MIN_SR_NS[0], 'o-', label="MIN-SR-NS (nill. stiff)") 

58plt.semilogy(nNodes, nil_MIN_SR_NS[1], 'o--', label="MIN-SR-NS (nill. non-stiff)") 

59plt.semilogy(nNodes, nil_MIN_SR_S[0], 's-', label="MIN-SR-S (nill. stiff)") 

60plt.semilogy(nNodes, nil_MIN_SR_S[1], 's--', label="MIN-SR-S (nill. non-stiff)") 

61plt.legend() 

62plt.semilogy(nNodes, 14**nNodes * 1e-17, ':', c="gray") 

63plt.grid(True) 

64plt.xlabel("M") 

65plt.ylabel("nilpotency") 

66fig.set_size_inches(8.65, 5.33) 

67plt.tight_layout()