# Coverage for pySDC/implementations/problem_classes/Boussinesq_2D_FD_imex.py: 0%

## 66 statements

, created at 2024-09-20 17:10 +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,

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',

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(),

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:

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