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 ptype, 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(ptype): 

15 r""" 

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

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

18 

19 .. math:: 

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

21 

22 .. math:: 

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

24 

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

26 

27 .. math:: 

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

29 

30 .. math:: 

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

32 

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

34 

35 Parameters 

36 ---------- 

37 nvars : int, optional 

38 Number of degrees of freedom. 

39 cs : float, optional 

40 Sound velocity :math:`c_s`. 

41 cadv : float, optional 

42 Advection speed :math:`U`. 

43 order_adv : int, optional 

44 Order of which the advective derivative is discretized. 

45 waveno : int, optional 

46 The wave number. 

47 

48 Attributes 

49 ---------- 

50 mesh : np.1darray 

51 1d mesh. 

52 dx : float 

53 Mesh size. 

54 Dx : scipy.csc_matrix 

55 Matrix for the advection operator. 

56 Id : scipy.csc_matrix 

57 Sparse identity matrix. 

58 A : scipy.csc_matrix 

59 Matrix for the wave operator. 

60 

61 References 

62 ---------- 

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

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

65 """ 

66 

67 dtype_u = mesh 

68 dtype_f = imex_mesh 

69 

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

71 """Initialization routine""" 

72 

73 if nvars is None: 

74 nvars = (2, 300) 

75 

76 # invoke super init, passing number of dofs 

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

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

79 

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

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

82 

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

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

85 self.A = -self.cs * A 

86 

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

88 

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

90 r""" 

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

92 

93 Parameters 

94 ---------- 

95 rhs : dtype_f 

96 Right-hand side for the linear system. 

97 factor : float 

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

99 u0 : dtype_u 

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

101 t : float 

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

103 

104 Returns 

105 ------- 

106 me : dtype_u 

107 The solution as mesh. 

108 """ 

109 

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

111 

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

113 

114 sol = spsolve(M, b) 

115 

116 me = self.dtype_u(self.init) 

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

118 

119 return me 

120 

121 def __eval_fexpl(self, u, t): 

122 """ 

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

124 

125 Parameters 

126 ---------- 

127 u : dtype_u 

128 Current values of the numerical solution. 

129 t : float 

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

131 

132 Returns 

133 ------- 

134 fexpl : dtype_f 

135 Explicit part of the right-hand side. 

136 """ 

137 

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

139 sol = self.Dx.dot(b) 

140 

141 fexpl = self.dtype_u(self.init) 

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

143 

144 return fexpl 

145 

146 def __eval_fimpl(self, u, t): 

147 """ 

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

149 

150 Parameters 

151 ---------- 

152 u : dtype_u 

153 Current values of the numerical solution. 

154 t : float 

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

156 

157 Returns 

158 ------- 

159 fimpl : dtype_f 

160 Implicit part of the right-hand side. 

161 """ 

162 

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

164 sol = self.A.dot(b) 

165 

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

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

168 

169 return fimpl 

170 

171 def eval_f(self, u, t): 

172 """ 

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

174 

175 Parameters 

176 ---------- 

177 u : dtype_u 

178 Current values of the numerical solution. 

179 t : float 

180 Current time of the numerical solution is computed. 

181 

182 Returns 

183 ------- 

184 f : dtype_f 

185 The right-hand side divided into two parts. 

186 """ 

187 

188 f = self.dtype_f(self.init) 

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

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

191 

192 self.work_counters['rhs']() 

193 return f 

194 

195 def u_exact(self, t): 

196 r""" 

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

198 

199 Parameters 

200 ---------- 

201 t : float 

202 Time of the exact solution. 

203 

204 Returns 

205 ------- 

206 me : dtype_u 

207 The exact solution. 

208 """ 

209 

210 def u_initial(x, k): 

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

212 

213 me = self.dtype_u(self.init) 

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

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

216 ) 

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

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

219 ) 

220 return me