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

1import numpy as np 

2from scipy.sparse.linalg import spsolve 

3 

4from pySDC.core.errors import ParameterError 

5from pySDC.core.problem import Problem, WorkCounter 

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

7from pySDC.implementations.problem_classes.acoustic_helpers.buildWave1DMatrix import ( 

8 getWave1DMatrix, 

9 getWave1DAdvectionMatrix, 

10) 

11 

12 

13# noinspection PyUnusedLocal 

14class acoustic_1d_imex(Problem): 

15 r""" 

16 Periodic 1D acoustic-advection system with finite differences, acoustic waves implicit and advection explicit. 

17 

18 This class implements the one-dimensional acoustics advection equation on a periodic domain :math:`[0, 1]` 

19 fully investigated in [1]_. The equations are given by 

20 

21 .. math:: 

22 \frac{\partial u}{\partial t} + c_s \frac{\partial p}{\partial x} + U \frac{\partial u}{\partial x} = 0, 

23 

24 .. math:: 

25 \frac{\partial p}{\partial t} + c_s \frac{\partial u}{\partial x} + U \frac{\partial p}{\partial x} = 0. 

26 

27 For initial data :math:`u(x, 0) \equiv 0` and :math:`p(x, 0) = p_0 (x)` the analytical solution is 

28 

29 .. math:: 

30 u(x, t) = \frac{1}{2} p_0 (x - (U + c_s) t) - \frac{1}{2} p_0 (x - (U - c_s) t), 

31 

32 .. math:: 

33 p(x, t) = \frac{1}{2} p_0 (x - (U + c_s) t) + \frac{1}{2} p_0 (x - (U - c_s) t). 

34 

35 The problem is implemented in the way that is used for IMEX time-stepping. 

36 

37 Parameters 

38 ---------- 

39 nvars : int, optional 

40 Number of degrees of freedom. 

41 cs : float, optional 

42 Sound velocity :math:`c_s`. 

43 cadv : float, optional 

44 Advection speed :math:`U`. 

45 order_adv : int, optional 

46 Order of which the advective derivative is discretized. 

47 waveno : int, optional 

48 The wave number. 

49 

50 Attributes 

51 ---------- 

52 mesh : np.1darray 

53 1d mesh. 

54 dx : float 

55 Mesh size. 

56 Dx : scipy.csc_matrix 

57 Matrix for the advection operator. 

58 Id : scipy.csc_matrix 

59 Sparse identity matrix. 

60 A : scipy.csc_matrix 

61 Matrix for the wave operator. 

62 

63 References 

64 ---------- 

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

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

67 """ 

68 

69 dtype_u = mesh 

70 dtype_f = imex_mesh 

71 

72 def __init__(self, nvars=None, cs=0.5, cadv=0.1, order_adv=5, waveno=5): 

73 """Initialization routine""" 

74 

75 if nvars is None: 

76 nvars = (2, 300) 

77 

78 # invoke super init, passing number of dofs 

79 super().__init__((nvars, None, np.dtype('float64'))) 

80 self._makeAttributeAndRegister('nvars', 'cs', 'cadv', 'order_adv', 'waveno', localVars=locals(), readOnly=True) 

81 

82 self.mesh = np.linspace(0.0, 1.0, self.nvars[1], endpoint=False) 

83 self.dx = self.mesh[1] - self.mesh[0] 

84 

85 self.Dx = -self.cadv * getWave1DAdvectionMatrix(self.nvars[1], self.dx, self.order_adv) 

86 self.Id, A = getWave1DMatrix(self.nvars[1], self.dx, ['periodic', 'periodic'], ['periodic', 'periodic']) 

87 self.A = -self.cs * A 

88 

89 self.work_counters['rhs'] = WorkCounter() 

90 

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

92 r""" 

93 Simple linear solver for :math:`(I-factor\cdot A)\vec{u}=\vec{rhs}`. 

94 

95 Parameters 

96 ---------- 

97 rhs : dtype_f 

98 Right-hand side for the linear system. 

99 factor : float 

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

101 u0 : dtype_u 

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

103 t : float 

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

105 

106 Returns 

107 ------- 

108 me : dtype_u 

109 The solution as mesh. 

110 """ 

111 

112 M = self.Id - factor * self.A 

113 

114 b = np.concatenate((rhs[0, :], rhs[1, :])) 

115 

116 sol = spsolve(M, b) 

117 

118 me = self.dtype_u(self.init) 

119 me[0, :], me[1, :] = np.split(sol, 2) 

120 

121 return me 

122 

123 def __eval_fexpl(self, u, t): 

124 """ 

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

126 

127 Parameters 

128 ---------- 

129 u : dtype_u 

130 Current values of the numerical solution. 

131 t : float 

132 Current time of the numerical solution is computed (not used here). 

133 

134 Returns 

135 ------- 

136 fexpl : dtype_f 

137 Explicit part of the right-hand side. 

138 """ 

139 

140 b = np.concatenate((u[0, :], u[1, :])) 

141 sol = self.Dx.dot(b) 

142 

143 fexpl = self.dtype_u(self.init) 

144 fexpl[0, :], fexpl[1, :] = np.split(sol, 2) 

145 

146 return fexpl 

147 

148 def __eval_fimpl(self, u, t): 

149 """ 

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

151 

152 Parameters 

153 ---------- 

154 u : dtype_u 

155 Current values of the numerical solution. 

156 t : float 

157 Current time of the numerical solution is computed (not used here). 

158 

159 Returns 

160 ------- 

161 fimpl : dtype_f 

162 Implicit part of the right-hand side. 

163 """ 

164 

165 b = np.concatenate((u[:][0, :], u[:][1, :])) 

166 sol = self.A.dot(b) 

167 

168 fimpl = self.dtype_u(self.init, val=0.0) 

169 fimpl[0, :], fimpl[1, :] = np.split(sol, 2) 

170 

171 return fimpl 

172 

173 def eval_f(self, u, t): 

174 """ 

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

176 

177 Parameters 

178 ---------- 

179 u : dtype_u 

180 Current values of the numerical solution. 

181 t : float 

182 Current time of the numerical solution is computed. 

183 

184 Returns 

185 ------- 

186 f : dtype_f 

187 The right-hand side divided into two parts. 

188 """ 

189 

190 f = self.dtype_f(self.init) 

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

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

193 

194 self.work_counters['rhs']() 

195 return f 

196 

197 def u_exact(self, t): 

198 r""" 

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

200 

201 Parameters 

202 ---------- 

203 t : float 

204 Time of the exact solution. 

205 

206 Returns 

207 ------- 

208 me : dtype_u 

209 The exact solution. 

210 """ 

211 

212 def u_initial(x, k): 

213 return np.sin(k * 2.0 * np.pi * x) + np.sin(2.0 * np.pi * x) 

214 

215 me = self.dtype_u(self.init) 

216 me[0, :] = 0.5 * u_initial(self.mesh - (self.cadv + self.cs) * t, self.waveno) - 0.5 * u_initial( 

217 self.mesh - (self.cadv - self.cs) * t, self.waveno 

218 ) 

219 me[1, :] = 0.5 * u_initial(self.mesh - (self.cadv + self.cs) * t, self.waveno) + 0.5 * u_initial( 

220 self.mesh - (self.cadv - self.cs) * t, self.waveno 

221 ) 

222 return me