Coverage for pySDC/projects/DAE/sweepers/semiImplicitDAE.py: 98%
52 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-09-04 15:08 +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, level):
87 """Initialization routine"""
89 if 'QI' not in params:
90 params['QI'] = 'IE'
92 super().__init__(params, level)
94 msg = f"Quadrature type {self.params.quad_type} is not implemented yet. Use 'RADAU-RIGHT' instead!"
95 if self.coll.left_is_node:
96 raise ParameterError(msg)
98 self.QI = self.get_Qdelta_implicit(qd_type=self.params.QI)
100 def integrate(self):
101 r"""
102 Returns the solution by integrating its gradient (fundamental theorem of calculus) at each collocation node.
103 ``level.f`` stores the gradient of solution ``level.u``.
105 Returns
106 -------
107 me : list of lists
108 Integral of the gradient at each collocation node.
109 """
111 # get current level and problem description
112 L = self.level
113 P = L.prob
114 M = self.coll.num_nodes
116 me = []
117 for m in range(1, M + 1):
118 # new instance of dtype_u, initialize values with 0
119 me.append(P.dtype_u(P.init, val=0.0))
120 for j in range(1, M + 1):
121 me[-1].diff[:] += L.dt * self.coll.Qmat[m, j] * L.f[j].diff[:]
123 return me
125 def update_nodes(self):
126 r"""
127 Updates the values of solution ``u`` and their gradient stored in ``f``.
128 """
130 L = self.level
131 P = L.prob
133 # only if the level has been touched before
134 assert L.status.unlocked
135 M = self.coll.num_nodes
137 integral = self.integrate()
138 # build the rest of the known solution u_0 + del_t(Q - Q_del)U_k
139 for m in range(1, M + 1):
140 for j in range(1, m + 1):
141 integral[m - 1].diff[:] -= L.dt * self.QI[m, j] * L.f[j].diff[:]
142 integral[m - 1].diff[:] += L.u[0].diff
144 # do the sweep
145 for m in range(1, M + 1):
146 u_approx = P.dtype_u(integral[m - 1])
147 for j in range(1, m):
148 u_approx.diff[:] += L.dt * self.QI[m, j] * L.f[j].diff[:]
150 u0 = P.dtype_u(P.init)
151 u0.diff[:], u0.alg[:] = L.f[m].diff[:], L.u[m].alg[:]
152 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 # ---- update U' and z ----
155 L.f[m].diff[:] = u_new.diff[:]
156 L.u[m].alg[:] = u_new.alg[:]
158 # Update solution approximation
159 integral = self.integrate()
160 for m in range(M):
161 L.u[m + 1].diff[:] = L.u[0].diff[:] + integral[m].diff[:]
163 # indicate presence of new values at this level
164 L.status.updated = True
166 return None
168 @staticmethod
169 def F(du, P, factor, u_approx, t):
170 r"""
171 This function builds the implicit system to be solved for a semi-explicit
172 DAE of the form
174 .. math::
175 u' = f(u, z, t),
177 .. math::
178 0 = g(u, z, t)
180 Applying semi-implicit SDC yields the (non)-linear system to be solved
182 .. math::
183 0 = U_m^{k+1} - f(u_0 + \sum_{j=1}^M (q_{mj}-\tilde{q}_{mj}) U_j^k
184 + \sum_{j=1}^m \tilde{q}_{mj} U_j^{k+1}, z_m^{k+1}, \tau_m)
185 0 = g(u_0 + \sum_{j=1}^M (q_{mj}-\tilde{q}_{mj}) U_j^k
186 + \sum_{j=1}^m \tilde{q}_{mj} U_j^{k+1}, z_m^{k+1}, \tau_m)
188 which is solved for the derivative of u and the algebraic variable z.
190 Parameters
191 ----------
192 du : dtype_u
193 Unknowns of the system (derivative of solution u).
194 P : pySDC.projects.DAE.misc.problemDAE.ProblemDAE
195 Problem class.
196 factor : float
197 Abbrev. for the node-to-node stepsize.
198 u_approx : dtype_u
199 Approximation to numerical solution :math:`u`.
200 t : float
201 Current time :math:`t`.
203 Returns
204 -------
205 sys : dtype_f
206 System to be solved.
207 """
209 local_du_approx = P.dtype_f(du)
211 local_u_approx = P.dtype_u(u_approx)
212 local_u_approx.diff[:] += factor * local_du_approx.diff[:]
213 local_u_approx.alg[:] = local_du_approx.alg[:]
215 sys = P.eval_f(local_u_approx, local_du_approx, t)
216 return sys