Coverage for pySDC/projects/parallelSDC/linearized_implicit_fixed_parallel.py: 97%

1import numpy as np 

2 

3from pySDC.projects.parallelSDC.linearized_implicit_parallel import linearized_implicit_parallel 

4 

5 

6class linearized_implicit_fixed_parallel(linearized_implicit_parallel): 

7 """ 

8 Custom sweeper class, implements Sweeper.py 

9 

10 Generic implicit sweeper, expecting lower triangular matrix QI as input 

11 

12 Attributes: 

13 D: eigenvalues of the QI 

14 """ 

15 

16 def __init__(self, params): 

17 """ 

18 Initialization routine for the custom sweeper 

19 

20 Args: 

21 params: parameters for the sweeper 

22 """ 

23 

24 if 'fixed_time_in_jacobian' not in params: 

25 params['fixed_time_in_jacobian'] = 0 

26 

27 # call parent's initialization routine 

28 super(linearized_implicit_fixed_parallel, self).__init__(params) 

29 

30 assert self.params.fixed_time_in_jacobian in range(self.coll.num_nodes + 1), ( 

31 "ERROR: fixed_time_in_jacobian is too small or too large, got %s" % self.params.fixed_time_in_jacobian 

32 ) 

33 

34 self.D, self.V = np.linalg.eig(self.coll.Qmat[1:, 1:]) 

35 self.Vi = np.linalg.inv(self.V) 

36 

37 def update_nodes(self): 

38 """ 

39 Update the u- and f-values at the collocation nodes -> corresponds to a single sweep over all nodes 

40 

41 Returns: 

42 None 

43 """ 

44 

45 # get current level and problem description 

46 L = self.level 

47 P = L.prob 

48 

49 # only if the level has been touched before 

50 assert L.status.unlocked 

51 

52 # get number of collocation nodes for easier access 

53 M = self.coll.num_nodes 

54 

55 # form Jacobian at fixed time 

56 jtime = self.params.fixed_time_in_jacobian 

57 dfdu = P.eval_jacobian(L.u[jtime]) 

58 

59 # form collocation problem 

60 Gu = self.integrate() 

61 for m in range(M): 

62 Gu[m] -= L.u[m + 1] - L.u[0] 

63 

64 # transform collocation problem forward 

65 Guv = [] 

66 for m in range(M): 

67 Guv.append(P.dtype_u((P.init[0], P.init[1], np.dtype('complex128')), val=0.0 + 0.0j)) 

68 for j in range(M): 

69 Guv[m] += self.Vi[m, j] * Gu[j] 

70 

71 # solve implicit system with Jacobian (just this one, does not change with the nodes) 

72 uv = [] 

73 for m in range(M): # hell yeah, this is parallel!! 

74 uv.append( 

75 P.solve_system_jacobian(dfdu, Guv[m], L.dt * self.D[m], L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) 

76 ) 

77 

78 # transform soultion backward 

79 for m in range(M): 

80 tmp = P.dtype_u((P.init[0], P.init[1], np.dtype('complex128')), val=0.0 + 0.0j) 

81 for j in range(M): 

82 tmp += self.V[m, j] * uv[j] 

83 L.u[m + 1][:] += np.real(tmp) 

84 

85 # evaluate f 

86 for m in range(M): # hell yeah, this is parallel!! 

87 L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) 

88 

89 # indicate presence of new values at this level 

90 L.status.updated = True 

91 

92 return None