Coverage for pySDC/implementations/problem_classes/HeatEquation_ND_FD.py: 66%

1import numpy as np 

2 

3from pySDC.implementations.problem_classes.generic_ND_FD import GenericNDimFinDiff 

4from pySDC.implementations.datatype_classes.mesh import imex_mesh 

5 

6 

7class heatNd_unforced(GenericNDimFinDiff): 

8 r""" 

9 This class implements the unforced :math:`N`-dimensional heat equation with periodic boundary conditions 

10 

11 .. math:: 

12 \frac{\partial u}{\partial t} = \nu 

13 \left( 

14 \frac{\partial^2 u}{\partial x^2_1} + .. + \frac{\partial^2 u}{\partial x^2_N} 

15 \right) 

16 

17 for :math:`(x_1,..,x_N) \in [0, 1]^{N}` with :math:`N \leq 3`. The initial solution is of the form 

18 

19 .. math:: 

20 u({\bf x},0) = \prod_{i=1}^N \sin(\pi k_i x_i). 

21 

22 The spatial term is discretized using finite differences. 

23 

24 Parameters 

25 ---------- 

26 nvars : int, optional 

27 Spatial resolution (same in all dimensions). Using a tuple allows to 

28 consider several dimensions, e.g ``nvars=(16,16)`` for a 2D problem. 

29 nu : float, optional 

30 Diffusion coefficient :math:`\nu`. 

31 freq : int, optional 

32 Spatial frequency :math:`k_i` of the initial conditions, can be tuple. 

33 stencil_type : str, optional 

34 Type of the finite difference stencil. 

35 order : int, optional 

36 Order of the finite difference discretization. 

37 lintol : float, optional 

38 Tolerance for spatial solver. 

39 liniter : int, optional 

40 Max. iterations number for spatial solver. 

41 solver_type : str, optional 

42 Solve the linear system directly or using CG. 

43 bc : str, optional 

44 Boundary conditions, either ``'periodic'`` or ``'dirichlet'``. 

45 sigma : float, optional 

46 If ``freq=-1`` and ``ndim=1``, uses a Gaussian initial solution of the form 

47 

48 .. math:: 

49 u(x,0) = e^{ 

50 \frac{\displaystyle 1}{\displaystyle 2} 

51 \left( 

52 \frac{\displaystyle x-1/2}{\displaystyle \sigma} 

53 \right)^2 

54 } 

55 

56 Attributes 

57 ---------- 

58 A : sparse matrix (CSC) 

59 FD discretization matrix of the ND operator. 

60 Id : sparse matrix (CSC) 

61 Identity matrix of the same dimension as A 

62 """ 

63 

64 def __init__( 

65 self, 

66 nvars=512, 

67 nu=0.1, 

68 freq=2, 

69 stencil_type='center', 

70 order=2, 

71 lintol=1e-12, 

72 liniter=10000, 

73 solver_type='direct', 

74 bc='periodic', 

75 sigma=6e-2, 

76 ): 

77 """Initialization routine""" 

78 super().__init__(nvars, nu, 2, freq, stencil_type, order, lintol, liniter, solver_type, bc) 

79 if solver_type == 'GMRES': 

80 self.logger.warning('GMRES is not usually used for heat equation') 

81 self._makeAttributeAndRegister('nu', localVars=locals(), readOnly=True) 

82 self._makeAttributeAndRegister('sigma', localVars=locals()) 

83 

84 def u_exact(self, t, **kwargs): 

85 r""" 

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

87 

88 Parameters 

89 ---------- 

90 t : float 

91 Time of the exact solution. 

92 

93 Returns 

94 ------- 

95 sol : dtype_u 

96 The exact solution. 

97 """ 

98 if 'u_init' in kwargs.keys() or 't_init' in kwargs.keys(): 

99 self.logger.warning( 

100 f'{type(self).__name__} uses an analytic exact solution from t=0. If you try to compute the local error, you will get the global error instead!' 

101 ) 

102 

103 ndim, freq, nu, sigma, dx, sol = self.ndim, self.freq, self.nu, self.sigma, self.dx, self.u_init 

104 

105 if ndim == 1: 

106 x = self.grids 

107 rho = (2.0 - 2.0 * np.cos(np.pi * freq[0] * dx)) / dx**2 

108 if freq[0] > 0: 

109 sol[:] = np.sin(np.pi * freq[0] * x) * np.exp(-t * nu * rho) 

110 elif freq[0] == -1: # Gaussian 

111 sol[:] = np.exp(-0.5 * ((x - 0.5) / sigma) ** 2) * np.exp(-t * nu * rho) 

112 elif ndim == 2: 

