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 Example implementing the unforced one-dimensional advection diffusion equation 

12 

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} 

15 

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 

20 

21 .. math:: 

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

23 

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

25 discretization. 

26 

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. 

39 

40 Attributes 

41 ---------- 

42 xvalues : np.1darray 

43 Contains the grid points in space. 

44 ddx : np.1darray 

45 Spectral operator for gradient. 

46 lap : np.1darray 

47 Spectral operator for Laplacian. 

48 """ 

49 

50 dtype_u = mesh 

51 dtype_f = imex_mesh 

52 

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

54 """Initialization routine""" 

55 

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) 

59 

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

63 

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

65 

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 

69 

70 self.ddx = kx * 1j 

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

72 

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

74 

75 def eval_f(self, u, t): 

76 """ 

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

78 

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. 

85 

86 Returns 

87 ------- 

88 f : dtype_f 

89 The right-hand side of the problem. 

90 """ 

91 

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) 

98 

99 self.work_counters['rhs']() 

100 return f 

101 

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

103 """ 

104 Simple FFT solver for the diffusion part. 

105 

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

116 

117 Returns 

118 ------- 

119 me : dtype_u 

120 The solution as mesh. 

121 """ 

122 

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) 

126 

127 return me 

128 

129 def u_exact(self, t): 

130 r""" 

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

132 

133 Parameters 

134 ---------- 

135 t : float 

136 Time of the exact solution. 

137 

138 Returns 

139 ------- 

140 me : dtype_u 

141 The exact solution. 

142 """ 

143 

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 

162 

163 

164class advectiondiffusion1d_implicit(advectiondiffusion1d_imex): 

165 r""" 

166 Example implementing the unforced one-dimensional advection diffusion equation 

167 

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} 

170 

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 

173 

174 .. math:: 

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

176 

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

178 discretization. 

179 

180 Note 

181 ---- 

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

183 """ 

184 

185 dtype_u = mesh 

186 dtype_f = mesh 

187 

188 def eval_f(self, u, t): 

189 """ 

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

191 

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. 

198 

199 Returns 

200 ------- 

201 f : dtype_f 

202 The right-hand side of the problem. 

203 """ 

204 

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) 

209 

210 self.work_counters['rhs'] 

211 return f 

212 

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

214 """ 

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

216 

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

227 

228 Returns 

229 ------- 

230 me : dtype_u 

231 The solution as mesh. 

232 """ 

233 

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) 

237 

238 return me