# Coverage for pySDC/projects/DAE/sweepers/fullyImplicitDAE.py: 97%

## 77 statements

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

1import numpy as np

2from scipy import optimize

4from pySDC.core.errors import ParameterError

5from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit

8class FullyImplicitDAE(generic_implicit):

9 r"""

10 Custom sweeper class to implement the fully-implicit SDC for solving DAEs. It solves fully-implicit DAE problems

11 of the form

13 .. math::

14 F(t, u, u') = 0.

16 It solves a collocation problem of the form

18 .. math::

19 F(\tau, \vec{U}_0 + \Delta t (\mathbf{Q} \otimes \mathbf{I}_n) \vec{U}, \vec{U}) = 0,

21 where

23 - :math:\tau=(\tau_1,..,\tau_M) in \mathbb{R}^M the vector of collocation nodes,

24 - :math:\vec{U}_0 = (u_0,..,u_0) \in \mathbb{R}^{Mn} the vector of initial condition spread to each node,

25 - spectral integration matrix :math:\mathbf{Q} \in \mathbb{R}^{M \times M},

26 - :math:\vec{U}=(U_1,..,U_M) \in \mathbb{R}^{Mn} the vector of unknown derivatives

27 :math:U_m \approx U(\tau_m) = u'(\tau_m) \in \mathbb{R}^n,

28 - and identity matrix :math:\mathbf{I}_n \in \mathbb{R}^{n \times n}.

30 The construction of this sweeper is based on the concepts outlined in [1]_.

32 Parameters

33 ----------

34 params : dict

35 Parameters passed to the sweeper.

37 Attributes

38 ----------

39 QI : np.2darray

40 Implicit Euler integration matrix.

42 References

43 ----------

44 .. [1] J. Huang, J. Jun, M. L. Minion. Arbitrary order Krylov deferred correction methods for differential algebraic equation.

45 J. Comput. Phys. Vol. 221 No. 2 (2007).

46 """

48 def __init__(self, params):

49 """Initialization routine"""

51 if 'QI' not in params:

52 params['QI'] = 'IE'

54 # call parent's initialization routine

55 super().__init__(params)

58 if self.coll.left_is_node:

59 raise ParameterError(msg)

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

63 def update_nodes(self):

64 r"""

65 Updates values of u and f at collocation nodes. This correspond to a single iteration of the

66 preconditioned Richardson iteration in **"ordinary"** SDC.

67 """

69 L = self.level

70 P = L.prob

72 # only if the level has been touched before

73 assert L.status.unlocked

75 M = self.coll.num_nodes

77 # get QU^k where U = u'

78 integral = self.integrate()

79 # build the rest of the known solution u_0 + del_t(Q - Q_del)U_k

80 for m in range(1, M + 1):

81 for j in range(1, M + 1):

82 integral[m - 1] -= L.dt * self.QI[m, j] * L.f[j]

83 integral[m - 1] += L.u[0]

85 # do the sweep

86 for m in range(1, M + 1):

87 # add the known components from current sweep del_t*Q_del*U_k+1

88 u_approx = P.dtype_u(integral[m - 1])

89 for j in range(1, m):

90 u_approx += L.dt * self.QI[m, j] * L.f[j]

92 # update gradient (recall L.f is being used to store the gradient)

93 L.f[m] = P.solve_system(

94 self.F, u_approx, L.dt * self.QI[m, m], L.f[m], L.time + L.dt * self.coll.nodes[m - 1]

95 )

97 # Update solution approximation

98 integral = self.integrate()

99 for m in range(M):

100 L.u[m + 1] = L.u[0] + integral[m]

102 # indicate presence of new values at this level

103 L.status.updated = True

105 return None

107 def predict(self):

108 r"""

109 Predictor to fill values at nodes before first sweep. It can decide whether the

111 - initial condition is spread to each node ('initial_guess' = 'spread'),

112 - zero values are spread to each node ('initial_guess' = 'zero'),

113 - or random values are spread to each collocation node ('initial_guess' = 'random').

115 Default prediction for the sweepers, only copies the values to all collocation nodes. This function

116 overrides the base implementation by always initialising level.f to zero. This is necessary since

117 level.f stores the solution derivative in the fully implicit case, which is not initially known.

118 """

119 # get current level and problem description

120 L = self.level

121 P = L.prob

122 # set initial guess for gradient to zero

123 L.f[0] = P.dtype_f(init=P.init, val=0.0)

