Coverage for pySDC/projects/DAE/sweepers/fullyImplicitDAE.py: 97%

1import numpy as np 

2from scipy import optimize 

3 

4from pySDC.core.errors import ParameterError 

5from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit 

6 

7 

8class FullyImplicitDAE(generic_implicit): 

9 r""" 

10 Custom sweeper class to implement the fully-implicit SDC for solving DAEs. It solves fully-implicit DAE problems 

11 of the form 

12 

13 .. math:: 

14 F(t, u, u') = 0. 

15 

16 It solves a collocation problem of the form 

17 

18 .. math:: 

19 F(\tau, \vec{U}_0 + \Delta t (\mathbf{Q} \otimes \mathbf{I}_n) \vec{U}, \vec{U}) = 0, 

20 

21 where 

22 

23 - :math:`\tau=(\tau_1,..,\tau_M) in \mathbb{R}^M` the vector of collocation nodes, 

24 - :math:`\vec{U}_0 = (u_0,..,u_0) \in \mathbb{R}^{Mn}` the vector of initial condition spread to each node, 

25 - spectral integration matrix :math:`\mathbf{Q} \in \mathbb{R}^{M \times M}`, 

26 - :math:`\vec{U}=(U_1,..,U_M) \in \mathbb{R}^{Mn}` the vector of unknown derivatives 

27 :math:`U_m \approx U(\tau_m) = u'(\tau_m) \in \mathbb{R}^n`, 

28 - and identity matrix :math:`\mathbf{I}_n \in \mathbb{R}^{n \times n}`. 

29 

30 The construction of this sweeper is based on the concepts outlined in [1]_. 

31 

32 Parameters 

33 ---------- 

34 params : dict 

35 Parameters passed to the sweeper. 

36 

37 Attributes 

38 ---------- 

39 QI : np.2darray 

40 Implicit Euler integration matrix. 

41 

42 References 

43 ---------- 

44 .. [1] J. Huang, J. Jun, M. L. Minion. Arbitrary order Krylov deferred correction methods for differential algebraic equation. 

45 J. Comput. Phys. Vol. 221 No. 2 (2007). 

46 """ 

47 

48 def __init__(self, params, level): 

49 """Initialization routine""" 

50 

51 if 'QI' not in params: 

52 params['QI'] = 'IE' 

53 

54 super().__init__(params, level) 

55 

56 msg = f"Quadrature type {self.params.quad_type} is not implemented yet. Use 'RADAU-RIGHT' instead!" 

57 if self.coll.left_is_node: 

58 raise ParameterError(msg) 

59 

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

61 

62 def update_nodes(self): 

63 r""" 

64 Updates values of ``u`` and ``f`` at collocation nodes. This correspond to a single iteration of the 

65 preconditioned Richardson iteration in **"ordinary"** SDC. 

66 """ 

67 

68 L = self.level 

69 P = L.prob 

70 

71 # only if the level has been touched before 

72 assert L.status.unlocked 

73 

74 M = self.coll.num_nodes 

75 

76 # get QU^k where U = u' 

77 integral = self.integrate() 

78 # build the rest of the known solution u_0 + del_t(Q - Q_del)U_k 

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

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

81 integral[m - 1] -= L.dt * self.QI[m, j] * L.f[j] 

82 integral[m - 1] += L.u[0] 

83 

84 # do the sweep 

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

86 # add the known components from current sweep del_t*Q_del*U_k+1 

87 u_approx = P.dtype_u(integral[m - 1]) 

88 for j in range(1, m): 

89 u_approx += L.dt * self.QI[m, j] * L.f[j] 

90 

91 # update gradient (recall L.f is being used to store the gradient) 

92 L.f[m] = P.solve_system( 

93 self.F, u_approx, L.dt * self.QI[m, m], L.f[m], L.time + L.dt * self.coll.nodes[m - 1] 

94 ) 

95 

96 # Update solution approximation 

97 integral = self.integrate() 

98 for m in range(M): 

99 L.u[m + 1] = L.u[0] + integral[m] 

100 

101 # indicate presence of new values at this level 

102 L.status.updated = True 

103 

104 return None 

105 

106 def predict(self): 

107 r""" 

108 Predictor to fill values at nodes before first sweep. It can decide whether the 

109 

110 - initial condition is spread to each node ('initial_guess' = 'spread'), 

111 - zero values are spread to each node ('initial_guess' = 'zero'), 

112 - or random values are spread to each collocation node ('initial_guess' = 'random'). 

113 

114 Default prediction for the sweepers, only copies the values to all collocation nodes. This function 

115 overrides the base implementation by always initialising ``level.f`` to zero. This is necessary since 

116 ``level.f`` stores the solution derivative in the fully implicit case, which is not initially known. 

117 """ 

118 # get current level and problem description 

119 L = self.level 

120 P = L.prob 

121 # set initial guess for gradient to zero 

122 L.f[0] = P.dtype_f(init=P.init, val=0.0) 

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

