Coverage for pySDC/implementations/problem_classes/FastWaveSlowWave_0D.py: 100%

1import numpy as np 

2 

3from pySDC.core.errors import ParameterError 

4from pySDC.core.problem import Problem 

5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh 

6 

7 

8# noinspection PyUnusedLocal 

9class swfw_scalar(Problem): 

10 r""" 

11 This class implements the fast-wave-slow-wave scalar problem fully investigated in [1]_. It is defined by 

12 

13 .. math:: 

14 \frac{d u(t)}{dt} = \lambda_f u(t) + \lambda_s u(t), 

15 

16 where :math:`\lambda_f` denotes the part of the fast wave, and :math:`\lambda_s` is the part of the slow wave with 

17 :math:`\lambda_f \gg \lambda_s`. Let :math:`u_0` be the initial condition to the problem, then the exact solution 

18 is given by 

19 

20 .. math:: 

21 u(t) = u_0 \exp((\lambda_f + \lambda_s) t). 

22 

23 Parameters 

24 ---------- 

25 lambda_s : np.1darray, optional 

26 Part of the slow wave :math:`\lambda_s`. 

27 lambda_f : np.1darray, optional 

28 Part of the fast wave :math:`\lambda_f`. 

29 u0 : np.1darray, optional 

30 Initial condition of the problem. 

31 

32 References 

33 ---------- 

34 .. [1] D. Ruprecht, R. Speck. Spectral deferred corrections with fast-wave slow-wave splitting. 

35 SIAM J. Sci. Comput. Vol. 38 No. 4 (2016). 

36 """ 

37 

38 dtype_u = mesh 

39 dtype_f = imex_mesh 

40 

41 def __init__(self, lambda_s=-1, lambda_f=-1000, u0=1): 

42 """Initialization routine""" 

43 

44 init = ([lambda_s.size, lambda_f.size], None, np.dtype('complex128')) 

45 super().__init__(init) 

46 self._makeAttributeAndRegister('lambda_s', 'lambda_f', 'u0', localVars=locals(), readOnly=True) 

47 

48 def solve_system(self, rhs, factor, u0, t): 

49 r""" 

50 Simple im=nversion of :math:`(1 - \Delta t \cdot \lambda)\vec{u} = \vec{rhs}`. 

51 

52 Parameters 

53 ---------- 

54 rhs : dtype_f 

55 Right-hand side for the nonlinear system. 

56 factor : float 

57 Abbrev. for the node-to-node stepsize (or any other factor required). 

58 u0 : dtype_u 

59 Initial guess for the iterative solver (not used here so far). 

60 t : float 

61 Current time (e.g. for time-dependent BCs). 

62 

63 Returns 

64 ------- 

65 me : dtype_u 

66 The solution as mesh. 

67 """ 

68 

69 me = self.dtype_u(self.init) 

70 for i in range(self.lambda_s.size): 

71 for j in range(self.lambda_f.size): 

72 me[i, j] = rhs[i, j] / (1.0 - factor * self.lambda_f[j]) 

73 

74 return me 

75 

76 def __eval_fexpl(self, u, t): 

77 """ 

78 Helper routine to evaluate the explicit part of the right-hand side. 

79 

80 Parameters 

81 ---------- 

82 u : dtype_u 

83 Current values of the numerical solution. 

84 t : float 

85 Current time at which the numerical solution is computed (not used here). 

86 

87 Returns 

88 ------- 

89 fexpl : dtype_u 

90 Explicit part of right-hand side. 

91 """ 

92 

93 fexpl = self.dtype_u(self.init) 

94 for i in range(self.lambda_s.size): 

95 for j in range(self.lambda_f.size): 

96 fexpl[i, j] = self.lambda_s[i] * u[i, j] 

97 return fexpl 

98 

99 def __eval_fimpl(self, u, t): 

100 """ 

101 Helper routine to evaluate the implicit part of the right-hand side. 

102 

103 Parameters 

104 ---------- 

105 u : dtype_u 

106 Current values of the numerical solution. 

107 t : float 

108 Current time at which the numerical solution is computed (not used here). 

109 

110 Returns 

111 ------- 

112 fimpl : dtype_u 

113 Implicit part of right-hand side. 

114 """ 

115 

116 fimpl = self.dtype_u(self.init) 

117 for i in range(self.lambda_s.size): 

118 for j in range(self.lambda_f.size): 

119 fimpl[i, j] = self.lambda_f[j] * u[i, j] 

120 

121 return fimpl 

122 

123 def eval_f(self, u, t): 

124 """ 

125 Routine to evaluate both parts of the right-hand side of the problem. 

126 

127 Parameters 

128 ---------- 

129 u : dtype_u 

130 Current values of the numerical solution. 

131 t : float 

132 Current time at which the numerical solution is computed. 

133 

134 Returns 

135 ------- 

136 f : dtype_f 

137 The right-hand side divided into two parts. 

138 """ 

139 

140 f = self.dtype_f(self.init) 

141 f.impl = self.__eval_fimpl(u, t) 

142 f.expl = self.__eval_fexpl(u, t) 

143 return f 

144 

145 def u_exact(self, t): 

146 r""" 

147 Routine to compute the exact solution at time :math:`t`. 

148 

149 Parameters 

150 ---------- 

151 t : float 

152 Time of the exact solution. 

153 

154 Returns 

155 ------- 

156 me : dtype_u 

157 The exact solution. 

158 """ 

159 

160 me = self.dtype_u(self.init) 

161 for i in range(self.lambda_s.size): 

162 for j in range(self.lambda_f.size): 

163 me[i, j] = self.u0 * np.exp((self.lambda_f[j] + self.lambda_s[i]) * t) 

164 return me