124 for m in range(1, self.coll.num_nodes + 1):

125 # copy u[0] to all collocation nodes and set f (the gradient) to zero

127 L.u[m] = P.dtype_u(L.u[0])

128 L.f[m] = P.dtype_f(init=P.init, val=0.0)

129 elif self.params.initial_guess == 'zero':

130 L.u[m] = P.dtype_u(init=P.init, val=0.0)

131 L.f[m] = P.dtype_f(init=P.init, val=0.0)

133 elif self.params.initial_guess == 'random':

134 L.u[m] = P.dtype_u(init=P.init, val=np.random.rand(1)[0])

135 L.f[m] = P.dtype_f(init=P.init, val=np.random.rand(1)[0])

136 else:

137 raise ParameterError(f'initial_guess option {self.params.initial_guess} not implemented')

139 # indicate that this level is now ready for sweeps

140 L.status.unlocked = True

141 L.status.updated = True

143 def compute_residual(self, stage=None):

144 r"""

145 Uses the absolute value of the DAE system

147 .. math::

148 ||F(t, u, u')||

150 for computing the residual in a chosen norm.

152 Parameters

153 ----------

154 stage : str, optional

155 The current stage of the step the level belongs to.

156 """

158 # get current level and problem description

159 L = self.level

160 P = L.prob

162 # Check if we want to skip the residual computation to gain performance

163 # Keep in mind that skipping any residual computation is likely to give incorrect outputs of the residual!

164 if stage in self.params.skip_residual_computation:

165 L.status.residual = 0.0 if L.status.residual is None else L.status.residual

166 return None

168 # compute the residual for each node

169 res_norm = []

170 for m in range(self.coll.num_nodes):

171 # use abs function from data type here

172 res_norm.append(abs(P.eval_f(L.u[m + 1], L.f[m + 1], L.time + L.dt * self.coll.nodes[m])))

174 # find maximal residual over the nodes

175 if L.params.residual_type == 'full_abs':

176 L.status.residual = max(res_norm)

177 elif L.params.residual_type == 'last_abs':

178 L.status.residual = res_norm[-1]

179 elif L.params.residual_type == 'full_rel':

180 L.status.residual = max(res_norm) / abs(L.u[0])

181 elif L.params.residual_type == 'last_rel':

182 L.status.residual = res_norm[-1] / abs(L.u[0])

183 else:

184 raise ParameterError(

185 f'residual_type = {L.params.residual_type} not implemented, choose '

186 f'full_abs, last_abs, full_rel or last_rel instead'

187 )

189 # indicate that the residual has seen the new values

190 L.status.updated = False

192 return None

194 def compute_end_point(self):

195 """

196 Compute u at the right point of the interval.

198 The value uend computed here is a full evaluation of the Picard formulation unless do_full_update==False

200 Returns:

201 None

202 """

204 if not self.coll.right_is_node or self.params.do_coll_update:

205 raise NotImplementedError()

207 super().compute_end_point()

209 @staticmethod

210 def F(du, P, factor, u_approx, t):

211 r"""

212 This function builds the implicit system to be solved for a DAE of the form

214 .. math::

215 0 = F(u, u', t)

217 Applying SDC yields the (non)-linear system to be solved

219 .. math::

220 0 = F(u_0 + \sum_{j=1}^M (q_{mj}-\tilde{q}_{mj}) U_j

221 + \sum_{j=1}^m \tilde{q}_{mj} U_j, U_m, \tau_m),

223 which is solved for the derivative of u.

225 Note

226 ----

227 This function is also used for Runge-Kutta methods since only

228 the argument u_approx differs. However, the implicit system to be solved

229 is the same.

231 Parameters

232 ----------

233 du : dtype_u

234 Unknowns of the system (derivative of solution u).

235 P : pySDC.projects.DAE.misc.problemDAE.ProblemDAE

236 Problem class.

237 factor : float

238 Abbrev. for the node-to-node stepsize.

239 u_approx : dtype_u

240 Approximation to numerical solution :math:u.

241 t : float

242 Current time :math:t.

244 Returns

245 -------

246 sys : dtype_f

247 System to be solved.

248 """

250 local_du_approx = P.dtype_f(du)

252 # build parameters to pass to implicit function

253 local_u_approx = P.dtype_f(u_approx)

255 # note that derivatives of algebraic variables are taken into account here too

256 # these do not directly affect the output of eval_f but rather indirectly via QI

257 local_u_approx += factor * local_du_approx

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

260 return sys