Coverage for pySDC/implementations/problem_classes/HeatEquation_ND_FD.py: 99%
68 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1from pySDC.implementations.problem_classes.generic_ND_FD import GenericNDimFinDiff
2from pySDC.implementations.datatype_classes.mesh import imex_mesh
5class heatNd_unforced(GenericNDimFinDiff):
6 r"""
7 This class implements the unforced :math:`N`-dimensional heat equation with periodic boundary conditions
9 .. math::
10 \frac{\partial u}{\partial t} = \nu
11 \left(
12 \frac{\partial^2 u}{\partial x^2_1} + .. + \frac{\partial^2 u}{\partial x^2_N}
13 \right)
15 for :math:`(x_1,..,x_N) \in [0, 1]^{N}` with :math:`N \leq 3`. The initial solution is of the form
17 .. math::
18 u({\bf x},0) = \prod_{i=1}^N \sin(\pi k_i x_i).
20 The spatial term is discretized using finite differences.
22 Parameters
23 ----------
24 nvars : int, optional
25 Spatial resolution (same in all dimensions). Using a tuple allows to
26 consider several dimensions, e.g ``nvars=(16,16)`` for a 2D problem.
27 nu : float, optional
28 Diffusion coefficient :math:`\nu`.
29 freq : int, optional
30 Spatial frequency :math:`k_i` of the initial conditions, can be tuple.
31 stencil_type : str, optional
32 Type of the finite difference stencil.
33 order : int, optional
34 Order of the finite difference discretization.
35 lintol : float, optional
36 Tolerance for spatial solver.
37 liniter : int, optional
38 Max. iterations number for spatial solver.
39 solver_type : str, optional
40 Solve the linear system directly or using CG.
41 bc : str, optional
42 Boundary conditions, either ``'periodic'`` or ``'dirichlet'``.
43 sigma : float, optional
44 If ``freq=-1`` and ``ndim=1``, uses a Gaussian initial solution of the form
46 .. math::
47 u(x,0) = e^{
48 \frac{\displaystyle 1}{\displaystyle 2}
49 \left(
50 \frac{\displaystyle x-1/2}{\displaystyle \sigma}
51 \right)^2
52 }
54 dtype : dtype-like, optional
55 Precision the state is stored at, ``float64`` by default. See :class:`GenericNDimFinDiff` for how the
56 operators follow it.
57 useGPU : bool, optional
58 Run on the GPU with CuPy instead of on the CPU with NumPy.
59 solve_precision : dtype-like or None, optional
60 Precision the implicit solve runs at, independently of ``dtype``. See :class:`GenericNDimFinDiff`.
62 Attributes
63 ----------
64 A : sparse matrix (CSC)
65 FD discretization matrix of the ND operator.
66 Id : sparse matrix (CSC)
67 Identity matrix of the same dimension as A
68 """
70 def __init__(
71 self,
72 nvars=512,
73 nu=0.1,
74 freq=2,
75 stencil_type='center',
76 order=2,
77 lintol=1e-12,
78 liniter=10000,
79 solver_type='direct',
80 bc='periodic',
81 sigma=6e-2,
82 dtype='float64',
83 useGPU=False,
84 solve_precision=None,
85 ):
86 """Initialization routine"""
87 super().__init__(
88 nvars,
89 nu,
90 2,
91 freq,
92 stencil_type,
93 order,
94 lintol,
95 liniter,
96 solver_type,
97 bc,
98 dtype=dtype,
99 useGPU=useGPU,
100 solve_precision=solve_precision,
101 )
102 if solver_type == 'GMRES':
103 self.logger.warning('GMRES is not usually used for heat equation')
104 self._makeAttributeAndRegister('nu', localVars=locals(), readOnly=True)
105 self._makeAttributeAndRegister('sigma', localVars=locals())
107 def u_exact(self, t, **kwargs):
108 r"""
109 Routine to compute the exact solution at time :math:`t`.
111 Parameters
112 ----------
113 t : float
114 Time of the exact solution.
116 Returns
117 -------
118 sol : dtype_u
119 The exact solution.
120 """
121 if 'u_init' in kwargs.keys() or 't_init' in kwargs.keys():
122 self.logger.warning(
123 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!'
124 )
126 ndim, freq, nu, sigma, dx, sol = self.ndim, self.freq, self.nu, self.sigma, self.dx, self.u_init
128 if ndim == 1:
129 x = self.grids
130 rho = (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2
131 if freq[0] > 0:
132 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.exp(-t * nu * rho)
133 elif freq[0] == -1: # Gaussian, spreading with the heat kernel, plus the periodic images that
134 # are still above machine precision at the current width
135 s = (sigma**2 + 2 * nu * t) ** 0.5
136 K = int(9 * s) + 1
137 sol[:] = sum(sigma / s * self.xp.exp(-0.5 * ((x - 0.5 - k) / s) ** 2) for k in range(-K, K + 1))
138 elif ndim == 2:
139 rho = (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2 + (
140 2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[1] * dx)
141 ) / dx**2
142 x, y = self.grids
143 sol[:] = (
144 self.xp.sin(self.xp.pi * freq[0] * x)
145 * self.xp.sin(self.xp.pi * freq[1] * y)
146 * self.xp.exp(-t * nu * rho)
147 )
148 elif ndim == 3:
149 rho = (
150 (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2
151 + (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[1] * dx)) / dx**2
152 + (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[2] * dx)) / dx**2
153 )
154 x, y, z = self.grids
155 sol[:] = (
156 self.xp.sin(self.xp.pi * freq[0] * x)
157 * self.xp.sin(self.xp.pi * freq[1] * y)
158 * self.xp.sin(self.xp.pi * freq[2] * z)
159 * self.xp.exp(-t * nu * rho)
160 )
162 return sol
165class heatNd_forced(heatNd_unforced):
166 r"""
167 This class implements the forced :math:`N`-dimensional heat equation with periodic boundary conditions
169 .. math::
170 \frac{\partial u}{\partial t} = \nu
171 \left(
172 \frac{\partial^2 u}{\partial x^2_1} + .. + \frac{\partial^2 u}{\partial x^2_N}
173 \right) + f({\bf x}, t)
175 for :math:`(x_1,..,x_N) \in [0, 1]^{N}` with :math:`N \leq 3`, and forcing term
177 .. math::
178 f({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \left(
179 \nu \pi^2 \sum_{i=1}^N k_i^2 \cos(t) - \sin(t)
180 \right),
182 where :math:`k_i` denotes the frequency in the :math:`i^{th}` dimension. The exact solution is
184 .. math::
185 u({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \cos(t).
187 The spatial term is discretized using finite differences.
188 """
190 dtype_f = imex_mesh
192 def setup_GPU(self):
193 """
194 Switch to GPU modules, keeping the split right-hand side this class needs
195 """
196 from pySDC.implementations.datatype_classes.cupy_mesh import imex_cupy_mesh
198 super().setup_GPU()
199 self.dtype_f = imex_cupy_mesh
201 def eval_f_increment(self, base, delta, t):
202 """
203 Evaluate the right-hand side increment, split the way :meth:`eval_f` splits it.
205 The forcing does not depend on ``u``, so the explicit part of the increment is zero. Without
206 this override the linear one inherited from :class:`GenericNDimFinDiff` would return an
207 unsplit right-hand side, which is the wrong type here and silently the wrong answer.
209 Parameters
210 ----------
211 base : dtype_u
212 The base state, unused: the implicit part is linear.
213 delta : dtype_u
214 The correction.
215 t : float
216 Current time, accepted for interface compatibility.
218 Returns
219 -------
220 f : dtype_f
221 The increment, with a zero explicit part.
222 """
223 f = self.dtype_f(self.init)
224 f.impl[:] = self.A.dot(delta.flatten()).reshape(self.nvars)
225 f.expl[:] = 0.0
226 return f
228 def eval_f(self, u, t):
229 """
230 Routine to evaluate the right-hand side of the problem.
232 Parameters
233 ----------
234 u : dtype_u
235 Current values of the numerical solution.
236 t : float
237 Current time of the numerical solution is computed.
239 Returns
240 -------
241 f : dtype_f
242 The right-hand side of the problem.
243 """
245 f = self.f_init
246 f.impl[:] = self.A.dot(u.flatten()).reshape(self.nvars)
248 ndim, freq, nu = self.ndim, self.freq, self.nu
249 if ndim == 1:
250 x = self.grids
251 f.expl[:] = self.xp.sin(self.xp.pi * freq[0] * x) * (
252 nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t)
253 )
254 elif ndim == 2:
255 x, y = self.grids
256 f.expl[:] = (
257 self.xp.sin(self.xp.pi * freq[0] * x)
258 * self.xp.sin(self.xp.pi * freq[1] * y)
259 * (nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t))
260 )
261 elif ndim == 3:
262 x, y, z = self.grids
263 f.expl[:] = (
264 self.xp.sin(self.xp.pi * freq[0] * x)
265 * self.xp.sin(self.xp.pi * freq[1] * y)
266 * self.xp.sin(self.xp.pi * freq[2] * z)
267 * (nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t))
268 )
270 return f
272 def u_exact(self, t):
273 r"""
274 Routine to compute the exact solution at time :math:`t`.
276 Parameters
277 ----------
278 t : float
279 Time of the exact solution.
281 Returns
282 -------
283 sol : dtype_u
284 The exact solution.
285 """
286 ndim, freq, sol = self.ndim, self.freq, self.u_init
287 if ndim == 1:
288 x = self.grids
289 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.cos(t)
290 elif ndim == 2:
291 x, y = self.grids
292 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.sin(self.xp.pi * freq[1] * y) * self.xp.cos(t)
293 elif ndim == 3:
294 x, y, z = self.grids
295 sol[:] = (
296 self.xp.sin(self.xp.pi * freq[0] * x)
297 * self.xp.sin(self.xp.pi * freq[1] * y)
298 * self.xp.sin(self.xp.pi * freq[2] * z)
299 * self.xp.cos(t)
300 )
301 return sol