Coverage for pySDC/projects/DAE/sweepers/rungeKuttaDAE.py: 100%
54 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +0000
1from pySDC.projects.DAE.sweepers.fullyImplicitDAE import FullyImplicitDAE
2from pySDC.implementations.sweeper_classes.Runge_Kutta import (
3 RungeKutta,
4 BackwardEuler,
5 CrankNicolson,
6 EDIRK4,
7 DIRK43_2,
8)
11class RungeKuttaDAE(RungeKutta):
12 r"""
13 Custom sweeper class to implement Runge-Kutta (RK) methods for general differential-algebraic equations (DAEs)
14 of the form
16 .. math::
17 0 = F(u, u', t).
19 RK methods for general DAEs have the form
21 .. math::
22 0 = F(u_0 + \Delta t \sum_{j=1}^M a_{i,j} U_j, U_m),
24 .. math::
25 u_M = u_0 + \Delta t \sum_{j=1}^M b_j U_j.
27 In pySDC, RK methods are implemented in the way that the coefficient matrix :math:`A` in the Butcher
28 tableau is a lower triangular matrix so that the stages are updated node-by-node. This class therefore only supports
29 RK methods with lower triangular matrices in the tableau.
31 Parameters
32 ----------
33 params : dict
34 Parameters passed to the sweeper.
36 Attributes
37 ----------
38 du_init : dtype_f
39 Stores the initial condition for each step.
41 Note
42 ----
43 When using a RK sweeper to simulate a problem make sure the DAE problem class has a ``du_exact`` method since RK methods need an initial
44 condition for :math:`u'(t_0)` as well.
46 In order to implement a new RK method for DAEs a new tableau can be added in ``pySDC.implementations.sweeper_classes.Runge_Kutta.py``.
47 For example, a new method called ``newRungeKuttaMethod`` with nodes :math:`c=(c_1, c_2, c_3)`, weights :math:`b=(b_1, b_2, b_3)` and
48 coefficient matrix
50 ..math::
51 \begin{eqnarray}
52 A = \begin{pmatrix}
53 a_{11} & 0 & 0 \\
54 a_{21} & a_{22} & 0 \\
55 a_{31} & a_{32} & & 0 \\
56 \end{pmatrix}
57 \end{eqnarray}
59 can be implemented as follows:
61 >>> class newRungeKuttaMethod(RungeKutta):
62 >>> nodes = np.array([c1, c2, c3])
63 >>> weights = np.array([b1, b2, b3])
64 >>> matrix = np.zeros((3, 3))
65 >>> matrix[0, 0] = a11
66 >>> matrix[1, :2] = [a21, a22]
67 >>> matrix[2, :] = [a31, a32, a33]
68 >>> ButcherTableauClass = ButcherTableau
70 The new class ``newRungeKuttaMethodDAE`` can then be used by defining the DAE class inheriting from both, this base class and class containing
71 the Butcher tableau:
73 >>> class newRungeKuttaMethodDAE(RungeKuttaDAE, newRungeKuttaMethod):
74 >>> pass
76 More details can be found [here](https://github.com/Parallel-in-Time/pySDC/blob/master/pySDC/implementations/sweeper_classes/Runge_Kutta.py).
77 """
79 def __init__(self, params, level):
80 super().__init__(params, level)
81 self.du_init = None
82 self.fully_initialized = False
84 def predict(self):
85 """
86 Predictor to fill values with zeros at nodes before first sweep.
87 """
89 # get current level and problem
90 lvl = self.level
91 prob = lvl.prob
93 if not self.fully_initialized:
94 self.du_init = prob.du_exact(lvl.time)
95 self.fully_initialized = True
97 lvl.f[0] = prob.dtype_f(self.du_init)
98 for m in range(1, self.coll.num_nodes + 1):
99 lvl.u[m] = prob.dtype_u(init=prob.init, val=0.0)
100 lvl.f[m] = prob.dtype_f(init=prob.init, val=0.0)
102 lvl.status.unlocked = True
103 lvl.status.updated = True
105 def integrate(self):
106 r"""
107 Returns the solution by integrating its gradient (fundamental theorem of calculus) at each collocation node.
108 ``level.f`` stores the gradient of solution ``level.u``.
110 Returns
111 -------
112 me : list of lists
113 Integral of the gradient at each collocation node.
114 """
116 # get current level and problem
117 lvl = self.level
118 prob = lvl.prob
120 # integrate RHS over all collocation nodes
121 me = []
122 for m in range(1, self.coll.num_nodes + 1):
123 # new instance of dtype_u, initialize values with 0
124 me.append(prob.dtype_u(prob.init, val=0.0))
125 for j in range(1, self.coll.num_nodes + 1):
126 me[-1] += lvl.dt * self.coll.Qmat[m, j] * lvl.f[j]
128 return me
130 def update_nodes(self):
131 r"""
132 Updates the values of solution ``u`` and their gradient stored in ``f``.
133 """
135 # get current level and problem description
136 lvl = self.level
137 prob = lvl.prob
139 # only if the level has been touched before
140 assert lvl.status.unlocked
141 assert lvl.status.sweep <= 1, "RK schemes are direct solvers. Please perform only 1 iteration!"
143 M = self.coll.num_nodes
144 for m in range(M):
145 # An explicit first stage (a_11 = 0, as in EDIRK4 and the trapezoidal rule) has no unknown to
146 # solve for: the stage value is just u_0, so the stage derivative is the du we came in with.
147 # Solving for it anyway is not merely wasteful, it is ill-posed. With a zero factor the
148 # algebraic part of F does not depend on the unknown at all, so the derivative of the algebraic
149 # variable is left completely undetermined and the solver returns whatever its trust region
150 # happened to wander to -- reporting success, since the residual is satisfied either way. That
151 # value feeds the following stages through u_approx and can pull them onto the wrong branch of
152 # the algebraic constraint, which for EDIRK4 on DiscontinuousTestDAE flipped the sign of z and
153 # cost six orders of magnitude of accuracy at one dt out of seven.
154 # Only the *first* stage can be short-circuited like this: a zero diagonal further down still
155 # has a non-empty sum below it, so its stage value is not u_0 and it needs solving some other way.
156 if self.QI[m + 1, m + 1] == 0:
157 assert m == 0, 'Only an explicit *first* stage is supported'
158 lvl.f[m + 1][:] = lvl.f[0][:]
159 continue
161 u_approx = prob.dtype_u(lvl.u[0])
162 for j in range(1, m + 1):
163 u_approx += lvl.dt * self.QI[m + 1, j] * lvl.f[j][:]
165 finit = lvl.f[m].flatten()
166 lvl.f[m + 1][:] = prob.solve_system(
167 FullyImplicitDAE.F,
168 u_approx,
169 lvl.dt * self.QI[m + 1, m + 1],
170 finit,
171 lvl.time + lvl.dt * self.coll.nodes[m + 1],
172 )
174 # Update numerical solution
175 integral = self.integrate()
176 for m in range(M):
177 lvl.u[m + 1][:] = lvl.u[0][:] + integral[m][:]
179 self.du_init = prob.dtype_f(lvl.f[-1])
181 lvl.status.updated = True
183 return None
186class BackwardEulerDAE(RungeKuttaDAE, BackwardEuler):
187 pass
190class TrapezoidalRuleDAE(RungeKuttaDAE, CrankNicolson):
191 pass
194class EDIRK4DAE(RungeKuttaDAE, EDIRK4):
195 pass
198class DIRK43_2DAE(RungeKuttaDAE, DIRK43_2):
199 pass