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

## 58 statements

, created at 2024-09-09 14:59 +0000

1import numpy as np

3from pySDC.implementations.problem_classes.generic_ND_FD import GenericNDimFinDiff

4from pySDC.implementations.datatype_classes.mesh import imex_mesh

7class heatNd_unforced(GenericNDimFinDiff):

8 r"""

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

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)

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

19 .. math::

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

22 The spatial term is discretized using finite differences.

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

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 }

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

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

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

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

85 r"""

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

88 Parameters

89 ----------

90 t : float

91 Time of the exact solution.

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 )

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

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 )

132 return sol

135class heatNd_forced(heatNd_unforced):

136 r"""

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

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)

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

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

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

154 .. math::

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

157 The spatial term is discretized using finite differences.

158 """

160 dtype_f = imex_mesh

162 def eval_f(self, u, t):

163 """

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

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.

173 Returns

174 -------

175 f : dtype_f

176 The right-hand side of the problem.

177 """

179 f = self.f_init

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

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 )

204 return f

206 def u_exact(self, t):

207 r"""

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

210 Parameters

211 ----------

212 t : float

213 Time of the exact solution.

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