Coverage for pySDC/implementations/sweeper_classes/imex_1st_order_mass.py: 100%

1from pySDC.implementations.sweeper_classes.generic_implicit_mass import generic_implicit_mass 

2from pySDC.implementations.sweeper_classes.imex_1st_order import imex_1st_order 

3 

4 

5class imex_1st_order_mass(imex_1st_order): 

6 """ 

7 IMEX-SDC sweeper for problems M u' = f(u) with a mass matrix M, as they arise from finite elements. 

8 

9 First-order IMEX sweeper using implicit/explicit Euler as base integrator, with mass or weighting matrix 

10 

11 The end point is only available as a copy of the last node: ``compute_end_point`` raises unless the right end of 

12 the interval is a node and ``do_coll_update`` is off. 

13 """ 

14 

15 def update_nodes(self): 

16 """ 

17 Update the u- and f-values at the collocation nodes -> corresponds to a single sweep over all nodes 

18 

19 Returns: 

20 None 

21 """ 

22 

23 # get current level and problem description 

24 L = self.level 

25 P = L.prob 

26 

27 # only if the level has been touched before 

28 assert L.status.unlocked 

29 

30 # get number of collocation nodes for easier access 

31 M = self.coll.num_nodes 

32 

33 # gather all terms which are known already (e.g. from the previous iteration) 

34 # this corresponds to u0 + QF(u^k) - QIFI(u^k) - QEFE(u^k) + tau 

35 

36 # get QF(u^k) 

37 integral = self.integrate() 

38 

39 # This is somewhat ugly, but we have to apply the mass matrix on u0 only on the finest level 

40 if L.level_index == 0: 

41 u0 = P.apply_mass_matrix(L.u[0]) 

42 else: 

43 u0 = L.u[0] 

44 

45 for m in range(M): 

46 # subtract QIFI(u^k)_m + QEFE(u^k)_m 

47 for j in range(M + 1): 

48 integral[m] -= L.dt * (self.QI[m + 1, j] * L.f[j].impl + self.QE[m + 1, j] * L.f[j].expl) 

49 # add initial value 

50 integral[m] += u0 

51 # add tau if associated 

52 if L.tau[m] is not None: 

53 integral[m] += L.tau[m] 

54 

55 # do the sweep 

56 for m in range(0, M): 

57 # build rhs, consisting of the known values from above and new values from previous nodes (at k+1) 

58 rhs = P.dtype_u(integral[m]) 

59 for j in range(m + 1): 

60 rhs += L.dt * (self.QI[m + 1, j] * L.f[j].impl + self.QE[m + 1, j] * L.f[j].expl) 

61 

62 # implicit solve with prefactor stemming from QI 

63 L.u[m + 1] = P.solve_system( 

64 rhs, L.dt * self.QI[m + 1, m + 1], L.u[m + 1], L.time + L.dt * self.coll.nodes[m] 

65 ) 

66 # update function values 

67 L.f[m + 1] = P.eval_f(L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) 

68 

69 # indicate presence of new values at this level 

70 L.status.updated = True 

71 

72 return None 

73 

74 def compute_end_point(self): 

75 """ 

76 Compute u at the right point of the interval 

77 

78 This is a copy of the last node. The mass matrix sweeper does not support the collocation update, so this 

79 raises unless the right end of the interval is a node and do_coll_update is False. 

80 

81 Returns: 

82 None 

83 """ 

84 

85 # get current level and problem description 

86 L = self.level 

87 P = L.prob 

88 

89 # check if Mth node is equal to right point and do_coll_update is false, perform a simple copy 

90 if self.coll.right_is_node and not self.params.do_coll_update: 

91 # a copy is sufficient 

92 L.uend = P.dtype_u(L.u[-1]) 

93 # On coarse levels u[0] is carried in the dual space: the transfer hands down P^T M u0 and 

94 # the sweep consumes it as-is. PFASST copies uend straight into the next step's u[0] 

95 # (controller.recv_full), so what we hand over has to be dual as well -- otherwise the step 

96 # boundary feeds a primal value into a slot that is read as M u0. 

97 if L.level_index > 0: 

98 L.uend = P.apply_mass_matrix(L.uend) 

99 else: 

100 raise NotImplementedError('Mass matrix sweeper expect u_M = u_end') 

101 

102 return None 

103 

104 # same residual as the implicit mass sweeper: the mass matrix enters it the same way 

105 compute_residual = generic_implicit_mass.compute_residual