## 52 statements

, created at 2024-09-20 17:10 +0000

1import numpy as np

2from scipy.sparse.linalg import spsolve

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,

10)

13# noinspection PyUnusedLocal

14class acoustic_1d_imex(Problem):

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

19 .. math::

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

22 .. math::

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

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

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),

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).

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

35 Parameters

36 ----------

37 nvars : int, optional

38 Number of degrees of freedom.

39 cs : float, optional

40 Sound velocity :math:c_s.

42 Advection speed :math:U.

44 Order of which the advective derivative is discretized.

45 waveno : int, optional

46 The wave number.

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.

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 """

67 dtype_u = mesh

68 dtype_f = imex_mesh

71 """Initialization routine"""

73 if nvars is None:

74 nvars = (2, 300)

76 # invoke super init, passing number of dofs

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

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

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

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

85 self.A = -self.cs * A

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

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}.

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).

104 Returns

105 -------

106 me : dtype_u

107 The solution as mesh.

108 """

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

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

114 sol = spsolve(M, b)

116 me = self.dtype_u(self.init)

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

119 return me

121 def __eval_fexpl(self, u, t):

122 """

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

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).

132 Returns

133 -------

134 fexpl : dtype_f

135 Explicit part of the right-hand side.

136 """

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

139 sol = self.Dx.dot(b)

141 fexpl = self.dtype_u(self.init)

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

144 return fexpl

146 def __eval_fimpl(self, u, t):

147 """

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

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).

157 Returns

158 -------

159 fimpl : dtype_f

160 Implicit part of the right-hand side.

161 """

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

164 sol = self.A.dot(b)

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

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

169 return fimpl

171 def eval_f(self, u, t):

172 """

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

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.

182 Returns

183 -------

184 f : dtype_f

185 The right-hand side divided into two parts.

186 """

188 f = self.dtype_f(self.init)

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

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

192 self.work_counters['rhs']()

193 return f

195 def u_exact(self, t):

196 r"""

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

199 Parameters

200 ----------

201 t : float

202 Time of the exact solution.

204 Returns

205 -------

206 me : dtype_u

207 The exact solution.

208 """

210 def u_initial(x, k):

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

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