Coverage for pySDC/implementations/sweeper_classes/multi_implicit.py: 87%

1from pySDC.core.sweeper import Sweeper 

2 

3 

4class multi_implicit(Sweeper): 

5 """ 

6 SDC sweeper for right hand sides with two parts, each treated implicitly with its own solver and preconditioner. 

7 

8 First-order multi-implicit sweeper for two components 

9 

10 Attributes: 

11 Q1: implicit integration matrix for the first component 

12 Q2: implicit integration matrix for the second component 

13 """ 

14 

15 def __init__(self, params, level): 

16 """ 

17 Initialization routine for the custom sweeper 

18 

19 Args: 

20 params: parameters for the sweeper 

21 level (pySDC.Level.level): the level that uses this sweeper 

22 """ 

23 

24 # Default choice: implicit Euler 

25 if 'Q1' not in params: 

26 params['Q1'] = 'IE' 

27 if 'Q2' not in params: 

28 params['Q2'] = 'IE' 

29 

30 # call parent's initialization routine 

31 super(multi_implicit, self).__init__(params, level) 

32 

33 # Integration matrices 

34 self.Q1 = self.get_Qdelta_implicit(qd_type=self.params.Q1) 

35 self.Q2 = self.get_Qdelta_implicit(qd_type=self.params.Q2) 

36 

37 def integrate(self): 

38 """ 

39 Integrates the right-hand side (two components) 

40 

41 Returns: 

42 list of dtype_u: containing the integral as values 

43 """ 

44 

45 # get current level and problem description 

46 L = self.level 

47 P = L.prob 

48 M = self.coll.num_nodes 

49 

50 # sum the components once per node, not once per entry of Q 

51 F = [L.f[j].comp1 + L.f[j].comp2 for j in range(1, M + 1)] 

52 

53 me = [] 

54 

55 # integrate RHS over all collocation nodes 

56 for m in range(1, M + 1): 

57 # new instance of dtype_u, initialize values with 0 

58 me.append(P.dtype_u(P.init, val=0.0)) 

59 for j in range(1, M + 1): 

60 me[-1] += L.dt * self.coll.Qmat[m, j] * F[j - 1] 

61 

62 return me 

63 

64 def update_nodes(self): 

65 """ 

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

67 

68 Returns: 

69 None 

70 """ 

71 

72 # get current level and problem description 

73 L = self.level 

74 P = L.prob 

75 

76 # only if the level has been touched before 

77 assert L.status.unlocked 

78 

79 # get number of collocation nodes for easier access 

80 M = self.coll.num_nodes 

81 Q, Q1, Q2 = self.coll.Qmat, self.Q1, self.Q2 

82 

83 # the components of F(u^k), taken once: the sweep replaces L.f[m], so these stay the old values 

84 F1 = [L.f[j].comp1 for j in range(1, M + 1)] 

85 F2 = [L.f[j].comp2 for j in range(1, M + 1)] 

86 

87 # do the sweep 

88 for m in range(1, M + 1): 

89 t = L.time + L.dt * self.coll.nodes[m - 1] 

90 

91 # u0 + QF(u^k) - Q1F1(u^k) + tau, plus Q1F1(u^k+1) from the nodes already updated 

92 rhs = P.dtype_u(L.u[0]) 

93 for j in range(1, M + 1): 

94 rhs += L.dt * (Q[m, j] - Q1[m, j]) * F1[j - 1] 

95 rhs += L.dt * Q[m, j] * F2[j - 1] 

96 if L.tau[m - 1] is not None: 

97 rhs += L.tau[m - 1] 

98 for j in range(1, m): 

99 rhs += L.dt * Q1[m, j] * L.f[j].comp1 

100 

101 # implicit solve with prefactor stemming from Q1 

102 L.u[m] = P.solve_system_1(rhs, L.dt * Q1[m, m], L.u[m], t) 

103 

104 # substract Q2F2(u^k) and add Q2F2(u^k+1) 

105 rhs = P.dtype_u(L.u[m]) 

106 for j in range(1, M + 1): 

107 rhs -= L.dt * Q2[m, j] * F2[j - 1] 

108 for j in range(1, m): 

109 rhs += L.dt * Q2[m, j] * L.f[j].comp2 

110 

111 L.u[m] = P.solve_system_2(rhs, L.dt * Q2[m, m], L.u[m], t) # TODO: is L.u[m] a good guess? 

112 

113 # update function values 

114 L.f[m] = P.eval_f(L.u[m], t) 

115 

116 # indicate presence of new values at this level 

117 L.status.updated = True 

118 

119 return None 

120 

121 def compute_end_point(self): 

122 """ 

123 Compute u at the right point of the interval 

124 

125 The value uend computed here is a full evaluation of the Picard formulation unless do_full_update==False 

126 

127 Returns: 

128 None 

129 """ 

130 

131 # get current level and problem description 

132 L = self.level 

133 P = L.prob 

134 

135 # check if Mth node is equal to right point and do_coll_update is false, perform a simple copy 

136 if self.coll.right_is_node and not self.params.do_coll_update: 

137 # a copy is sufficient 

138 L.uend = P.dtype_u(L.u[-1]) 

139 else: 

140 # start with u0 and add integral over the full interval (using coll.weights) 

141 L.uend = P.dtype_u(L.u[0]) 

142 for m in range(self.coll.num_nodes): 

143 L.uend += L.dt * self.coll.weights[m] * (L.f[m + 1].comp1 + L.f[m + 1].comp2) 

144 # add up tau correction of the full interval (last entry) 

145 if L.tau[-1] is not None: 

146 L.uend += L.tau[-1] 

147 

148 return None