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

## 108 statements

, created at 2024-09-09 14:59 +0000

1import logging

3import dolfin as df

4import numpy as np

6from pySDC.core.problem import Problem

7from pySDC.implementations.datatype_classes.fenics_mesh import fenics_mesh, rhs_fenics_mesh

10# noinspection PyUnusedLocal

11class fenics_vortex_2d(Problem):

12 r"""

13 This class implements the vorticity-velocity problem in two dimensions with periodic boundary conditions

14 in :math:[0, 1]^2

16 .. math::

17 \frac{\partial w}{\partial t} = \nu \Delta w

19 for some parameter :math:\nu. In this class the problem is implemented that the spatial part is solved

20 using FEniCS [1]_. Hence, the problem is reformulated to the *weak formulation*

22 .. math::

23 \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx

25 This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.

26 The mass matrix needs inversion for this type of problem class, see the derived one for the mass-matrix version without inversion.

28 Parameters

29 ----------

30 c_nvars : List of int tuple, optional

31 Spatial resolution, i.e., numbers of degrees of freedom in space, e.g. c_nvars=[(128, 128)].

32 family : str, optional

33 Indicates the family of elements used to create the function space

34 for the trail and test functions. The default is 'CG', which are the class

35 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.

36 order : int, optional

37 Defines the order of the elements in the function space.

38 refinements : int, optional

39 Denotes the refinement of the mesh. refinements=2 refines the mesh by factor :math:2.

40 nu : float, optional

41 Diffusion coefficient :math:\nu.

42 rho : int, optional

43 Problem parameter.

44 delta : float, optional

45 Problem parameter.

47 Attributes

48 ----------

49 V : FunctionSpace

50 Defines the function space of the trial and test functions.

51 M : scalar, vector, matrix or higher rank tensor

52 Mass matrix for FENiCS.

53 K : scalar, vector, matrix or higher rank tensor

54 Stiffness matrix including diffusion coefficient (and correct sign).

56 References

57 ----------

58 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,

59 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).

60 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.

61 Wells and others. Springer (2012).

62 """

64 dtype_u = fenics_mesh

65 dtype_f = rhs_fenics_mesh

67 def __init__(self, c_nvars=None, family='CG', order=4, refinements=None, nu=0.01, rho=50, delta=0.05):

68 """

69 Initialization routine

71 Args:

72 problem_params (dict): custom parameters for the example

73 dtype_u: FEniCS mesh data type (will be passed to parent class)

74 dtype_f: FEniCS mesh data data type with implicit and explicit parts (will be passed to parent class)

75 """

77 if c_nvars is None:

78 c_nvars = [(32, 32)]

80 if refinements is None:

81 refinements = 1

83 # Subdomain for Periodic boundary condition

84 class PeriodicBoundary(df.SubDomain):

85 # Left boundary is "target domain" G

86 def inside(self, x, on_boundary):

87 # return True if on left or bottom boundary AND NOT on one of the two corners (0, 1) and (1, 0)

88 return bool(

89 (df.near(x[0], 0) or df.near(x[1], 0))

90 and (not ((df.near(x[0], 0) and df.near(x[1], 1)) or (df.near(x[0], 1) and df.near(x[1], 0))))

91 and on_boundary

92 )

94 def map(self, x, y):

95 if df.near(x[0], 1) and df.near(x[1], 1):

96 y[0] = x[0] - 1.0

97 y[1] = x[1] - 1.0

98 elif df.near(x[0], 1):

99 y[0] = x[0] - 1.0

100 y[1] = x[1]

101 else: # near(x[1], 1)

102 y[0] = x[0]

103 y[1] = x[1] - 1.0

105 # set logger level for FFC and dolfin

106 logging.getLogger('FFC').setLevel(logging.WARNING)

107 logging.getLogger('UFL').setLevel(logging.WARNING)

109 # set solver and form parameters

110 df.parameters["form_compiler"]["optimize"] = True

111 df.parameters["form_compiler"]["cpp_optimize"] = True

113 # set mesh and refinement (for multilevel)

114 mesh = df.UnitSquareMesh(c_nvars[0], c_nvars[1])

115 for _ in range(refinements):

116 mesh = df.refine(mesh)

118 self.mesh = df.Mesh(mesh)

120 # define function space for future reference

121 self.V = df.FunctionSpace(mesh, family, order, constrained_domain=PeriodicBoundary())

122 tmp = df.Function(self.V)

123 print('DoFs on this level:', len(tmp.vector()[:]))

124 self.fix_bc_for_residual = False

126 # invoke super init, passing number of dofs, dtype_u and dtype_f

127 super(fenics_vortex_2d, self).__init__(self.V)

128 self._makeAttributeAndRegister(

129 'c_nvars', 'family', 'order', 'refinements', 'nu', 'rho', 'delta', localVars=locals(), readOnly=True

130 )

132 w = df.TrialFunction(self.V)

133 v = df.TestFunction(self.V)

135 # Stiffness term (diffusion)

138 # Mass term

139 a_M = w * v * df.dx

141 self.M = df.assemble(a_M)

142 self.K = df.assemble(a_K)