113 rho = (2.0 - 2.0 * np.cos(np.pi * freq[0] * dx)) / dx**2 + ( 

114 2.0 - 2.0 * np.cos(np.pi * freq[1] * dx) 

115 ) / dx**2 

116 x, y = self.grids 

117 sol[:] = np.sin(np.pi * freq[0] * x) * np.sin(np.pi * freq[1] * y) * np.exp(-t * nu * rho) 

118 elif ndim == 3: 

119 rho = ( 

120 (2.0 - 2.0 * np.cos(np.pi * freq[0] * dx)) / dx**2 

121 + (2.0 - 2.0 * np.cos(np.pi * freq[1] * dx)) 

122 + (2.0 - 2.0 * np.cos(np.pi * freq[2] * dx)) / dx**2 

123 ) 

124 x, y, z = self.grids 

125 sol[:] = ( 

126 np.sin(np.pi * freq[0] * x) 

127 * np.sin(np.pi * freq[1] * y) 

128 * np.sin(np.pi * freq[2] * z) 

129 * np.exp(-t * nu * rho) 

130 ) 

131 

132 return sol 

133 

134 

135class heatNd_forced(heatNd_unforced): 

136 r""" 

137 This class implements the forced :math:`N`-dimensional heat equation with periodic boundary conditions 

138 

139 .. math:: 

140 \frac{\partial u}{\partial t} = \nu 

141 \left( 

142 \frac{\partial^2 u}{\partial x^2_1} + .. + \frac{\partial^2 u}{\partial x^2_N} 

143 \right) + f({\bf x}, t) 

144 

145 for :math:`(x_1,..,x_N) \in [0, 1]^{N}` with :math:`N \leq 3`, and forcing term 

146 

147 .. math:: 

148 f({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \left( 

149 \nu \pi^2 \sum_{i=1}^N k_i^2 \cos(t) - \sin(t) 

150 \right), 

151 

152 where :math:`k_i` denotes the frequency in the :math:`i^{th}` dimension. The exact solution is 

153 

154 .. math:: 

155 u({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \cos(t). 

156 

157 The spatial term is discretized using finite differences. 

158 """ 

159 

160 dtype_f = imex_mesh 

161 

162 def eval_f(self, u, t): 

163 """ 

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

165 

166 Parameters 

167 ---------- 

168 u : dtype_u 

169 Current values of the numerical solution. 

170 t : float 

171 Current time of the numerical solution is computed. 

172 

173 Returns 

174 ------- 

175 f : dtype_f 

176 The right-hand side of the problem. 

177 """ 

178 

179 f = self.f_init 

180 f.impl[:] = self.A.dot(u.flatten()).reshape(self.nvars) 

181 

182 ndim, freq, nu = self.ndim, self.freq, self.nu 

183 if ndim == 1: 

184 x = self.grids 

185 f.expl[:] = np.sin(np.pi * freq[0] * x) * ( 

186 nu * np.pi**2 * sum([freq**2 for freq in freq]) * np.cos(t) - np.sin(t) 

187 ) 

188 elif ndim == 2: 

189 x, y = self.grids 

190 f.expl[:] = ( 

191 np.sin(np.pi * freq[0] * x) 

192 * np.sin(np.pi * freq[1] * y) 

193 * (nu * np.pi**2 * sum([freq**2 for freq in freq]) * np.cos(t) - np.sin(t)) 

194 ) 

195 elif ndim == 3: 

196 x, y, z = self.grids 

197 f.expl[:] = ( 

198 np.sin(np.pi * freq[0] * x) 

199 * np.sin(np.pi * freq[1] * y) 

200 * np.sin(np.pi * freq[2] * z) 

201 * (nu * np.pi**2 * sum([freq**2 for freq in freq]) * np.cos(t) - np.sin(t)) 

202 ) 

203 

204 return f 

205 

206 def u_exact(self, t): 

207 r""" 

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

209 

210 Parameters 

211 ---------- 

212 t : float 

213 Time of the exact solution. 

214 

215 Returns 

216 ------- 

217 sol : dtype_u 

218 The exact solution. 

219 """ 

220 ndim, freq, sol = self.ndim, self.freq, self.u_init 

221 if ndim == 1: 

222 x = self.grids 

223 sol[:] = np.sin(np.pi * freq[0] * x) * np.cos(t) 

224 elif ndim == 2: 

225 x, y = self.grids 

226 sol[:] = np.sin(np.pi * freq[0] * x) * np.sin(np.pi * freq[1] * y) * np.cos(t) 

227 elif ndim == 3: 

228 x, y, z = self.grids 

229 sol[:] = np.sin(np.pi * freq[0] * x) * np.sin(np.pi * freq[1] * y) * np.sin(np.pi * freq[2] * z) * np.cos(t) 

230 return sol