Coverage for pySDC/implementations/sweeper_classes/generic_implicit_mass.py: 93%
61 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1from pySDC.core.errors import ParameterError
2from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit
5class generic_implicit_mass(generic_implicit):
6 """Fully implicit SDC sweeper for the mass-matrix formulation (no M^-1 anywhere)."""
8 def update_nodes(self):
9 """
10 Fully implicit sweep for problems posed in mass-matrix form, M u' = f(u), where eval_f
11 returns the load vector f and solve_system solves (M - factor J) u = rhs with rhs already
12 in the dual space. u0 carries the mass matrix on the finest level only; coarse levels
13 receive it already applied from the space-time transfer.
15 Returns:
16 None
17 """
19 L = self.level
20 P = L.prob
22 # only if the level has been touched before
23 assert L.status.unlocked
25 # get number of collocation nodes for easier access
26 M = self.coll.num_nodes
28 # gather all terms which are known already (e.g. from the previous iteration)
29 # this corresponds to u0 + QF(u^k) - QdF(u^k) + tau
31 # get QF(u^k)
32 integral = self.integrate()
34 # This is somewhat ugly, but we have to apply the mass matrix on u0 only on the finest level
35 if L.level_index == 0:
36 u0 = P.apply_mass_matrix(L.u[0])
37 else:
38 u0 = L.u[0]
40 for m in range(M):
41 # get -QdF(u^k)_m
42 for j in range(1, M + 1):
43 integral[m] -= L.dt * self.QI[m + 1, j] * L.f[j]
45 # add initial value
46 integral[m] += u0
47 # add tau if associated
48 if L.tau[m] is not None:
49 integral[m] += L.tau[m]
51 # do the sweep
52 for m in range(0, M):
53 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1)
54 rhs = P.dtype_u(integral[m])
55 for j in range(1, m + 1):
56 rhs += L.dt * self.QI[m + 1, j] * L.f[j]
58 # implicit solve with prefactor stemming from the diagonal of Qd
59 L.u[m + 1] = P.solve_system(
60 rhs, L.dt * self.QI[m + 1, m + 1], L.u[m + 1], L.time + L.dt * self.coll.nodes[m]
61 )
63 # update function values
64 L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m])
66 # indicate presence of new values at this level
67 L.status.updated = True
69 return None
71 def compute_end_point(self):
72 """
73 Compute u at the right point of the interval
75 The value uend computed here is a full evaluation of the Picard formulation unless do_full_update==False
77 Returns:
78 None
79 """
81 # get current level and problem description
82 L = self.level
83 P = L.prob
85 # check if Mth node is equal to right point and do_coll_update is false, perform a simple copy
86 if self.coll.right_is_node and not self.params.do_coll_update:
87 # a copy is sufficient
88 L.uend = P.dtype_u(L.u[-1])
89 # On coarse levels u[0] is carried in the dual space: the transfer hands down P^T M u0 and
90 # the sweep consumes it as-is. PFASST copies uend straight into the next step's u[0]
91 # (controller.recv_full), so what we hand over has to be dual as well -- otherwise the step
92 # boundary feeds a primal value into a slot that is read as M u0.
93 if L.level_index > 0:
94 L.uend = P.apply_mass_matrix(L.uend)
95 else:
96 raise NotImplementedError('Mass matrix sweeper expect u_M = u_end')
98 return None
100 def compute_residual(self, stage=None):
101 """
102 Computation of the residual using the collocation matrix Q
104 Args:
105 stage (str): The current stage of the step the level belongs to
106 """
108 # get current level and problem description
109 L = self.level
110 P = L.prob
112 # Check if we want to skip the residual computation to gain performance
113 # Keep in mind that skipping any residual computation is likely to give incorrect outputs of the residual!
114 if stage in self.params.skip_residual_computation:
115 L.status.residual = 0.0 if L.status.residual is None else L.status.residual
116 return None
118 # check if there are new values (e.g. from a sweep)
119 # assert L.status.updated
121 # compute the residual for each node
123 # build QF(u)
124 res_norm = []
125 res = self.integrate()
126 for m in range(self.coll.num_nodes):
128 # This is somewhat ugly, but we have to apply the mass matrix on u0 only on the finest level
130 if L.level_index == 0:
131 res[m] += P.apply_mass_matrix(L.u[0] - L.u[m + 1])
132 else:
133 res[m] += L.u[0] - P.apply_mass_matrix(L.u[m + 1])
134 # add tau if associated
135 if L.tau[m] is not None:
136 res[m] += L.tau[m]
138 # Due to different boundary conditions we might have to fix the residual
139 if getattr(L.prob, 'fix_bc_for_residual', False):
140 L.prob.fix_residual(res[m])
141 # use abs function from data type here
142 res_norm.append(abs(res[m]))
144 # find maximal residual over the nodes
145 if L.params.residual_type == 'full_abs':
146 L.status.residual = max(res_norm)
147 elif L.params.residual_type == 'last_abs':
148 L.status.residual = res_norm[-1]
149 elif L.params.residual_type == 'full_rel':
150 L.status.residual = max(res_norm) / abs(L.u[0])
151 elif L.params.residual_type == 'last_rel':
152 L.status.residual = res_norm[-1] / abs(L.u[0])
153 else:
154 raise ParameterError(
155 f'residual_type = {L.params.residual_type} not implemented, choose '
156 f'full_abs, last_abs, full_rel or last_rel instead'
157 )
159 # indicate that the residual has seen the new values
160 L.status.updated = False
162 return None