Coverage for pySDC/implementations/problem_classes/HeatEquation_ND_FD.py: 99%
68 statements
« prev ^ index » next coverage.py v7.17.3, created at 2026-11-12 18:41 +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 How to solve the linear systems: ``'direct'``, ``'GMRES'``, ``'CG'``, ``'FFT'`` or ``'MG'`` (multigrid).
41 See :class:`GenericNDimFinDiff`.
42 bc : str, optional
43 Boundary conditions, either ``'periodic'`` or ``'dirichlet'``.
44 sigma : float, optional
45 If ``freq=-1`` and ``ndim=1``, uses a Gaussian initial solution of the form
47 .. math::
48 u(x,0) = e^{
49 \frac{\displaystyle 1}{\displaystyle 2}
50 \left(
51 \frac{\displaystyle x-1/2}{\displaystyle \sigma}
52 \right)^2
53 }
55 dtype : dtype-like, optional
56 Precision the state is stored at, ``float64`` by default. See :class:`GenericNDimFinDiff` for how the
57 operators follow it.
58 useGPU : bool, optional
59 Run on the GPU with CuPy instead of on the CPU with NumPy.
60 solve_precision : dtype-like or None, optional
61 Precision the implicit solve runs at, independently of ``dtype``. See :class:`GenericNDimFinDiff`.
63 Attributes
64 ----------
65 A : sparse matrix (CSC)
66 FD discretization matrix of the ND operator.
67 Id : sparse matrix (CSC)
68 Identity matrix of the same dimension as A
69 """
71 def __init__(
72 self,
73 nvars=512,
74 nu=0.1,
75 freq=2,
76 stencil_type='center',
77 order=2,
78 lintol=1e-12,
79 liniter=10000,
80 solver_type='direct',
81 bc='periodic',
82 sigma=6e-2,
83 dtype='float64',
84 useGPU=False,
85 solve_precision=None,
86 ):
87 """Initialization routine"""
88 super().__init__(
89 nvars,
90 nu,
91 2,
92 freq,
93 stencil_type,
94 order,
95 lintol,
96 liniter,
97 solver_type,
98 bc,
99 dtype=dtype,
100 useGPU=useGPU,
101 solve_precision=solve_precision,
102 )
103 if solver_type == 'GMRES':
104 self.logger.warning('GMRES is not usually used for heat equation')
105 self._makeAttributeAndRegister('nu', localVars=locals(), readOnly=True)
106 self._makeAttributeAndRegister('sigma', localVars=locals())
108 def u_exact(self, t, **kwargs):
109 r"""
110 Routine to compute the exact solution at time :math:`t`.
112 Parameters
113 ----------
114 t : float
115 Time of the exact solution.
117 Returns
118 -------
119 sol : dtype_u
120 The exact solution.
121 """
122 if 'u_init' in kwargs.keys() or 't_init' in kwargs.keys():
123 self.logger.warning(
124 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!'
125 )
127 ndim, freq, nu, sigma, dx, sol = self.ndim, self.freq, self.nu, self.sigma, self.dx, self.u_init
129 if ndim == 1:
130 x = self.grids
131 rho = (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2
132 if freq[0] > 0:
133 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.exp(-t * nu * rho)
134 elif freq[0] == -1: # Gaussian, spreading with the heat kernel, plus the periodic images that
135 # are still above machine precision at the current width
136 s = (sigma**2 + 2 * nu * t) ** 0.5
137 K = int(9 * s) + 1
138 sol[:] = sum(sigma / s * self.xp.exp(-0.5 * ((x - 0.5 - k) / s) ** 2) for k in range(-K, K + 1))
139 elif ndim == 2:
140 rho = (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2 + (
141 2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[1] * dx)
142 ) / dx**2
143 x, y = self.grids
144 sol[:] = (
145 self.xp.sin(self.xp.pi * freq[0] * x)
146 * self.xp.sin(self.xp.pi * freq[1] * y)
147 * self.xp.exp(-t * nu * rho)
148 )
149 elif ndim == 3:
150 rho = (
151 (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2
152 + (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[1] * dx)) / dx**2
153 + (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[2] * dx)) / dx**2
154 )
155 x, y, z = self.grids
156 sol[:] = (
157 self.xp.sin(self.xp.pi * freq[0] * x)
158 * self.xp.sin(self.xp.pi * freq[1] * y)
159 * self.xp.sin(self.xp.pi * freq[2] * z)
160 * self.xp.exp(-t * nu * rho)
161 )
163 return sol
166class heatNd_forced(heatNd_unforced):
167 r"""
168 This class implements the forced :math:`N`-dimensional heat equation with periodic boundary conditions
170 .. math::
171 \frac{\partial u}{\partial t} = \nu
172 \left(
173 \frac{\partial^2 u}{\partial x^2_1} + .. + \frac{\partial^2 u}{\partial x^2_N}
174 \right) + f({\bf x}, t)
176 for :math:`(x_1,..,x_N) \in [0, 1]^{N}` with :math:`N \leq 3`, and forcing term
178 .. math::
179 f({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \left(
180 \nu \pi^2 \sum_{i=1}^N k_i^2 \cos(t) - \sin(t)
181 \right),
183 where :math:`k_i` denotes the frequency in the :math:`i^{th}` dimension. The exact solution is
185 .. math::
186 u({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \cos(t).
188 The spatial term is discretized using finite differences.
189 """
191 dtype_f = imex_mesh
193 def setup_GPU(self):
194 """
195 Switch to GPU modules, keeping the split right-hand side this class needs
196 """
197 from pySDC.implementations.datatype_classes.cupy_mesh import imex_cupy_mesh
199 super().setup_GPU()
200 self.dtype_f = imex_cupy_mesh
202 def eval_f_increment(self, base, delta, t):
203 """
204 Evaluate the right-hand side increment, split the way :meth:`eval_f` splits it.
206 The forcing does not depend on ``u``, so the explicit part of the increment is zero. Without
207 this override the linear one inherited from :class:`GenericNDimFinDiff` would return an
208 unsplit right-hand side, which is the wrong type here and silently the wrong answer.
210 Parameters
211 ----------
212 base : dtype_u
213 The base state, unused: the implicit part is linear.
214 delta : dtype_u
215 The correction.
216 t : float
217 Current time, accepted for interface compatibility.
219 Returns
220 -------
221 f : dtype_f
222 The increment, with a zero explicit part.
223 """
224 f = self.dtype_f(self.init)
225 f.impl[:] = self.A.dot(delta.flatten()).reshape(self.nvars)
226 f.expl[:] = 0.0
227 return f
229 def eval_f(self, u, t):
230 """
231 Routine to evaluate the right-hand side of the problem.
233 Parameters
234 ----------
235 u : dtype_u
236 Current values of the numerical solution.
237 t : float
238 Current time of the numerical solution is computed.
240 Returns
241 -------
242 f : dtype_f
243 The right-hand side of the problem.
244 """
246 f = self.f_init
247 f.impl[:] = self.A.dot(u.flatten()).reshape(self.nvars)
249 ndim, freq, nu = self.ndim, self.freq, self.nu
250 if ndim == 1:
251 x = self.grids
252 f.expl[:] = self.xp.sin(self.xp.pi * freq[0] * x) * (
253 nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t)
254 )
255 elif ndim == 2:
256 x, y = self.grids
257 f.expl[:] = (
258 self.xp.sin(self.xp.pi * freq[0] * x)
259 * self.xp.sin(self.xp.pi * freq[1] * y)
260 * (nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t))
261 )
262 elif ndim == 3:
263 x, y, z = self.grids
264 f.expl[:] = (
265 self.xp.sin(self.xp.pi * freq[0] * x)
266 * self.xp.sin(self.xp.pi * freq[1] * y)
267 * self.xp.sin(self.xp.pi * freq[2] * z)
268 * (nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t))
269 )
271 return f
273 def u_exact(self, t):
274 r"""
275 Routine to compute the exact solution at time :math:`t`.
277 Parameters
278 ----------
279 t : float
280 Time of the exact solution.
282 Returns
283 -------
284 sol : dtype_u
285 The exact solution.
286 """
287 ndim, freq, sol = self.ndim, self.freq, self.u_init
288 if ndim == 1:
289 x = self.grids
290 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.cos(t)
291 elif ndim == 2:
292 x, y = self.grids
293 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.sin(self.xp.pi * freq[1] * y) * self.xp.cos(t)
294 elif ndim == 3:
295 x, y, z = self.grids
296 sol[:] = (
297 self.xp.sin(self.xp.pi * freq[0] * x)
298 * self.xp.sin(self.xp.pi * freq[1] * y)
299 * self.xp.sin(self.xp.pi * freq[2] * z)
300 * self.xp.cos(t)
301 )
302 return sol