Coverage for pySDC/projects/DAE/sweepers/rungeKuttaDAE.py: 100%
50 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +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):
80 super().__init__(params)
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 u_approx = prob.dtype_u(lvl.u[0])
146 for j in range(1, m + 1):
147 u_approx += lvl.dt * self.QI[m + 1, j] * lvl.f[j][:]
149 finit = lvl.f[m].flatten()
150 lvl.f[m + 1][:] = prob.solve_system(
151 FullyImplicitDAE.F,
152 u_approx,
153 lvl.dt * self.QI[m + 1, m + 1],
154 finit,
155 lvl.time + lvl.dt * self.coll.nodes[m + 1],
156 )
158 # Update numerical solution
159 integral = self.integrate()
160 for m in range(M):
161 lvl.u[m + 1][:] = lvl.u[0][:] + integral[m][:]
163 self.du_init = prob.dtype_f(lvl.f[-1])
165 lvl.status.updated = True
167 return None
170class BackwardEulerDAE(RungeKuttaDAE, BackwardEuler):
171 pass
174class TrapezoidalRuleDAE(RungeKuttaDAE, CrankNicolson):
175 pass
178class EDIRK4DAE(RungeKuttaDAE, EDIRK4):
179 pass
182class DIRK43_2DAE(RungeKuttaDAE, DIRK43_2):
183 pass