## 50 statements

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

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

2from pySDC.implementations.sweeper_classes.Runge_Kutta import (

3 RungeKutta,

4 BackwardEuler,

5 CrankNicholson,

6 EDIRK4,

7 DIRK43_2,

8)

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:

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

171 pass