Coverage for pySDC/implementations/sweeper_classes/imex_1st_order.py: 100%

1import numpy as np 

2 

3from pySDC.core.sweeper import Sweeper 

4 

5 

6class imex_1st_order(Sweeper): 

7 """ 

8 SDC sweeper for right hand sides split into an implicitly and an explicitly treated part (IMEX-SDC). 

9 

10 First-order IMEX sweeper using implicit/explicit Euler as base integrator 

11 

12 Attributes: 

13 QI: implicit Euler integration matrix 

14 QE: explicit Euler integration matrix 

15 """ 

16 

17 def __init__(self, params, level): 

18 """ 

19 Initialization routine for the custom sweeper 

20 

21 Args: 

22 params: parameters for the sweeper 

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

24 """ 

25 

26 if 'QI' not in params: 

27 params['QI'] = 'IE' 

28 if 'QE' not in params: 

29 params['QE'] = 'EE' 

30 

31 super().__init__(params, level) 

32 

33 # IMEX integration matrices 

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

35 self.QE = self.get_Qdelta_explicit(qd_type=self.params.QE) 

36 

37 def integrate(self): 

38 """ 

39 Integrates the right-hand side (here impl + expl) 

40 

41 Returns: 

42 list of dtype_u: containing the integral as values 

43 """ 

44 

45 L = self.level 

46 P = L.prob 

47 

48 # sum the two parts once per node, not once per entry of Q 

49 F = [L.f[j].impl + L.f[j].expl for j in range(1, self.coll.num_nodes + 1)] 

50 

51 me = [] 

52 # integrate RHS over all collocation nodes 

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

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

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

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

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 # the terms which are known already (from the previous iteration) are u0 + QF(u^k) - QIFI(u^k) - QEFE(u^k) + tau, 

79 # built with Q - QI and Q - QE in one pass over the parts of F(u^k), kept here because the sweep replaces L.f[m] 

80 FI = [L.f[j].impl for j in range(1, M + 1)] 

81 FE = [L.f[j].expl for j in range(1, M + 1)] 

82 diff_weights_impl = self.coll.Qmat - self.QI 

83 diff_weights_expl = self.coll.Qmat - self.QE 

84 

85 # do the sweep 

86 for m in range(0, M): 

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

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

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

90 rhs += L.dt * diff_weights_impl[m + 1, j] * FI[j - 1] 

91 rhs += L.dt * diff_weights_expl[m + 1, j] * FE[j - 1] 

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

93 rhs += L.tau[m] 

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

95 rhs += L.dt * self.QI[m + 1, j] * L.f[j].impl 

96 rhs += L.dt * self.QE[m + 1, j] * L.f[j].expl 

97 

98 # implicit solve with prefactor stemming from QI 

99 L.u[m + 1] = P.solve_system( 

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

101 ) 

102 

103 # update function values 

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

105 

106 # indicate presence of new values at this level 

107 L.status.updated = True 

108 

109 return None 

110 

111 def compute_end_point(self): 

112 """ 

113 Compute u at the right point of the interval 

114 

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

116 

117 Returns: 

118 None 

119 """ 

120 

121 # get current level and problem description 

122 L = self.level 

123 P = L.prob 

124 

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

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

127 # a copy is sufficient 

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

129 else: 

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

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

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

133 L.uend += L.dt * self.coll.weights[m] * (L.f[m + 1].impl + L.f[m + 1].expl) 

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

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

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

137 

138 return None 

139 

140 def get_sweeper_mats(self): 

141 """ 

142 Returns the three matrices Q, QI, QE which define the sweeper. 

143 

144 The first row and column, corresponding to the left starting value, are removed to correspond to the notation 

145 introduced in Ruprecht & Speck, Spectral deferred corrections with fast-wave slow-wave splitting, 2016 

146 """ 

147 QE = self.QE[1:, 1:] 

148 QI = self.QI[1:, 1:] 

149 Q = self.coll.Qmat[1:, 1:] 

150 return QE, QI, Q 

151 

152 def get_scalar_problems_sweeper_mats(self, lambdas=None): 

153 """ 

154 This function returns the corresponding matrices of an IMEX-SDC sweep in matrix formulation. 

155 

156 See Ruprecht & Speck, Spectral deferred corrections with fast-wave slow-wave splitting, 2016 for the derivation. 

157 

158 Args: 

159 lambdas (numpy.ndarray): the first entry in lambdas is lambda_fast, the second is lambda_slow. 

160 """ 

161 QE, QI, Q = self.get_sweeper_mats() 

162 if lambdas is None: 

163 pass 

164 # should use lambdas from attached problem and make sure it is a scalar IMEX 

165 raise NotImplementedError("At the moment, the values for lambda have to be provided") 

166 else: 

167 lambda_fast = lambdas[0] 

168 lambda_slow = lambdas[1] 

169 nnodes = self.coll.num_nodes 

170 dt = self.level.dt 

171 LHS = np.eye(nnodes) - dt * (lambda_fast * QI + lambda_slow * QE) 

172 RHS = dt * ((lambda_fast + lambda_slow) * Q - (lambda_fast * QI + lambda_slow * QE)) 

173 return LHS, RHS 

174 

175 def get_scalar_problems_manysweep_mat(self, nsweeps, lambdas=None): 

176 """ 

177 For a scalar problem, K sweeps of IMEX-SDC can be written in matrix form. 

178 

179 Args: 

180 nsweeps (int): number of sweeps 

181 lambdas (numpy.ndarray): the first entry in lambdas is lambda_fast, the second is lambda_slow. 

182 """ 

183 LHS, RHS = self.get_scalar_problems_sweeper_mats(lambdas=lambdas) 

184 Pinv = np.linalg.inv(LHS) 

185 Mat_sweep = np.linalg.matrix_power(Pinv.dot(RHS), nsweeps) 

186 for k in range(0, nsweeps): 

187 Mat_sweep += np.linalg.matrix_power(Pinv.dot(RHS), k).dot(Pinv) 

188 return Mat_sweep