Coverage for pySDC/projects/DAE/sweepers/semiImplicitDAE.py: 98%

1from pySDC.core.errors import ParameterError 

2from pySDC.projects.DAE.sweepers.fullyImplicitDAE import FullyImplicitDAE 

3 

4 

5class SemiImplicitDAE(FullyImplicitDAE): 

6 r""" 

7 Custom sweeper class to implement SDC for solving semi-explicit DAEs of the form 

8 

9 .. math:: 

10 u' = f(u, z, t), 

11 

12 .. math:: 

13 0 = g(u, z, t) 

14 

15 with :math:`u(t), u'(t) \in\mathbb{R}^{N_d}` the differential variables and their derivates, 

16 algebraic variables :math:`z(t) \in\mathbb{R}^{N_a}`, :math:`f(u, z, t) \in \mathbb{R}^{N_d}`, 

17 and :math:`g(u, z, t) \in \mathbb{R}^{N_a}`. :math:`N = N_d + N_a` is the dimension of the whole 

18 system of DAEs. 

19 

20 It solves a collocation problem of the form 

21 

22 .. math:: 

23 U = f(\vec{U}_0 + \Delta t (\mathbf{Q} \otimes \mathbf{I}_{n_d}) \vec{U}, \vec{z}, \tau), 

24 

25 .. math:: 

26 0 = g(\vec{U}_0 + \Delta t (\mathbf{Q} \otimes \mathbf{I}_{n_d}) \vec{U}, \vec{z}, \tau), 

27 

28 where 

29 

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

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

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

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

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

35 - :math:`\vec{z}=(z_1,..,z_M) \in \mathbb{R}^{MN_a}` the vector of unknown algebraic variables 

36 :math:`z_m \approx z(\tau_m) \in \mathbb{R}^{N_a}`, 

37 - and identity matrix :math:`\mathbf{I}_{N_d} \in \mathbb{R}^{N_d \times N_d}`. 

38 

39 This sweeper treats the differential and the algebraic variables differently by only integrating the differential 

40 components. Solving the nonlinear system, :math:`{U,z}` are the unknowns. 

41 

42 The sweeper implementation is based on the ideas mentioned in the KDC publication [1]_. 

43 

44 Parameters 

45 ---------- 

46 params : dict 

47 Parameters passed to the sweeper. 

48 

49 Attributes 

50 ---------- 

51 QI : np.2darray 

52 Implicit Euler integration matrix. 

53 

54 References 

55 ---------- 

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

57 equation. J. Comput. Phys. Vol. 221 No. 2 (2007). 

58 

59 Note 

60 ---- 

61 The right-hand side of the problem DAE classes using this sweeper has to be exactly implemented in the way, the 

62 semi-explicit DAE is defined. Define :math:`\vec{x}=(y, z)^T`, :math:`F(\vec{x})=(f(\vec{x}), g(\vec{x}))`, and the 

63 matrix 

64 

65 .. math:: 

66 A = \begin{matrix} 

67 I & 0 \\ 

68 0 & 0 

69 \end{matrix} 

70 

71 then, the problem can be reformulated as 

72 

73 .. math:: 

74 A\vec{x}' = F(\vec{x}). 

75 

76 Then, setting :math:`F_{new}(\vec{x}, \vec{x}') = A\vec{x}' - F(\vec{x})` defines a DAE of fully-implicit form 

77 

78 .. math:: 

79 0 = F_{new}(\vec{x}, \vec{x}'). 

80 

81 Hence, the method ``eval_f`` of problem DAE classes of semi-explicit form implements the right-hand side in the way of 

82 returning :math:`F(\vec{x})`, whereas ``eval_f`` of problem classes of fully-implicit form return the right-hand side 

83 :math:`F_{new}(\vec{x}, \vec{x}')`. 

84 """ 

85 

86 def __init__(self, params): 

87 """Initialization routine""" 

88 

89 if 'QI' not in params: 

