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 preconditionners MIN-SR-S and MIN-SR-NS 

7with increasing number of nodes. 

8""" 

9import numpy as np 

10import matplotlib.pyplot as plt 

11 

12from pySDC.core.Sweeper import sweeper 

13 

14quadType = "LOBATTO" 

15nodeType = "LEGENDRE" 

16 

17 

18def nilpotencyS(d, Q): 

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

20 d = d[1:] 

21 Q = Q[1:, 1:] 

22 M = d.size 

23 D = np.diag(1 / d) 

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

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

26 

27 

28def nilpotencyNS(d, Q): 

29 M = d.size 

30 D = np.diag(d) 

31 K = D - Q 

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

33 

34 

35nil_MIN_SR_S = [] 

36nil_MIN_SR_NS = [] 

37nNodes = range(2, 20) 

38for m in nNodes: 

39 s = sweeper({"num_nodes": m, "quad_type": quadType, "node_type": nodeType}) 

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

41 nodes = s.coll.nodes 

42 

43 qDelta = s.get_Qdelta_implicit(s.coll, qd_type="MIN-SR-S") 

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

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

46 

47 qDelta = s.get_Qdelta_implicit(s.coll, qd_type="MIN-SR-NS") 

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

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

50 

51nil_MIN_SR_NS = np.array(nil_MIN_SR_NS).T 

52nil_MIN_SR_S = np.array(nil_MIN_SR_S).T 

53 

54fig = plt.figure("nilpotency") 

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

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

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

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

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

60plt.legend() 

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

62plt.grid(True) 

63plt.xlabel("M") 

64plt.ylabel("nilpotency") 

65fig.set_size_inches(8.65, 5.33) 

66plt.tight_layout()