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