90 params['QI'] = 'IE' 

91 

92 # call parent's initialization routine 

93 super().__init__(params) 

94 

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

96 if self.coll.left_is_node: 

97 raise ParameterError(msg) 

98 

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

100 

101 def integrate(self): 

102 r""" 

103 Returns the solution by integrating its gradient (fundamental theorem of calculus) at each collocation node. 

104 ``level.f`` stores the gradient of solution ``level.u``. 

105 

106 Returns 

107 ------- 

108 me : list of lists 

109 Integral of the gradient at each collocation node. 

110 """ 

111 

112 # get current level and problem description 

113 L = self.level 

114 P = L.prob 

115 M = self.coll.num_nodes 

116 

117 me = [] 

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

119 # new instance of dtype_u, initialize values with 0 

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

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

122 me[-1].diff[:] += L.dt * self.coll.Qmat[m, j] * L.f[j].diff[:] 

123 

124 return me 

125 

126 def update_nodes(self): 

127 r""" 

128 Updates the values of solution ``u`` and their gradient stored in ``f``. 

129 """ 

130 

131 L = self.level 

132 P = L.prob 

133 

134 # only if the level has been touched before 

135 assert L.status.unlocked 

136 M = self.coll.num_nodes 

137 

138 integral = self.integrate() 

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

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

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

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

143 integral[m - 1].diff[:] += L.u[0].diff 

144 

145 # do the sweep 

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

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

148 for j in range(1, m): 

149 u_approx.diff[:] += L.dt * self.QI[m, j] * L.f[j].diff[:] 

150 

151 u0 = P.dtype_u(P.init) 

152 u0.diff[:], u0.alg[:] = L.f[m].diff[:], L.u[m].alg[:] 

153 u_new = P.solve_system(self.F, u_approx, L.dt * self.QI[m, m], u0, L.time + L.dt * self.coll.nodes[m - 1]) 

154 

155 # ---- update U' and z ---- 

156 L.f[m].diff[:] = u_new.diff[:] 

157 L.u[m].alg[:] = u_new.alg[:] 

158 

159 # Update solution approximation 

160 integral = self.integrate() 

161 for m in range(M): 

162 L.u[m + 1].diff[:] = L.u[0].diff[:] + integral[m].diff[:] 

163 

164 # indicate presence of new values at this level 

165 L.status.updated = True 

166 

167 return None 

168 

169 @staticmethod 

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

171 r""" 

172 This function builds the implicit system to be solved for a semi-explicit 

173 DAE of the form 

174 

175 .. math:: 

176 u' = f(u, z, t), 

177 

178 .. math:: 

179 0 = g(u, z, t) 

180 

181 Applying semi-implicit SDC yields the (non)-linear system to be solved 

182 

183 .. math:: 

184 0 = U_m^{k+1} - f(u_0 + \sum_{j=1}^M (q_{mj}-\tilde{q}_{mj}) U_j^k 

185 + \sum_{j=1}^m \tilde{q}_{mj} U_j^{k+1}, z_m^{k+1}, \tau_m) 

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

187 + \sum_{j=1}^m \tilde{q}_{mj} U_j^{k+1}, z_m^{k+1}, \tau_m) 

188 

189 which is solved for the derivative of u and the algebraic variable z. 

190 

191 Parameters 

192 ---------- 

193 du : dtype_u 

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

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

196 Problem class. 

197 factor : float 

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

199 u_approx : dtype_u 

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

201 t : float 

202 Current time :math:`t`. 

203 

204 Returns 

205 ------- 

206 sys : dtype_f 

207 System to be solved. 

208 """ 

209 

210 local_du_approx = P.dtype_f(du) 

211 

212 local_u_approx = P.dtype_u(u_approx) 

213 local_u_approx.diff[:] += factor * local_du_approx.diff[:] 

214 local_u_approx.alg[:] = local_du_approx.alg[:] 

215 

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

217 return sys