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

1import numpy as np 

2 

3from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit 

4 

5 

6class linearized_implicit_parallel(generic_implicit): 

7 """ 

8 Parallel sweeper using Newton for linearization 

9 """ 

10 

11 def __init__(self, params): 

12 """ 

13 Initialization routine for the custom sweeper 

14 

15 Args: 

16 params: parameters for the sweeper 

17 """ 

18 

19 if 'fixed_time_in_jacobian' not in params: 

20 params['fixed_time_in_jacobian'] = 0 

21 

22 # call parent's initialization routine 

23 super(linearized_implicit_parallel, self).__init__(params) 

24 

25 self.D, self.V = np.linalg.eig(self.QI[1:, 1:]) 

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

27 

28 def update_nodes(self): 

29 """ 

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

31 

32 Returns: 

33 None 

34 """ 

35 

36 # get current level and problem description 

37 L = self.level 

38 P = L.prob 

39 

40 # only if the level has been touched before 

41 assert L.status.unlocked 

42 

43 # get number of collocation nodes for easier access 

44 M = self.coll.num_nodes 

45 

46 # form Jacobian on each node 

47 dfdu = [] 

48 for m in range(M + 1): 

49 dfdu.append(P.eval_jacobian(L.u[m])) 

50 

51 # form collocation problem 

52 Gu = self.integrate() 

53 for m in range(M): 

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

55 

56 # transform collocation problem forward 

57 Guv = [] 

58 for m in range(M): 

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

60 for j in range(M): 

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

62 

63 # solve implicit system with Jacobians 

64 uv = [] 

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

66 uv.append( 

67 P.solve_system_jacobian( 

68 dfdu[m], Guv[m], L.dt * self.D[m], L.u[m + 1], L.time + L.dt * self.coll.nodes[m] 

69 ) 

70 ) 

71 

72 # transform solution backward 

73 for m in range(M): 

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

75 for j in range(M): 

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

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

78 

79 # evaluate f 

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

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

82 

83 # indicate presence of new values at this level 

84 L.status.updated = True 

85 

86 return None