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 Scalar test equation with a fast and a slow wave, IMEX with the fast part implicit and the slow part explicit. 

12 

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

14 

15 .. math:: 

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

17 

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

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

20 is given by 

21 

22 .. math:: 

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

24 

25 Parameters 

26 ---------- 

27 lambda_s : np.1darray, optional 

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

29 lambda_f : np.1darray, optional 

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

31 u0 : np.1darray, optional 

32 Initial condition of the problem. 

33 

34 References 

35 ---------- 

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

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

38 """ 

39 

40 dtype_u = mesh 

41 dtype_f = imex_mesh 

42 

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

44 """Initialization routine""" 

45 

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

47 super().__init__(init) 

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

49 

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

51 r""" 

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

53 

54 Parameters 

55 ---------- 

56 rhs : dtype_f 

57 Right-hand side for the nonlinear system. 

58 factor : float 

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

60 u0 : dtype_u 

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

62 t : float 

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

64 

65 Returns 

66 ------- 

67 me : dtype_u 

68 The solution as mesh. 

69 """ 

70 

71 me = self.dtype_u(self.init) 

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

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

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

75 

76 return me 

77 

78 def __eval_fexpl(self, u, t): 

79 """ 

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

81 

82 Parameters 

83 ---------- 

84 u : dtype_u 

85 Current values of the numerical solution. 

86 t : float 

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

88 

89 Returns 

90 ------- 

91 fexpl : dtype_u 

92 Explicit part of right-hand side. 

93 """ 

94 

95 fexpl = self.dtype_u(self.init) 

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

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

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

99 return fexpl 

100 

101 def __eval_fimpl(self, u, t): 

102 """ 

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

104 

105 Parameters 

106 ---------- 

107 u : dtype_u 

108 Current values of the numerical solution. 

109 t : float 

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

111 

112 Returns 

113 ------- 

114 fimpl : dtype_u 

115 Implicit part of right-hand side. 

116 """ 

117 

118 fimpl = self.dtype_u(self.init) 

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

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

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

122 

123 return fimpl 

124 

125 def eval_f(self, u, t): 

126 """ 

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

128 

129 Parameters 

130 ---------- 

131 u : dtype_u 

132 Current values of the numerical solution. 

133 t : float 

134 Current time at which the numerical solution is computed. 

135 

136 Returns 

137 ------- 

138 f : dtype_f 

139 The right-hand side divided into two parts. 

140 """ 

141 

142 f = self.dtype_f(self.init) 

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

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

145 return f 

146 

147 def u_exact(self, t): 

148 r""" 

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

150 

151 Parameters 

152 ---------- 

153 t : float 

154 Time of the exact solution. 

155 

156 Returns 

157 ------- 

158 me : dtype_u 

159 The exact solution. 

160 """ 

161 

162 me = self.dtype_u(self.init) 

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

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

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

166 return me