Coverage for pySDC/implementations/problem_classes/Boussinesq_2D_FD_imex.py: 99%
67 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1import numpy as np
2from scipy.sparse.linalg import gmres
4from pySDC.core.problem import Problem
5from pySDC.implementations.datatype_classes.mesh import mesh, imex_mesh
6from pySDC.implementations.problem_classes.boussinesq_helpers.build2DFDMatrix import get2DMesh
7from pySDC.implementations.problem_classes.boussinesq_helpers.buildBoussinesq2DMatrix import getBoussinesq2DMatrix
8from pySDC.implementations.problem_classes.boussinesq_helpers.buildBoussinesq2DMatrix import getBoussinesq2DUpwindMatrix
9from pySDC.implementations.problem_classes.boussinesq_helpers.helper_classes import Callback, logging
10from pySDC.implementations.problem_classes.boussinesq_helpers.unflatten import unflatten
13# noinspection PyUnusedLocal
14class boussinesq_2d_imex(Problem):
15 r"""
16 Linearized 2D Boussinesq equations with finite differences, IMEX with waves implicit (GMRES), advection explicit.
18 The equations, for different boundary conditions, are
20 .. math::
21 \frac{\partial u}{\partial t} + U \frac{\partial u}{\partial x} + \frac{\partial p}{\partial x} = 0,
23 .. math::
24 \frac{\partial w}{\partial t} + U \frac{\partial w}{\partial x} + \frac{\partial p}{\partial z} = 0,
26 .. math::
27 \frac{\partial b}{\partial t} + U \frac{\partial b}{\partial x} + N^2 w = 0,
29 .. math::
30 \frac{\partial p}{\partial t} + U \frac{\partial p}{\partial x} + c^2 (\frac{\partial u}{\partial x} + \frac{\partial w}{\partial z}) = 0.
32 They can be derived from the linearized Euler equations by a transformation of variables [1]_.
34 Parameters
35 ----------
36 nvars : list of tuple, optional
37 List of number of unknowns nvars, e.g. ``nvars=[(4, 300, 3)]``.
38 c_s : float, optional
39 Acoustic velocity :math:`c_s`.
40 u_adv : float, optional
41 Advection speed :math:`U`.
42 Nfreq : float, optional
43 Stability frequency.
44 x_bounds : list, optional
45 Domain in x-direction.
46 z_bounds : list, optional
47 Domain in z-direction.
48 order_upw : int, optional
49 Order of upwind scheme for discretization.
50 order : int, optional
51 Order for discretization.
52 gmres_maxiter : int, optional
53 Maximum number of iterations for GMRES solver.
54 gmres_restart : int, optional
55 Number of iterations between restarts in GMRES solver.
56 gmres_tol_limit : float, optional
57 Tolerance for GMRES solver to terminate.
59 Attributes
60 ----------
61 N : list
62 List of number of unknowns.
63 bc_hor : list
64 Contains type of boundary conditions for both boundaries for both dimensions.
65 bc_ver :
66 Contains type of boundary conditions for both boundaries for both dimemsions, e.g. ``'neumann'`` or ``'dirichlet'``.
67 xx : np.ndarray
68 List of np.ndarrays for mesh in x-direction.
69 zz : np.ndarray
70 List of np.ndarrays for mesh in z-direction.
71 h : float
72 Mesh size.
73 Id : sp.sparse.eye
74 Identity matrix for the equation of appropriate size.
75 M : np.ndarray
76 Boussinesq 2D Matrix.
77 D_upwind : sp.csc_matrix
78 Boussinesq 2D Upwind matrix for discretization.
79 gmres_logger : object
80 Logger for GMRES solver.
82 References
83 ----------
84 .. [1] D. R. Durran. Numerical Methods for Fluid Dynamics. Texts Appl. Math. 32. Springer-Verlag, New York (2010)
85 http://dx.doi.org/10.1007/978-1-4419-6412-0.
86 """
88 dtype_u = mesh
89 dtype_f = imex_mesh
91 def __init__(
92 self,
93 nvars=None,
94 c_s=0.3,
95 u_adv=0.02,
96 Nfreq=0.01,
97 x_bounds=None,
98 z_bounds=None,
99 order_upw=5,
100 order=4,
101 gmres_maxiter=500,
102 gmres_restart=10,
103 gmres_tol_limit=1e-5,
104 ):
105 """Initialization routine"""
107 if nvars is None:
108 nvars = (4, 300, 30)
110 if x_bounds is None:
111 x_bounds = (-150.0, 150.0)
113 if z_bounds is None:
114 z_bounds = (0.0, 10.0)
116 # invoke super init, passing number of dofs, dtype_u and dtype_f
117 super().__init__(init=(nvars, None, np.dtype('float64')))
118 self._makeAttributeAndRegister(
119 'nvars',
120 'c_s',
121 'u_adv',
122 'Nfreq',
123 'x_bounds',
124 'z_bounds',
125 'order_upw',
126 'order',
127 'gmres_maxiter',
128 'gmres_restart',
129 localVars=locals(),
130 readOnly=True,
131 )
132 # not read-only: rungmrescounter_boussinesq tightens it for its reference solution
133 self._makeAttributeAndRegister('gmres_tol_limit', localVars=locals())
135 self.N = [self.nvars[1], self.nvars[2]]
137 self.bc_hor = [
138 ['periodic', 'periodic'],
139 ['periodic', 'periodic'],
140 ['periodic', 'periodic'],
141 ['periodic', 'periodic'],
142 ]
143 self.bc_ver = [
144 ['neumann', 'neumann'],
145 ['dirichlet', 'dirichlet'],
146 ['dirichlet', 'dirichlet'],
147 ['neumann', 'neumann'],
148 ]
150 self.xx, self.zz, self.h = get2DMesh(self.N, self.x_bounds, self.z_bounds, self.bc_hor[0], self.bc_ver[0])
152 self.Id, self.M = getBoussinesq2DMatrix(
153 self.N, self.h, self.bc_hor, self.bc_ver, self.c_s, self.Nfreq, self.order
154 )
155 self.D_upwind = getBoussinesq2DUpwindMatrix(self.N, self.h[0], self.u_adv, self.order_upw)
157 self.gmres_logger = logging()
159 def solve_system(self, rhs, factor, u0, t):
160 r"""
161 Simple linear solver for :math:`(I - factor \cdot A) \vec{u} = \vec{rhs}` using GMRES.
163 Parameters
164 ----------
165 rhs : dtype_f
166 Right-hand side for the nonlinear system.
167 factor : float
168 Abbrev. for the node-to-node stepsize (or any other factor required).
169 u0 : dtype_u
170 Initial guess for the iterative solver (not used here so far).
171 t : float
172 Current time (e.g. for time-dependent BCs).
174 Returns
175 -------
176 me : dtype_u
177 The solution as mesh.
178 """
180 b = rhs.flatten()
181 cb = Callback()
183 sol, info = gmres(
184 self.Id - factor * self.M,
185 b,
186 x0=u0.flatten(),
187 rtol=self.gmres_tol_limit,
188 restart=self.gmres_restart,
189 maxiter=self.gmres_maxiter,
190 atol=0,
191 callback=cb,
192 )
193 # If this is a dummy call with factor==0.0, do not log because it should not be counted as a solver call
194 if factor != 0.0:
195 self.gmres_logger.add(cb.getcounter())
196 me = self.dtype_u(self.init)
197 me[:] = unflatten(sol, 4, self.N[0], self.N[1])
199 return me
201 def __eval_fexpl(self, u, t):
202 """
203 Helper routine to evaluate the explicit part of the right-hand side.
205 Parameters
206 ----------
207 u : dtype_u
208 Current values of the numerical solution.
209 t : float
210 Current time at which the numerical solution is computed (not used here).
212 Returns
213 -------
214 fexpl : dtype_u
215 Explicit part of right-hand side.
216 """
218 # Evaluate right hand side
219 fexpl = self.dtype_u(self.init)
220 temp = u.flatten()
221 temp = self.D_upwind.dot(temp)
222 fexpl[:] = unflatten(temp, 4, self.N[0], self.N[1])
224 return fexpl
226 def __eval_fimpl(self, u, t):
227 """
228 Helper routine to evaluate the implicit part of the right-hand side.
230 Parameters
231 ----------
232 u : dtype_u
233 Current values of the numerical solution.
234 t : float
235 Current time at which the numerical solution is computed (not used here).
237 Returns
238 -------
239 fexpl : dtype_u
240 Implicit part of right-hand side.
241 """
243 temp = u.flatten()
244 temp = self.M.dot(temp)
245 fimpl = self.dtype_u(self.init)
246 fimpl[:] = unflatten(temp, 4, self.N[0], self.N[1])
248 return fimpl
250 def eval_f(self, u, t):
251 """
252 Routine to evaluate both parts of the right-hand side.
254 Parameters
255 ----------
256 u : dtype_u
257 Current values of the numerical solution.
258 t : float
259 Current time the numerical solution is computed.
261 Returns
262 -------
263 f : dtype_f
264 Right-hand side divided into two parts.
265 """
267 f = self.dtype_f(self.init)
268 f.impl[:] = self.__eval_fimpl(u, t)
269 f.expl[:] = self.__eval_fexpl(u, t)
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 me : dtype_u
284 The exact solution.
285 """
286 assert t == 0, 'ERROR: u_exact only valid for t=0'
288 dtheta = 0.01
289 H = 10.0
290 a = 5.0
291 x_c = -50.0
293 me = self.dtype_u(self.init)
294 me[0, :, :] = 0.0 * self.xx
295 me[1, :, :] = 0.0 * self.xx
296 # me[2,:,:] = 0.0*self.xx
297 # me[3,:,:] = np.exp(-0.5*(self.xx-0.0)**2.0/0.15**2.0)*np.exp(-0.5*(self.zz-0.5)**2/0.15**2)
298 # me[2,:,:] = np.exp(-0.5*(self.xx-0.0)**2.0/0.05**2.0)*np.exp(-0.5*(self.zz-0.5)**2/0.2**2)
299 me[2, :, :] = dtheta * np.sin(np.pi * self.zz / H) / (1.0 + np.square(self.xx - x_c) / (a * a))
300 me[3, :, :] = 0.0 * self.xx
301 return me