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

1from pySDC.core.sweeper import Sweeper 

2 

3 

4class multi_implicit(Sweeper): 

5 """ 

6 Custom sweeper class, implements Sweeper.py 

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 

49 me = [] 

50 

51 # integrate RHS over all collocation nodes 

52 for m in range(1, self.coll.num_nodes + 1): 

53 # new instance of dtype_u, initialize values with 0 

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

55 for j in range(1, self.coll.num_nodes + 1): 

56 me[-1] += L.dt * self.coll.Qmat[m, j] * (L.f[j].comp1 + L.f[j].comp2) 

57 

58 return me 

59 

60 def update_nodes(self): 

61 """ 

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

63 

64 Returns: 

65 None 

66 """ 

67 

68 # get current level and problem description 

69 L = self.level 

70 P = L.prob 

71 

72 # only if the level has been touched before 

73 assert L.status.unlocked 

74 

75 # get number of collocation nodes for easier access 

76 M = self.coll.num_nodes 

77 

78 # gather all terms which are known already (e.g. from the previous iteration) 

79 

80 # get QF(u^k) 

81 integral = self.integrate() 

82 for m in range(M): 

83 # subtract Q1F1(u^k)_m 

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

85 integral[m] -= L.dt * self.Q1[m + 1, j] * L.f[j].comp1 

86 # add initial value 

87 integral[m] += L.u[0] 

88 # add tau if associated 

89 if L.tau[m] is not None: 

90 integral[m] += L.tau[m] 

91 

92 # store Q2F2(u^k) for later usage 

93 Q2int = [] 

94 for m in range(M): 

95 Q2int.append(P.dtype_u(P.init, val=0.0)) 

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

97 Q2int[-1] += L.dt * self.Q2[m + 1, j] * L.f[j].comp2 

98 

99 # do the sweep 

100 for m in range(0, M): 

101 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

102 rhs = P.dtype_u(integral[m]) 

103 for j in range(1, m + 1): 

104 rhs += L.dt * self.Q1[m + 1, j] * L.f[j].comp1 

105 

106 # implicit solve with prefactor stemming from Q1 

107 L.u[m + 1] = P.solve_system_1( 

108 rhs, L.dt * self.Q1[m + 1, m + 1], L.u[m + 1], L.time + L.dt * self.coll.nodes[m] 

109 ) 

110 

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

112 rhs = L.u[m + 1] - Q2int[m] 

113 for j in range(1, m + 1): 

114 rhs += L.dt * self.Q2[m + 1, j] * L.f[j].comp2 

115 

116 L.u[m + 1] = P.solve_system_2( 

117 rhs, 

118 L.dt * self.Q2[m + 1, m + 1], 

119 L.u[m + 1], # TODO: is this a good guess? 

120 L.time + L.dt * self.coll.nodes[m], 

121 ) 

122 

123 # update function values 

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

125 

126 # indicate presence of new values at this level 

127 L.status.updated = True 

128 

129 return None 

130 

131 def compute_end_point(self): 

132 """ 

133 Compute u at the right point of the interval 

134 

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

136 

137 Returns: 

138 None 

139 """ 

140 

141 # get current level and problem description 

142 L = self.level 

143 P = L.prob 

144 

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

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

147 # a copy is sufficient 

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

149 else: 

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

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

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

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

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

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

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

157 

158 return None