124 # copy u[0] to all collocation nodes and set f (the gradient) to zero 

125 if self.params.initial_guess == 'spread': 

126 L.u[m] = P.dtype_u(L.u[0]) 

127 L.f[m] = P.dtype_f(init=P.init, val=0.0) 

128 elif self.params.initial_guess == 'zero': 

129 L.u[m] = P.dtype_u(init=P.init, val=0.0) 

130 L.f[m] = P.dtype_f(init=P.init, val=0.0) 

131 # start with random initial guess 

132 elif self.params.initial_guess == 'random': 

133 L.u[m] = P.dtype_u(init=P.init, val=np.random.rand(1)[0]) 

134 L.f[m] = P.dtype_f(init=P.init, val=np.random.rand(1)[0]) 

135 else: 

136 raise ParameterError(f'initial_guess option {self.params.initial_guess} not implemented') 

137 

138 # indicate that this level is now ready for sweeps 

139 L.status.unlocked = True 

140 L.status.updated = True 

141 

142 def compute_residual(self, stage=None): 

143 r""" 

144 Uses the absolute value of the DAE system 

145 

146 .. math:: 

147 ||F(t, u, u')|| 

148 

149 for computing the residual in a chosen norm. 

150 

151 Parameters 

152 ---------- 

153 stage : str, optional 

154 The current stage of the step the level belongs to. 

155 """ 

156 

157 # get current level and problem description 

158 L = self.level 

159 P = L.prob 

160 

161 # Check if we want to skip the residual computation to gain performance 

162 # Keep in mind that skipping any residual computation is likely to give incorrect outputs of the residual! 

163 if stage in self.params.skip_residual_computation: 

164 L.status.residual = 0.0 if L.status.residual is None else L.status.residual 

165 return None 

166 

167 # compute the residual for each node 

168 res_norm = [] 

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

170 # use abs function from data type here 

171 res_norm.append(abs(P.eval_f(L.u[m + 1], L.f[m + 1], L.time + L.dt * self.coll.nodes[m]))) 

172 

173 # find maximal residual over the nodes 

174 if L.params.residual_type == 'full_abs': 

175 L.status.residual = max(res_norm) 

176 elif L.params.residual_type == 'last_abs': 

177 L.status.residual = res_norm[-1] 

178 elif L.params.residual_type == 'full_rel': 

179 L.status.residual = max(res_norm) / abs(L.u[0]) 

180 elif L.params.residual_type == 'last_rel': 

181 L.status.residual = res_norm[-1] / abs(L.u[0]) 

182 else: 

183 raise ParameterError( 

184 f'residual_type = {L.params.residual_type} not implemented, choose ' 

185 f'full_abs, last_abs, full_rel or last_rel instead' 

186 ) 

187 

188 # indicate that the residual has seen the new values 

189 L.status.updated = False 

190 

191 return None 

192 

193 def compute_end_point(self): 

194 """ 

195 Compute u at the right point of the interval. 

196 

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

198 

199 Returns: 

200 None 

201 """ 

202 

203 if not self.coll.right_is_node or self.params.do_coll_update: 

204 raise NotImplementedError() 

205 

206 super().compute_end_point() 

207 

208 @staticmethod 

209 def F(du, P, factor, u_approx, t): 

210 r""" 

211 This function builds the implicit system to be solved for a DAE of the form 

212 

213 .. math:: 

214 0 = F(u, u', t) 

215 

216 Applying SDC yields the (non)-linear system to be solved 

217 

218 .. math:: 

219 0 = F(u_0 + \sum_{j=1}^M (q_{mj}-\tilde{q}_{mj}) U_j 

220 + \sum_{j=1}^m \tilde{q}_{mj} U_j, U_m, \tau_m), 

221 

222 which is solved for the derivative of u. 

223 

224 Note 

225 ---- 

226 This function is also used for Runge-Kutta methods since only 

227 the argument ``u_approx`` differs. However, the implicit system to be solved 

228 is the same. 

229 

230 Parameters 

231 ---------- 

232 du : dtype_u 

233 Unknowns of the system (derivative of solution u). 

234 P : pySDC.projects.DAE.misc.problemDAE.ProblemDAE 

235 Problem class. 

236 factor : float 

237 Abbrev. for the node-to-node stepsize. 

238 u_approx : dtype_u 

239 Approximation to numerical solution :math:`u`. 

240 t : float 

241 Current time :math:`t`. 

242 

243 Returns 

244 ------- 

245 sys : dtype_f 

246 System to be solved. 

247 """ 

248 

249 local_du_approx = P.dtype_f(du) 

250 

251 # build parameters to pass to implicit function 

252 local_u_approx = P.dtype_f(u_approx) 

253 

254 # note that derivatives of algebraic variables are taken into account here too 

255 # these do not directly affect the output of eval_f but rather indirectly via QI 

256 local_u_approx += factor * local_du_approx 

257 

258 sys = P.eval_f(local_u_approx, local_du_approx, t) 

259 return sys