Coverage for pySDC/projects/DAE/sweepers/fully_implicit_DAE.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 

6from pySDC.projects.DAE.misc.DAEMesh import DAEMesh 

7 

8 

9class fully_implicit_DAE(generic_implicit): 

10 r""" 

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

12 of the form 

13 

14 .. math:: 

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

16 

17 It solves a collocation problem of the form 

18 

19 .. math:: 

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

21 

22 where 

23 

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

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

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

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

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

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

30 

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

32 

33 Parameters 

34 ---------- 

35 params : dict 

36 Parameters passed to the sweeper. 

37 

38 Attributes 

39 ---------- 

40 QI : np.2darray 

41 Implicit Euler integration matrix. 

42 

43 References 

44 ---------- 

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

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

47 """ 

48 

49 def __init__(self, params): 

50 """Initialization routine""" 

51 

52 if 'QI' not in params: 

53 params['QI'] = 'IE' 

54 

55 # call parent's initialization routine 

56 super().__init__(params) 

57 

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

59 if self.coll.left_is_node: 

60 raise ParameterError(msg) 

61 

62 self.QI = self.get_Qdelta_implicit(coll=self.coll, qd_type=self.params.QI) 

63 

64 def update_nodes(self): 

65 r""" 

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

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

68 """ 

69 

70 L = self.level 

71 P = L.prob 

72 

73 # only if the level has been touched before 

74 assert L.status.unlocked 

75 

76 M = self.coll.num_nodes 

77 

78 # get QU^k where U = u' 

79 integral = self.integrate() 

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

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

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

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

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

85 

86 # do the sweep 

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

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

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

90 for j in range(1, m): 

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

92 

93 # params contains U = u' 

94 def implSystem(params): 

95 """ 

96 Build implicit system to solve in order to find the unknowns. 

97 

98 Parameters 

99 ---------- 

100 params : dtype_u 

101 Unknowns of the system. 

102 

103 Returns 

104 ------- 

105 sys : 

106 System to be solved as implicit function. 

107 """ 

108 

109 params_mesh = P.dtype_f(params) 

110 

111 # build parameters to pass to implicit function 

112 local_u_approx = P.dtype_f(u_approx) 

113 

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

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

116 local_u_approx += L.dt * self.QI[m, m] * params_mesh 

117 

118 sys = P.eval_f(local_u_approx, params_mesh, L.time + L.dt * self.coll.nodes[m - 1]) 

119 return sys 

120 

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

122 L.f[m] = P.solve_system(implSystem, L.f[m], L.time + L.dt * self.coll.nodes[m - 1]) 

123 

124 # Update solution approximation 

125 integral = self.integrate() 

126 for m in range(M): 

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

128 

129 # indicate presence of new values at this level 

130 L.status.updated = True 

131 

132 return None 

133 

134 def predict(self): 

135 r""" 

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

137 

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

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

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

141 

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

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

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

145 """ 

146 # get current level and problem description 

147 L = self.level 

148 P = L.prob 

149 # set initial guess for gradient to zero 

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

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

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

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

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

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

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

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

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

159 # start with random initial guess 

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

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

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

163 else: 

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

165 

166 # indicate that this level is now ready for sweeps 

167 L.status.unlocked = True 

168 L.status.updated = True 

169 

170 def compute_residual(self, stage=None): 

171 r""" 

172 Uses the absolute value of the DAE system 

173 

174 .. math:: 

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

176 

177 for computing the residual. 

178 

179 Parameters 

180 ---------- 

181 stage : str, optional 

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

183 """ 

184 

185 # get current level and problem description 

186 L = self.level 

187 P = L.prob 

188 

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

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

191 if stage in self.params.skip_residual_computation: 

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

193 return None 

194 

195 # compute the residual for each node 

196 res_norm = [] 

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

198 # use abs function from data type here 

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

200 

201 # find maximal residual over the nodes 

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

203 L.status.residual = max(res_norm) 

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

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

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

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

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

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

210 else: 

211 raise ParameterError( 

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

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

214 ) 

215 

216 # indicate that the residual has seen the new values 

217 L.status.updated = False 

218 

219 return None 

220 

221 def compute_end_point(self): 

222 """ 

223 Compute u at the right point of the interval 

224 

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

226 

227 Returns: 

228 None 

229 """ 

230 

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

232 raise NotImplementedError() 

233 

234 super().compute_end_point()