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

## 52 statements

, created at 2024-09-09 14:59 +0000

1from pySDC.core.errors import ParameterError

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

5class SemiImplicitDAE(FullyImplicitDAE):

6 r"""

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

9 .. math::

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

12 .. math::

13 0 = g(u, z, t)

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.

20 It solves a collocation problem of the form

22 .. math::

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

25 .. math::

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

28 where

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}.

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.

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

44 Parameters

45 ----------

46 params : dict

47 Parameters passed to the sweeper.

49 Attributes

50 ----------

51 QI : np.2darray

52 Implicit Euler integration matrix.

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

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

65 .. math::

66 A = \begin{matrix}

67 I & 0 \\

68 0 & 0

69 \end{matrix}

71 then, the problem can be reformulated as

73 .. math::

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

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

78 .. math::

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

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 """

86 def __init__(self, params):

87 """Initialization routine"""

89 if 'QI' not in params:

90 params['QI'] = 'IE'

92 # call parent's initialization routine

93 super().__init__(params)

96 if self.coll.left_is_node:

97 raise ParameterError(msg)

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

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.

106 Returns

107 -------

108 me : list of lists

109 Integral of the gradient at each collocation node.

110 """

112 # get current level and problem description

113 L = self.level

114 P = L.prob

115 M = self.coll.num_nodes

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[:]

124 return me

126 def update_nodes(self):

127 r"""

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

129 """

131 L = self.level

132 P = L.prob

134 # only if the level has been touched before

135 assert L.status.unlocked

136 M = self.coll.num_nodes

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

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[:]

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

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

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

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

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[:]

164 # indicate presence of new values at this level

165 L.status.updated = True

167 return None

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

175 .. math::

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

178 .. math::

179 0 = g(u, z, t)

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

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)

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

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.

204 Returns

205 -------

206 sys : dtype_f

207 System to be solved.

208 """

210 local_du_approx = P.dtype_f(du)

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[:]

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

217 return sys