Coverage for pySDC/projects/Monodomain/sweeper_classes/runge_kutta/imexexp_1st_order.py: 96%

1import numpy as np 

2 

3from pySDC.core.Sweeper import sweeper 

4from pySDC.core.Errors import CollocationError 

5 

6 

7class imexexp_1st_order(sweeper): 

8 """ 

9 Custom sweeper class, implements Sweeper.py 

10 

11 First-order IMEXEXP sweeper using implicit/explicit/exponential Euler as base integrator 

12 In the cardiac electrphysiology community this is known as Rush-Larsen scheme. 

13 

14 """ 

15 

16 def __init__(self, params): 

17 """ 

18 Initialization routine for the custom sweeper 

19 

20 Args: 

21 params: parameters for the sweeper 

22 """ 

23 

24 if "QI" not in params: 

25 params["QI"] = "IE" 

26 

27 # call parent's initialization routine 

28 super(imexexp_1st_order, self).__init__(params) 

29 

30 # IMEX integration matrices 

31 self.QI = self.get_Qdelta_implicit(coll=self.coll, qd_type=self.params.QI) 

32 self.delta = np.diagonal(self.QI)[1:] 

33 

34 def eval_phi_f_exp(self, u, factor): 

35 """ 

36 Evaluates the exponential part of the right-hand side f_exp(u)=lambda(u)*(u-y_inf(u)) multiplied by the exponential factor phi_1(factor*lambda) 

37 Since phi_1(z)=(e^z-1)/z then phi_1(factor*lambda) * f_exp(u) = ((e^(factor*lambda)-1)/factor) *(u-y_inf(u)) 

38 """ 

39 L = self.level 

40 P = L.prob 

41 self.lmbda = P.dtype_u(init=P.init, val=0.0) 

42 self.yinf = P.dtype_u(init=P.init, val=0.0) 

43 P.eval_lmbda_yinf_exp(u, self.lmbda, self.yinf) 

44 phi_f_exp = P.dtype_u(init=P.init, val=0.0) 

45 for i in P.rhs_exp_indeces: 

46 phi_f_exp[i][:] = u[i] - self.yinf[i][:] 

47 phi_f_exp[i][:] *= (np.exp(factor * self.lmbda[i]) - 1.0) / factor 

48 

49 return phi_f_exp 

50 

51 def integrate(self): 

52 """ 

53 Integrates the right-hand side (here impl + expl + exp) 

54 

55 Returns: 

56 list of dtype_u: containing the integral as values 

57 """ 

58 

59 # get current level and problem description 

60 L = self.level 

61 

62 me = [] 

63 

64 # integrate RHS over all collocation nodes 

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

66 me.append(L.dt * self.coll.Qmat[m, 1] * (L.f[1].impl + L.f[1].expl + L.f[1].exp)) 

67 for j in range(2, self.coll.num_nodes + 1): 

68 me[m - 1] += L.dt * self.coll.Qmat[m, j] * (L.f[j].impl + L.f[j].expl + L.f[j].exp) 

69 

70 return me 

71 

72 def update_nodes(self): 

73 """ 

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

75 

76 Returns: 

77 None 

78 """ 

79 

80 # get current level and problem description 

81 L = self.level 

82 P = L.prob 

83 

84 # only if the level has been touched before 

85 assert L.status.unlocked 

86 

87 # get number of collocation nodes for easier access 

88 M = self.coll.num_nodes 

89 

90 integral = self.integrate() 

91 for m in range(M): 

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

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

94 for i in range(1, M): 

95 integral[M - i] -= integral[M - i - 1] 

96 

97 # do the sweep 

98 for m in range(M): 

99 integral[m] -= ( 

100 L.dt 

101 * self.delta[m] 

102 * (L.f[m].expl + L.f[m + 1].impl + self.eval_phi_f_exp(L.u[m], L.dt * self.delta[m])) 

103 ) 

104 for m in range(M): 

105 rhs = ( 

106 L.u[m] 

107 + integral[m] 

108 + L.dt * self.delta[m] * (L.f[m].expl + self.eval_phi_f_exp(L.u[m], L.dt * self.delta[m])) 

109 ) 

110 

111 # implicit solve with prefactor stemming from QI 

112 L.u[m + 1] = P.solve_system(rhs, L.dt * self.delta[m], L.u[m + 1], L.time + L.dt * self.coll.nodes[m]) 

113 

114 # update function values 

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

116 

117 # indicate presence of new values at this level 

118 L.status.updated = True 

119 

120 return None 

121 

122 def compute_end_point(self): 

123 """ 

124 Compute u at the right point of the interval 

125 

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

127 

128 Returns: 

129 None 

130 """ 

131 

132 # get current level and problem description 

133 L = self.level 

134 P = L.prob 

135 

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

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

138 # a copy is sufficient 

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

140 else: 

141 raise CollocationError( 

142 "In this sweeper we expect the right point to be a collocation node and do_coll_update==False" 

143 ) 

144 

145 return None