Coverage for pySDC/implementations/problem_classes/Boussinesq_2D_FD_imex.py: 0%
66 statements
« prev ^ index » next coverage.py v7.6.9, created at 2024-12-20 14:51 +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 This class implements the two-dimensional Boussinesq equations for different boundary conditions with
18 .. math::
19 \frac{\partial u}{\partial t} + U \frac{\partial u}{\partial x} + \frac{\partial p}{\partial x} = 0,
21 .. math::
22 \frac{\partial w}{\partial t} + U \frac{\partial w}{\partial x} + \frac{\partial p}{\partial z} = 0,
24 .. math::
25 \frac{\partial b}{\partial t} + U \frac{\partial b}{\partial x} + N^2 w = 0,
27 .. math::
28 \frac{\partial p}{\partial t} + U \frac{\partial p}{\partial x} + c^2 (\frac{\partial u}{\partial x} + \frac{\partial w}{\partial z}) = 0.
30 They can be derived from the linearized Euler equations by a transformation of variables [1]_.
32 Parameters
33 ----------
34 nvars : list of tuple, optional
35 List of number of unknowns nvars, e.g. ``nvars=[(4, 300, 3)]``.
36 c_s : float, optional
37 Acoustic velocity :math:`c_s`.
38 u_adv : float, optional
39 Advection speed :math:`U`.
40 Nfreq : float, optional
41 Stability frequency.
42 x_bounds : list, optional
43 Domain in x-direction.
44 z_bounds : list, optional
45 Domain in z-direction.
46 order_upwind : int, optional
47 Order of upwind scheme for discretization.
48 order : int, optional
49 Order for discretization.
50 gmres_maxiter : int, optional
51 Maximum number of iterations for GMRES solver.
52 gmres_restart : int, optional
53 Number of iterations between restarts in GMRES solver.
54 gmres_tol_limit : float, optional
55 Tolerance for GMRES solver to terminate.
57 Attributes
58 ----------
59 N : list
60 List of number of unknowns.
61 bc_hor : list
62 Contains type of boundary conditions for both boundaries for both dimensions.
63 bc_ver :
64 Contains type of boundary conditions for both boundaries for both dimemsions, e.g. ``'neumann'`` or ``'dirichlet'``.
65 xx : np.ndarray
66 List of np.ndarrays for mesh in x-direction.
67 zz : np.ndarray
68 List of np.ndarrays for mesh in z-direction.
69 h : float
70 Mesh size.
71 Id : sp.sparse.eye
72 Identity matrix for the equation of appropriate size.
73 M : np.ndarray
74 Boussinesq 2D Matrix.
75 D_upwind : sp.csc_matrix
76 Boussinesq 2D Upwind matrix for discretization.
77 gmres_logger : object
78 Logger for GMRES solver.
80 References
81 ----------
82 .. [1] D. R. Durran. Numerical Methods for Fluid Dynamics. Texts Appl. Math. 32. Springer-Verlag, New York (2010)
83 http://dx.doi.org/10.1007/978-1-4419-6412-0.
84 """
86 dtype_u = mesh
87 dtype_f = imex_mesh
89 def __init__(
90 self,
91 nvars=None,
92 c_s=0.3,
93 u_adv=0.02,
94 Nfreq=0.01,
95 x_bounds=None,
96 z_bounds=None,
97 order_upw=5,
98 order=4,
99 gmres_maxiter=500,
100 gmres_restart=10,
101 gmres_tol_limit=1e-5,
102 ):
103 """Initialization routine"""
105 if nvars is None:
106 nvars = [(4, 300, 30)]
108 if x_bounds is None:
109 x_bounds = [(-150.0, 150.0)]
111 if z_bounds is None:
112 z_bounds = [(0.0, 10.0)]
114 # invoke super init, passing number of dofs, dtype_u and dtype_f
115 super().__init__(init=(nvars, None, np.dtype('float64')))
116 self._makeAttributeAndRegister(
117 'nvars',
118 'c_s',
119 'u_adv',
120 'Nfreq',
121 'x_bounds',
122 'z_bounds',
123 'order_upw',
124 'order',
125 'gmres_maxiter',
126 'gmres_restart',
127 'gmres_tol_limit',
128 localVars=locals(),
129 readOnly=True,
130 )
132 self.N = [self.nvars[1], self.nvars[2]]
134 self.bc_hor = [
135 ['periodic', 'periodic'],
136 ['periodic', 'periodic'],
137 ['periodic', 'periodic'],
138 ['periodic', 'periodic'],
139 ]
140 self.bc_ver = [
141 ['neumann', 'neumann'],
142 ['dirichlet', 'dirichlet'],
143 ['dirichlet', 'dirichlet'],
144 ['neumann', 'neumann'],
145 ]
147 self.xx, self.zz, self.h = get2DMesh(self.N, self.x_bounds, self.z_bounds, self.bc_hor[0], self.bc_ver[0])
149 self.Id, self.M = getBoussinesq2DMatrix(
150 self.N, self.h, self.bc_hor, self.bc_ver, self.c_s, self.Nfreq, self.order
151 )
152 self.D_upwind = getBoussinesq2DUpwindMatrix(self.N, self.h[0], self.u_adv, self.order_upw)
154 self.gmres_logger = logging()
156 def solve_system(self, rhs, factor, u0, t):
157 r"""
158 Simple linear solver for :math:`(I - factor \cdot A) \vec{u} = \vec{rhs}` using GMRES.
160 Parameters
161 ----------
162 rhs : dtype_f
163 Right-hand side for the nonlinear system.
164 factor : float
165 Abbrev. for the node-to-node stepsize (or any other factor required).
166 u0 : dtype_u
167 Initial guess for the iterative solver (not used here so far).
168 t : float
169 Current time (e.g. for time-dependent BCs).
171 Returns
172 -------
173 me : dtype_u
174 The solution as mesh.
175 """
177 b = rhs.flatten()
178 cb = Callback()
180 sol, info = gmres(
181 self.Id - factor * self.M,
182 b,
183 x0=u0.flatten(),
184 rtol=self.gmres_tol_limit,
185 restart=self.gmres_restart,
186 maxiter=self.gmres_maxiter,
187 atol=0,
188 callback=cb,
189 )
190 # If this is a dummy call with factor==0.0, do not log because it should not be counted as a solver call
191 if factor != 0.0:
192 self.gmres_logger.add(cb.getcounter())
193 me = self.dtype_u(self.init)
194 me[:] = unflatten(sol, 4, self.N[0], self.N[1])
196 return me
198 def __eval_fexpl(self, u, t):
199 """
200 Helper routine to evaluate the explicit part of the right-hand side.
202 Parameters
203 ----------
204 u : dtype_u
205 Current values of the numerical solution.
206 t : float
207 Current time at which the numerical solution is computed (not used here).
209 Returns
210 -------
211 fexpl : dtype_u
212 Explicit part of right-hand side.
213 """
215 # Evaluate right hand side
216 fexpl = self.dtype_u(self.init)
217 temp = u.flatten()
218 temp = self.D_upwind.dot(temp)
219 fexpl[:] = unflatten(temp, 4, self.N[0], self.N[1])
221 return fexpl
223 def __eval_fimpl(self, u, t):
224 """
225 Helper routine to evaluate the implicit part of the right-hand side.
227 Parameters
228 ----------
229 u : dtype_u
230 Current values of the numerical solution.
231 t : float
232 Current time at which the numerical solution is computed (not used here).
234 Returns
235 -------
236 fexpl : dtype_u
237 Implicit part of right-hand side.
238 """
240 temp = u.flatten()
241 temp = self.M.dot(temp)
242 fimpl = self.dtype_u(self.init)
243 fimpl[:] = unflatten(temp, 4, self.N[0], self.N[1])
245 return fimpl
247 def eval_f(self, u, t):
248 """
249 Routine to evaluate both parts of the right-hand side.
251 Parameters
252 ----------
253 u : dtype_u
254 Current values of the numerical solution.
255 t : float
256 Current time the numerical solution is computed.
258 Returns
259 -------
260 f : dtype_f
261 Right-hand side divided into two parts.
262 """
264 f = self.dtype_f(self.init)
265 f.impl[:] = self.__eval_fimpl(u, t)
266 f.expl[:] = self.__eval_fexpl(u, t)
267 return f
269 def u_exact(self, t):
270 r"""
271 Routine to compute the exact solution at time :math:`t`.
273 Parameters
274 ----------
275 t : float
276 Time of the exact solution.
278 Returns
279 -------
280 me : dtype_u
281 The exact solution.
282 """
283 assert t == 0, 'ERROR: u_exact only valid for t=0'
285 dtheta = 0.01
286 H = 10.0
287 a = 5.0
288 x_c = -50.0
290 me = self.dtype_u(self.init)
291 me[0, :, :] = 0.0 * self.xx
292 me[1, :, :] = 0.0 * self.xx
293 # me[2,:,:] = 0.0*self.xx
294 # 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)
295 # 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)
296 me[2, :, :] = dtheta * np.sin(np.pi * self.zz / H) / (1.0 + np.square(self.xx - x_c) / (a * a))
297 me[3, :, :] = 0.0 * self.xx
298 return me