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

1from pySDC.core.Errors import ParameterError 

2from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE 

3 

4 

5class SemiImplicitDAE(fully_implicit_DAE): 

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(coll=self.coll, 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 def implSystem(unknowns): 

152 """ 

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

154 

155 Parameters 

156 ---------- 

157 unknowns : dtype_u 

158 Unknowns of the system. 

159 

160 Returns 

161 ------- 

162 sys : 

163 System to be solved as implicit function. 

164 """ 

165 

166 unknowns_mesh = P.dtype_f(unknowns) 

167 

168 local_u_approx = P.dtype_u(u_approx) 

169 local_u_approx.diff[:] += L.dt * self.QI[m, m] * unknowns_mesh.diff[:] 

170 local_u_approx.alg[:] = unknowns_mesh.alg[:] 

171 

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

173 return sys 

174 

175 u0 = P.dtype_u(P.init) 

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

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

178 # ---- update U' and z ---- 

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

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

181 

182 # Update solution approximation 

183 integral = self.integrate() 

184 for m in range(M): 

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

186 

187 # indicate presence of new values at this level 

188 L.status.updated = True 

189 

190 return None