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): 

49 """Initialization routine""" 

50 

51 if 'QI' not in params: 

52 params['QI'] = 'IE' 

53 

54 # call parent's initialization routine 

55 super().__init__(params) 

56 

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

58 if self.coll.left_is_node: 

59 raise ParameterError(msg) 

60 

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

62 

63 def update_nodes(self): 

64 r""" 

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

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

67 """ 

68 

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 M = self.coll.num_nodes 

76 

77 # get QU^k where U = u' 

78 integral = self.integrate() 

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

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

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

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

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

84 

85 # do the sweep 

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

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

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

89 for j in range(1, m): 

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

91 

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

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

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

95 ) 

96 

97 # Update solution approximation 

98 integral = self.integrate() 

99 for m in range(M): 

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

101 

102 # indicate presence of new values at this level 

103 L.status.updated = True 

104 

105 return None 

106 

107 def predict(self): 

108 r""" 

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

110 

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

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

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

114 

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

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

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

118 """ 

119 # get current level and problem description 

120 L = self.level 

121 P = L.prob 

122 # set initial guess for gradient to zero 

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

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

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

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

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

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

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

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

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

132 # start with random initial guess 

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

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

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

136 else: 

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

138 

139 # indicate that this level is now ready for sweeps 

140 L.status.unlocked = True 

141 L.status.updated = True 

142 

143 def compute_residual(self, stage=None): 

144 r""" 

145 Uses the absolute value of the DAE system 

146 

147 .. math:: 

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

149 

150 for computing the residual in a chosen norm. 

151 

152 Parameters 

153 ---------- 

154 stage : str, optional 

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

156 """ 

157 

158 # get current level and problem description 

159 L = self.level 

160 P = L.prob 

161 

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

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

164 if stage in self.params.skip_residual_computation: 

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

166 return None 

167 

168 # compute the residual for each node 

169 res_norm = [] 

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

171 # use abs function from data type here 

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

173 

174 # find maximal residual over the nodes 

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

176 L.status.residual = max(res_norm) 

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

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

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

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

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

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

183 else: 

184 raise ParameterError( 

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

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

187 ) 

188 

189 # indicate that the residual has seen the new values 

190 L.status.updated = False 

191 

192 return None 

193 

194 def compute_end_point(self): 

195 """ 

196 Compute u at the right point of the interval. 

197 

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

199 

200 Returns: 

201 None 

202 """ 

203 

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

205 raise NotImplementedError() 

206 

207 super().compute_end_point() 

208 

209 @staticmethod 

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

211 r""" 

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

213 

214 .. math:: 

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

216 

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

218 

219 .. math:: 

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

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

222 

223 which is solved for the derivative of u. 

224 

225 Note 

226 ---- 

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

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

229 is the same. 

230 

231 Parameters 

232 ---------- 

233 du : dtype_u 

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

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

236 Problem class. 

237 factor : float 

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

239 u_approx : dtype_u 

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

241 t : float 

242 Current time :math:`t`. 

243 

244 Returns 

245 ------- 

246 sys : dtype_f 

247 System to be solved. 

248 """ 

249 

250 local_du_approx = P.dtype_f(du) 

251 

252 # build parameters to pass to implicit function 

253 local_u_approx = P.dtype_f(u_approx) 

254 

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

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

257 local_u_approx += factor * local_du_approx 

258 

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

260 return sys