## 63 statements

, created at 2024-09-20 16:55 +0000

1import numpy as np

3from pySDC.core.errors import ParameterError, ProblemError

4from pySDC.core.problem import Problem, WorkCounter

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

8# noinspection PyUnusedLocal

10 r"""

11 Example implementing the unforced one-dimensional advection diffusion equation

13 .. math::

14 \frac{\partial u(x,t)}{\partial t} = - c \frac{\partial u(x,t)}{\partial x} + \nu \frac{\partial^2 u(x,t)}{\partial x^2}

16 with periodic boundary conditions in :math:[-\frac{L}{2}, \frac{L}{2}] in spectral space. The advection part

17 :math:- c \frac{\partial u(x,t)}{\partial x} is treated explicitly, whereas the diffusion part

18 :math:\nu \frac{\partial^2 u(x,t)}{\partial x^2} will be treated numerically in an implicit way. The exact solution is

19 given by

21 .. math::

22 u(x, t) = \sin(\omega (x - c t)) \exp(-t \nu \omega^2)

24 for :math:\omega=2 \pi k, where :math:k denotes the wave number. Fast Fourier transform is used for the spatial

25 discretization.

27 Parameters

28 ----------

29 nvars : int, optional

30 Number of points in spatial discretization.

31 c : float, optional

32 Advection speed :math:c.

33 freq : int, optional

34 Wave number :math:k.

35 nu : float, optional

36 Diffusion coefficient :math:\nu.

37 L : float, optional

38 Denotes the period of the function to be approximated for the Fourier transform.

40 Attributes

41 ----------

42 xvalues : np.1darray

43 Contains the grid points in space.

44 ddx : np.1darray

46 lap : np.1darray

47 Spectral operator for Laplacian.

48 """

50 dtype_u = mesh

51 dtype_f = imex_mesh

53 def __init__(self, nvars=256, c=1.0, freq=-1, nu=0.02, L=1.0):

54 """Initialization routine"""

56 # invoke super init, passing number of dofs, dtype_u and dtype_f

57 super().__init__(init=(nvars, None, np.dtype('float64')))

58 self._makeAttributeAndRegister('nvars', 'c', 'freq', 'nu', 'L', localVars=locals(), readOnly=True)

60 # we assert that nvars looks very particular here.. this will be necessary for coarsening in space later on

61 if (self.nvars) % 2 != 0:

62 raise ProblemError('setup requires nvars = 2^p')

64 self.xvalues = np.array([i * self.L / self.nvars - self.L / 2.0 for i in range(self.nvars)])

66 kx = np.zeros(self.init[0] // 2 + 1)

67 for i in range(0, len(kx)):

68 kx[i] = 2 * np.pi / self.L * i

70 self.ddx = kx * 1j

71 self.lap = -(kx**2)

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

75 def eval_f(self, u, t):

76 """

77 Routine to evaluate the right-hand side of the problem.

79 Parameters

80 ----------

81 u : dtype_u

82 Current values of the numerical solution.

83 t : float

84 Current time at which the numerical solution is computed.

86 Returns

87 -------

88 f : dtype_f

89 The right-hand side of the problem.

90 """

92 f = self.dtype_f(self.init)

93 tmp_u = np.fft.rfft(u)

94 tmp_impl = self.nu * self.lap * tmp_u

95 tmp_expl = -self.c * self.ddx * tmp_u

96 f.impl[:] = np.fft.irfft(tmp_impl)

97 f.expl[:] = np.fft.irfft(tmp_expl)

99 self.work_counters['rhs']()

100 return f

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

103 """

104 Simple FFT solver for the diffusion part.

106 Parameters

107 ----------

108 rhs : dtype_f

109 Right-hand side for the linear system.

110 factor : float

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

112 u0 : dtype_u

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

114 t : float

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

117 Returns

118 -------

119 me : dtype_u

120 The solution as mesh.

121 """

123 me = self.dtype_u(self.init)

124 tmp = np.fft.rfft(rhs) / (1.0 - self.nu * factor * self.lap)

125 me[:] = np.fft.irfft(tmp)

127 return me

129 def u_exact(self, t):

130 r"""

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

133 Parameters

134 ----------

135 t : float

136 Time of the exact solution.

138 Returns

139 -------

140 me : dtype_u

141 The exact solution.

142 """

144 me = self.dtype_u(self.init, val=0.0)

145 if self.freq > 0:

146 omega = 2.0 * np.pi * self.freq

147 me[:] = np.sin(omega * (self.xvalues - self.c * t)) * np.exp(-t * self.nu * omega**2)

148 elif self.freq == 0:

149 np.random.seed(1)

150 me[:] = np.random.rand(self.nvars)

151 else:

152 t00 = 0.08

153 if self.nu > 0:

154 nbox = int(np.ceil(np.sqrt(4.0 * self.nu * (t00 + t) * 37.0 / (self.L**2))))

155 for k in range(-nbox, nbox + 1):

156 for i in range(self.init[0]):

157 x = self.xvalues[i] - self.c * t + k * self.L

158 me[i] += np.sqrt(t00) / np.sqrt(t00 + t) * np.exp(-(x**2) / (4.0 * self.nu * (t00 + t)))

159 else:

160 raise ParameterError("There is no exact solution implemented for negative frequency and negative nu!")

161 return me

165 r"""

166 Example implementing the unforced one-dimensional advection diffusion equation

168 .. math::

169 \frac{\partial u(x,t)}{\partial t} = - c \frac{\partial u(x,t)}{\partial x} + \nu \frac{\partial^2 u(x,t)}{\partial x^2}

171 with periodic boundary conditions in :math:[-\frac{L}{2}, \frac{L}{2}] in spectral space. This class implements the

172 problem solving it with fully-implicit time-stepping. The exact solution is given by

174 .. math::

175 u(x, t) = \sin(\omega (x - c t)) \exp(-t \nu \omega^2)

177 for :math:\omega=2 \pi k, where :math:k denotes the wave number. Fast Fourier transform is used for the spatial

178 discretization.

180 Note

181 ----

182 This class has the same attributes as the class it inherits from.

183 """

185 dtype_u = mesh

186 dtype_f = mesh

188 def eval_f(self, u, t):

189 """

190 Routine to evaluate the right-hand side of the problem.

192 Parameters

193 ----------

194 u : dtype_u

195 Current values of the numerical solution.

196 t : float

197 Current time at which the numerical solution is computed.

199 Returns

200 -------

201 f : dtype_f

202 The right-hand side of the problem.

203 """

205 f = self.dtype_f(self.init)

206 tmp_u = np.fft.rfft(u)

207 tmp = self.nu * self.lap * tmp_u - self.c * self.ddx * tmp_u

208 f[:] = np.fft.irfft(tmp)

210 self.work_counters['rhs']

211 return f

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

214 """

215 Simple FFT solver for the diffusion and advection part (both are linear!).

217 Parameters

218 ----------

219 rhs : dtype_f

220 Right-hand side for the linear system.

221 factor : float

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

223 u0 : dtype_u

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

225 t : float

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

228 Returns

229 -------

230 me : dtype_u

231 The solution as mesh.

232 """

234 me = self.dtype_u(self.init)

235 tmp = np.fft.rfft(rhs) / (1.0 - factor * (self.nu * self.lap - self.c * self.ddx))

236 me[:] = np.fft.irfft(tmp)

238 return me