144 def solve_system(self, rhs, factor, u0, t):

145 r"""

146 Dolfin's linear solver for :math:(M - factor \cdot A)\vec{u} = \vec{rhs}.

148 Parameters

149 ----------

150 rhs : dtype_f

151 Right-hand side for the nonlinear system.

152 factor : float

153 Abbrev. for the node-to-node stepsize (or any other factor required).

154 u0 : dtype_u

155 Initial guess for the iterative solver (not used here so far).

156 t : float

157 Current time.

159 Returns

160 -------

161 u : dtype_u

162 The solution as mesh.

163 """

165 A = self.M + self.nu * factor * self.K

166 b = self.apply_mass_matrix(rhs)

167 u = self.dtype_u(u0)

168 df.solve(A, u.values.vector(), b.values.vector())

170 return u

172 def __eval_fexpl(self, u, t):

173 """

174 Helper routine to evaluate the explicit part of the right-hand side.

176 Parameters

177 ----------

178 u : dtype_u

179 Current values of the numerical solution.

180 t : float

181 Current time at which the numerical solution is computed.

183 Returns

184 -------

185 fexpl : dtype_u

186 Explicit part of the right-hand side.

187 """

189 b = self.apply_mass_matrix(u)

190 psi = self.dtype_u(self.V)

191 df.solve(self.K, psi.values.vector(), b.values.vector())

193 fexpl = self.dtype_u(self.V)

194 fexpl.values = df.project(

195 df.Dx(psi.values, 1) * df.Dx(u.values, 0) - df.Dx(psi.values, 0) * df.Dx(u.values, 1), self.V

196 )

198 return fexpl

200 def __eval_fimpl(self, u, t):

201 """

202 Helper routine to evaluate the implicit part of the right-hand side.

204 Parameters

205 ----------

206 u : dtype_u

207 Current values of the numerical solution.

208 t : float

209 Current time at which the numerical solution is computed.

211 Returns

212 -------

213 fimpl : dtype_u

214 Implicit part of the right-hand side.

215 """

217 A = -self.nu * self.K

218 fimpl = self.dtype_u(self.V)

219 A.mult(u.values.vector(), fimpl.values.vector())

220 fimpl = self.__invert_mass_matrix(fimpl)

222 return fimpl

224 def eval_f(self, u, t):

225 """

226 Routine to evaluate both parts 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.

235 Returns

236 -------

237 f : dtype_f

238 The right-hand side divided into two parts.

239 """

241 f = self.dtype_f(self.V)

242 f.impl = self.__eval_fimpl(u, t)

243 f.expl = self.__eval_fexpl(u, t)

244 return f

246 def apply_mass_matrix(self, u):

247 r"""

248 Routine to apply mass matrix.

250 Parameters

251 u : dtype_u

252 Current values of the numerical solution.

254 Returns

255 -------

256 me : dtype_u

257 The product :math: M\vec{u}.

258 """

260 me = self.dtype_u(self.V)

261 self.M.mult(u.values.vector(), me.values.vector())

263 return me

265 def __invert_mass_matrix(self, u):

266 """

267 Helper routine to invert mass matrix.

269 Parameters

270 ----------

271 u : dtype_u

272 Current values of the numerical solution.

274 Returns

275 -------

276 me : dtype_u

277 The product :math:M^{-1} \vec{u}.

278 """

280 me = self.dtype_u(self.V)

281 df.solve(self.M, me.values.vector(), u.values.vector())

282 return me

284 def u_exact(self, t):

285 r"""

286 Routine to compute the exact solution at time :math:t.

288 Parameters

289 ----------

290 t : float

291 Time of the exact solution.

293 Returns

294 -------

295 me : dtype_u

296 The exact solution.

297 """

298 assert t == 0, 'ERROR: u_exact only valid for t=0'

300 w = df.Expression(

301 'r*(1-pow(tanh(r*((0.75-4) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-4))),2)) - \

302 r*(1-pow(tanh(r*((0.75-3) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-3))),2)) - \

303 r*(1-pow(tanh(r*((0.75-2) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-2))),2)) - \

304 r*(1-pow(tanh(r*((0.75-1) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-1))),2)) - \

305 r*(1-pow(tanh(r*((0.75-0) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25-0))),2)) - \

306 r*(1-pow(tanh(r*((0.75+1) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+1))),2)) - \

307 r*(1-pow(tanh(r*((0.75+2) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+2))),2)) - \

308 r*(1-pow(tanh(r*((0.75+3) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+3))),2)) - \

309 r*(1-pow(tanh(r*((0.75+4) - x[1])),2)) + r*(1-pow(tanh(r*(x[1] - (0.25+4))),2)) - \

310 d*2*a*cos(2*a*(x[0]+0.25))',

311 d=self.delta,

312 r=self.rho,

313 a=np.pi,

314 degree=self.order,

315 )

317 me = self.dtype_u(self.V)

318 me.values = df.interpolate(w, self.V)

320 # df.plot(me.values)

321 # df.interactive()

322 # exit()

324 return me

327class fenics_vortex_2d_mass(fenics_vortex_2d):

