Coverage for pySDC/implementations/problem_classes/HeatEquation_ND_FD.py: 88%
66 statements
« prev ^ index » next coverage.py v7.16.1, created at 2026-09-25 20:28 +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 Attributes
55 ----------
56 A : sparse matrix (CSC)
57 FD discretization matrix of the ND operator.
58 Id : sparse matrix (CSC)
59 Identity matrix of the same dimension as A
60 """
62 def __init__(
63 self,
64 nvars=512,
65 nu=0.1,
66 freq=2,
67 stencil_type='center',
68 order=2,
69 lintol=1e-12,
70 liniter=10000,
71 solver_type='direct',
72 bc='periodic',
73 sigma=6e-2,
74 dtype='float64',
75 useGPU=False,
76 ):
77 """Initialization routine"""
78 super().__init__(
79 nvars, nu, 2, freq, stencil_type, order, lintol, liniter, solver_type, bc, dtype=dtype, useGPU=useGPU
80 )
81 if solver_type == 'GMRES':
82 self.logger.warning('GMRES is not usually used for heat equation')
83 self._makeAttributeAndRegister('nu', localVars=locals(), readOnly=True)
84 self._makeAttributeAndRegister('sigma', localVars=locals())
86 def u_exact(self, t, **kwargs):
87 r"""
88 Routine to compute the exact solution at time :math:`t`.
90 Parameters
91 ----------
92 t : float
93 Time of the exact solution.
95 Returns
96 -------
97 sol : dtype_u
98 The exact solution.
99 """
100 if 'u_init' in kwargs.keys() or 't_init' in kwargs.keys():
101 self.logger.warning(
102 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!'
103 )
105 ndim, freq, nu, sigma, dx, sol = self.ndim, self.freq, self.nu, self.sigma, self.dx, self.u_init
107 if ndim == 1:
108 x = self.grids
109 rho = (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2
110 if freq[0] > 0:
111 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.exp(-t * nu * rho)
112 elif freq[0] == -1: # Gaussian
113 sol[:] = self.xp.exp(-0.5 * ((x - 0.5) / sigma) ** 2) * self.xp.exp(-t * nu * rho)
114 elif ndim == 2:
115 rho = (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2 + (
116 2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[1] * dx)
117 ) / dx**2
118 x, y = self.grids
119 sol[:] = (
120 self.xp.sin(self.xp.pi * freq[0] * x)
121 * self.xp.sin(self.xp.pi * freq[1] * y)
122 * self.xp.exp(-t * nu * rho)
123 )
124 elif ndim == 3:
125 rho = (
126 (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[0] * dx)) / dx**2
127 + (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[1] * dx)) / dx**2
128 + (2.0 - 2.0 * self.xp.cos(self.xp.pi * freq[2] * dx)) / dx**2
129 )
130 x, y, z = self.grids
131 sol[:] = (
132 self.xp.sin(self.xp.pi * freq[0] * x)
133 * self.xp.sin(self.xp.pi * freq[1] * y)
134 * self.xp.sin(self.xp.pi * freq[2] * z)
135 * self.xp.exp(-t * nu * rho)
136 )
138 return sol
141class heatNd_forced(heatNd_unforced):
142 r"""
143 This class implements the forced :math:`N`-dimensional heat equation with periodic boundary conditions
145 .. math::
146 \frac{\partial u}{\partial t} = \nu
147 \left(
148 \frac{\partial^2 u}{\partial x^2_1} + .. + \frac{\partial^2 u}{\partial x^2_N}
149 \right) + f({\bf x}, t)
151 for :math:`(x_1,..,x_N) \in [0, 1]^{N}` with :math:`N \leq 3`, and forcing term
153 .. math::
154 f({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \left(
155 \nu \pi^2 \sum_{i=1}^N k_i^2 \cos(t) - \sin(t)
156 \right),
158 where :math:`k_i` denotes the frequency in the :math:`i^{th}` dimension. The exact solution is
160 .. math::
161 u({\bf x}, t) = \prod_{i=1}^N \sin(\pi k_i x_i) \cos(t).
163 The spatial term is discretized using finite differences.
164 """
166 dtype_f = imex_mesh
168 def setup_GPU(self):
169 """
170 Switch to GPU modules, keeping the split right-hand side this class needs
171 """
172 from pySDC.implementations.datatype_classes.cupy_mesh import imex_cupy_mesh
174 super().setup_GPU()
175 self.dtype_f = imex_cupy_mesh
177 def eval_f_increment(self, base, delta, t):
178 """
179 Evaluate the right-hand side increment, split the way :meth:`eval_f` splits it.
181 The forcing does not depend on ``u``, so the explicit part of the increment is zero. Without
182 this override the linear one inherited from :class:`GenericNDimFinDiff` would return an
183 unsplit right-hand side, which is the wrong type here and silently the wrong answer.
185 Parameters
186 ----------
187 base : dtype_u
188 The base state, unused: the implicit part is linear.
189 delta : dtype_u
190 The correction.
191 t : float
192 Current time, accepted for interface compatibility.
194 Returns
195 -------
196 f : dtype_f
197 The increment, with a zero explicit part.
198 """
199 f = self.dtype_f(self.init)
200 f.impl[:] = self.A.dot(delta.flatten()).reshape(self.nvars)
201 f.expl[:] = 0.0
202 return f
204 def eval_f(self, u, t):
205 """
206 Routine to evaluate the right-hand side of the problem.
208 Parameters
209 ----------
210 u : dtype_u
211 Current values of the numerical solution.
212 t : float
213 Current time of the numerical solution is computed.
215 Returns
216 -------
217 f : dtype_f
218 The right-hand side of the problem.
219 """
221 f = self.f_init
222 f.impl[:] = self.A.dot(u.flatten()).reshape(self.nvars)
224 ndim, freq, nu = self.ndim, self.freq, self.nu
225 if ndim == 1:
226 x = self.grids
227 f.expl[:] = self.xp.sin(self.xp.pi * freq[0] * x) * (
228 nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t)
229 )
230 elif ndim == 2:
231 x, y = self.grids
232 f.expl[:] = (
233 self.xp.sin(self.xp.pi * freq[0] * x)
234 * self.xp.sin(self.xp.pi * freq[1] * y)
235 * (nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t))
236 )
237 elif ndim == 3:
238 x, y, z = self.grids
239 f.expl[:] = (
240 self.xp.sin(self.xp.pi * freq[0] * x)
241 * self.xp.sin(self.xp.pi * freq[1] * y)
242 * self.xp.sin(self.xp.pi * freq[2] * z)
243 * (nu * self.xp.pi**2 * sum([freq**2 for freq in freq]) * self.xp.cos(t) - self.xp.sin(t))
244 )
246 return f
248 def u_exact(self, t):
249 r"""
250 Routine to compute the exact solution at time :math:`t`.
252 Parameters
253 ----------
254 t : float
255 Time of the exact solution.
257 Returns
258 -------
259 sol : dtype_u
260 The exact solution.
261 """
262 ndim, freq, sol = self.ndim, self.freq, self.u_init
263 if ndim == 1:
264 x = self.grids
265 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.cos(t)
266 elif ndim == 2:
267 x, y = self.grids
268 sol[:] = self.xp.sin(self.xp.pi * freq[0] * x) * self.xp.sin(self.xp.pi * freq[1] * y) * self.xp.cos(t)
269 elif ndim == 3:
270 x, y, z = self.grids
271 sol[:] = (
272 self.xp.sin(self.xp.pi * freq[0] * x)
273 * self.xp.sin(self.xp.pi * freq[1] * y)
274 * self.xp.sin(self.xp.pi * freq[2] * z)
275 * self.xp.cos(t)
276 )
277 return sol