Coverage for pySDC/implementations/sweeper_classes/generic_implicit.py: 98%
50 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-09-04 15:08 +0000
« prev ^ index » next coverage.py v7.10.6, created at 2025-09-04 15:08 +0000
1from pySDC.core.sweeper import Sweeper
4class generic_implicit(Sweeper):
5 """
6 Generic implicit sweeper, expecting lower triangular matrix type as input
8 Attributes:
9 QI: lower triangular matrix
10 """
12 def __init__(self, params, level):
13 """
14 Initialization routine for the custom sweeper
16 Args:
17 params: parameters for the sweeper
18 level (pySDC.Level.level): the level that uses this sweeper
19 """
21 if 'QI' not in params:
22 params['QI'] = 'IE'
24 super().__init__(params, level)
26 # get QI matrix
27 self.QI = self.get_Qdelta_implicit(qd_type=self.params.QI)
29 def integrate(self):
30 """
31 Integrates the right-hand side
33 Returns:
34 list of dtype_u: containing the integral as values
35 """
37 L = self.level
38 P = L.prob
40 me = []
42 # integrate RHS over all collocation nodes
43 for m in range(1, self.coll.num_nodes + 1):
44 # new instance of dtype_u, initialize values with 0
45 me.append(P.dtype_u(P.init, val=0.0))
46 for j in range(1, self.coll.num_nodes + 1):
47 me[-1] += L.dt * self.coll.Qmat[m, j] * L.f[j]
49 return me
51 def update_nodes(self):
52 """
53 Update the u- and f-values at the collocation nodes -> corresponds to a single sweep over all nodes
55 Returns:
56 None
57 """
59 L = self.level
60 P = L.prob
62 # only if the level has been touched before
63 assert L.status.unlocked
65 # get number of collocation nodes for easier access
66 M = self.coll.num_nodes
68 # gather all terms which are known already (e.g. from the previous iteration)
69 # this corresponds to u0 + QF(u^k) - QdF(u^k) + tau
71 # get QF(u^k)
72 integral = self.integrate()
73 for m in range(M):
74 # get -QdF(u^k)_m
75 for j in range(1, M + 1):
76 integral[m] -= L.dt * self.QI[m + 1, j] * L.f[j]
78 # add initial value
79 integral[m] += L.u[0]
80 # add tau if associated
81 if L.tau[m] is not None:
82 integral[m] += L.tau[m]
84 # do the sweep
85 for m in range(0, M):
86 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1)
87 rhs = P.dtype_u(integral[m])
88 for j in range(1, m + 1):
89 rhs += L.dt * self.QI[m + 1, j] * L.f[j]
91 # implicit solve with prefactor stemming from the diagonal of Qd
92 alpha = L.dt * self.QI[m + 1, m + 1]
93 if alpha == 0:
94 L.u[m + 1] = rhs
95 else:
96 L.u[m + 1] = P.solve_system(rhs, alpha, L.u[m + 1], L.time + L.dt * self.coll.nodes[m])
97 # update function values
98 L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m])
100 # indicate presence of new values at this level
101 L.status.updated = True
103 return None
105 def compute_end_point(self):
106 """
107 Compute u at the right point of the interval
109 The value uend computed here is a full evaluation of the Picard formulation unless do_full_update==False
111 Returns:
112 None
113 """
115 L = self.level
116 P = L.prob
118 # check if Mth node is equal to right point and do_coll_update is false, perform a simple copy
119 if self.coll.right_is_node and not self.params.do_coll_update:
120 # a copy is sufficient
121 L.uend = P.dtype_u(L.u[-1])
122 else:
123 # start with u0 and add integral over the full interval (using coll.weights)
124 L.uend = P.dtype_u(L.u[0])
125 for m in range(self.coll.num_nodes):
126 L.uend += L.dt * self.coll.weights[m] * L.f[m + 1]
127 # add up tau correction of the full interval (last entry)
128 if L.tau[-1] is not None:
129 L.uend += L.tau[-1]
131 return None