328 r"""

329 This class implements the vorticity-velocity problem in two dimensions with periodic boundary conditions

330 in :math:[0, 1]^2

332 .. math::

333 \frac{\partial w}{\partial t} = \nu \Delta w

335 for some parameter :math:\nu. In this class the problem is implemented that the spatial part is solved

336 using FEniCS [1]_. Hence, the problem is reformulated to the *weak formulation*

338 .. math::

339 \int_\Omega w_t v\,dx = - \nu \int_\Omega \nabla w \nabla v\,dx

341 This problem class treats the PDE in an IMEX way, with diffusion being the implicit part and everything else the explicit one.

342 No mass matrix inversion is needed here, i.e. using this problem class requires the imex_1st_order_mass sweeper.

344 Parameters

345 ----------

346 c_nvars : List of int tuple, optional

347 Spatial resolution, i.e., numbers of degrees of freedom in space, e.g. c_nvars=[(128, 128)].

348 family : str, optional

349 Indicates the family of elements used to create the function space

350 for the trail and test functions. The default is 'CG', which are the class

351 of Continuous Galerkin, a *synonym* for the Lagrange family of elements, see [2]_.

352 order : int, optional

353 Defines the order of the elements in the function space.

354 refinements : int, optional

355 Denotes the refinement of the mesh. refinements=2 refines the mesh by factor :math:2.

356 nu : float, optional

357 Diffusion coefficient :math:\nu.

358 rho : int, optional

359 Problem parameter.

360 delta : float, optional

361 Problem parameter.

363 Attributes

364 ----------

365 V : FunctionSpace

366 Defines the function space of the trial and test functions.

367 M : scalar, vector, matrix or higher rank tensor

368 Mass matrix for FENiCS.

369 K : scalar, vector, matrix or higher rank tensor

370 Stiffness matrix including diffusion coefficient (and correct sign).

372 References

373 ----------

374 .. [1] The FEniCS Project Version 1.5. M. S. Alnaes, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg,

375 C. Richardson, J. Ring, M. E. Rognes, G. N. Wells. Archive of Numerical Software (2015).

376 .. [2] Automated Solution of Differential Equations by the Finite Element Method. A. Logg, K.-A. Mardal, G. N.

377 Wells and others. Springer (2012).

378 """

380 def solve_system(self, rhs, factor, u0, t):

381 r"""

382 Dolfin's linear solver for :math:(M - factor \cdot A)\vec{u} = \vec{rhs}.

384 Parameters

385 ----------

386 rhs : dtype_f

387 Right-hand side for the nonlinear system.

388 factor : float

389 Abbrev. for the node-to-node stepsize (or any other factor required).

390 u0 : dtype_u

391 Initial guess for the iterative solver (not used here so far).

392 t : float

393 Current time.

395 Returns

396 -------

397 u : dtype_u

398 The solution as mesh.

399 """

401 A = self.M + self.nu * factor * self.K

402 b = self.dtype_u(rhs)

403 u = self.dtype_u(u0)

404 df.solve(A, u.values.vector(), b.values.vector())

406 return u

408 def __eval_fexpl(self, u, t):

409 """

410 Helper routine to evaluate the explicit part of the right-hand side.

412 Parameters

413 ----------

414 u : dtype_u

415 Current values of the numerical solution.

416 t : float

417 Current time at which the numerical solution is computed.

419 Returns

420 -------

421 fexpl : dtype_u

422 Explicit part of the right-hand side.

423 """

425 b = self.apply_mass_matrix(u)

426 psi = self.dtype_u(self.V)

427 df.solve(self.K, psi.values.vector(), b.values.vector())

429 fexpl = self.dtype_u(self.V)

430 fexpl.values = df.project(

431 df.Dx(psi.values, 1) * df.Dx(u.values, 0) - df.Dx(psi.values, 0) * df.Dx(u.values, 1), self.V

432 )

433 fexpl = self.apply_mass_matrix(fexpl)

435 return fexpl

437 def __eval_fimpl(self, u, t):

438 """

439 Helper routine to evaluate the implicit part of the right-hand side.

441 Parameters

442 ----------

443 u : dtype_u

444 Current values of the numerical solution.

445 t : float

446 Current time at which the numerical solution is computed.

448 Returns

449 -------

450 fimpl : dtype_u

451 Implicit part of the right-hand side.

452 """

454 A = -self.nu * self.K

455 fimpl = self.dtype_u(self.V)

456 A.mult(u.values.vector(), fimpl.values.vector())

458 return fimpl

460 def eval_f(self, u, t):

461 """

462 Routine to evaluate both parts of the right-hand side.

464 Note: Need to add this here, because otherwise the parent class will call the "local" functions __eval_*

465 and not the ones of the child class.

467 Parameters

468 ----------

469 u : dtype_u

470 Current values of the numerical solution.

471 t : float

472 Current time at which the numerical solution is computed.

474 Returns

475 -------

476 f : dtype_f

477 The right-hand side divided into two parts.

478 """

480 f = self.dtype_f(self.V)

481 f.impl = self.__eval_fimpl(u, t)

482 f.expl = self.__eval_fexpl(u, t)

483 return f