Coverage for pySDC/implementations/problem_classes/HeatEquation_ND_FD.py: 66%
58 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +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')
81 self._makeAttributeAndRegister('nu', localVars=locals(), readOnly=True)
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