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

1import numpy as np 

2 

3from pySDC.core.errors import ParameterError, ProblemError 

4from pySDC.core.problem import Problem, WorkCounter 

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

6 

7 

8# noinspection PyUnusedLocal 

9class advectiondiffusion1d_imex(Problem): 

10 r""" 

11 Periodic 1D advection-diffusion equation with FFTs, IMEX with diffusion implicit and advection explicit. 

12 

13 Example implementing the unforced one-dimensional advection diffusion equation 

14 

15 .. math:: 

16 \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} 

17 

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

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

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

21 given by 

22 

23 .. math:: 

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

25 

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

27 discretization. 

28 

29 Parameters 

30 ---------- 

31 nvars : int, optional 

32 Number of points in spatial discretization. 

33 c : float, optional 

34 Advection speed :math:`c`. 

35 freq : int, optional 

36 Wave number :math:`k`. 

37 nu : float, optional 

38 Diffusion coefficient :math:`\nu`. 

39 L : float, optional 

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

41 

42 Attributes 

43 ---------- 

44 xvalues : np.1darray 

45 Contains the grid points in space. 

46 ddx : np.1darray 

47 Spectral operator for gradient. 

48 lap : np.1darray 

49 Spectral operator for Laplacian. 

50 """ 

51 

52 dtype_u = mesh 

53 dtype_f = imex_mesh 

54 

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

56 """Initialization routine""" 

57 

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

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

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

61 

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

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

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

65 

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

67 

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

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

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

71 

72 self.ddx = kx * 1j 

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

74 

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

76 

77 def eval_f(self, u, t): 

78 """ 

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

80 

81 Parameters 

82 ---------- 

83 u : dtype_u 

84 Current values of the numerical solution. 

85 t : float 

86 Current time at which the numerical solution is computed. 

87 

88 Returns 

89 ------- 

90 f : dtype_f 

91 The right-hand side of the problem. 

92 """ 

93 

94 f = self.dtype_f(self.init) 

95 tmp_u = np.fft.rfft(u) 

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

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

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

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

100 

101 self.work_counters['rhs']() 

102 return f 

103 

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

105 """ 

106 Simple FFT solver for the diffusion part. 

107 

108 Parameters 

109 ---------- 

110 rhs : dtype_f 

111 Right-hand side for the linear system. 

112 factor : float 

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

114 u0 : dtype_u 

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

116 t : float 

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

118 

119 Returns 

120 ------- 

121 me : dtype_u 

122 The solution as mesh. 

123 """ 

124 

125 me = self.dtype_u(self.init) 

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

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

128 

129 return me 

130 

131 def u_exact(self, t): 

132 r""" 

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

134 

135 Parameters 

136 ---------- 

137 t : float 

138 Time of the exact solution. 

139 

140 Returns 

141 ------- 

142 me : dtype_u 

143 The exact solution. 

144 """ 

145 

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

147 if self.freq > 0: 

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

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

150 elif self.freq == 0: 

151 np.random.seed(1) 

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

153 else: 

154 t00 = 0.08 

155 if self.nu > 0: 

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

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

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

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

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

161 else: 

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

163 return me 

164 

165 

166class advectiondiffusion1d_implicit(advectiondiffusion1d_imex): 

167 r""" 

168 Periodic 1D advection-diffusion equation with FFTs, fully implicit in advection and diffusion. 

169 

170 Example implementing the unforced one-dimensional advection diffusion equation 

171 

172 .. math:: 

173 \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} 

174 

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

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

177 

178 .. math:: 

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

180 

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

182 discretization. 

183 

184 Note 

185 ---- 

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

187 """ 

188 

189 dtype_u = mesh 

190 dtype_f = mesh 

191 

192 def eval_f(self, u, t): 

193 """ 

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

195 

196 Parameters 

197 ---------- 

198 u : dtype_u 

199 Current values of the numerical solution. 

200 t : float 

201 Current time at which the numerical solution is computed. 

202 

203 Returns 

204 ------- 

205 f : dtype_f 

206 The right-hand side of the problem. 

207 """ 

208 

209 f = self.dtype_f(self.init) 

210 tmp_u = np.fft.rfft(u) 

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

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

213 

214 self.work_counters['rhs']() 

215 return f 

216 

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

218 """ 

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

220 

221 Parameters 

222 ---------- 

223 rhs : dtype_f 

224 Right-hand side for the linear system. 

225 factor : float 

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

227 u0 : dtype_u 

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

229 t : float 

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

231 

232 Returns 

233 ------- 

234 me : dtype_u 

235 The solution as mesh. 

236 """ 

237 

238 me = self.dtype_u(self.init) 

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

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

241 

242 return me