# Coverage for pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix_forced.py: 100%

## 112 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_heat(Problem):

12 r"""

13 Example implementing the forced one-dimensional heat equation with Dirichlet boundary conditions

15 .. math::

16 \frac{d u}{d t} = \nu \frac{d^2 u}{d x^2} + f

18 for :math:x \in \Omega:=[0,1], where the forcing term :math:f is defined by

20 .. math::

21 f(x, t) = -\sin(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).

23 For initial conditions with constant c and

25 .. math::

26 u(x, 0) = \sin(\pi x) + c

28 the exact solution of the problem is given by

30 .. math::

31 u(x, t) = \sin(\pi x)\cos(t) + c.

33 In this class the problem is implemented in the way that the spatial part is solved using FEniCS [1]_. Hence, the problem

34 is reformulated to the *weak formulation*

36 .. math:

37 \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.

39 The part containing the forcing term is treated explicitly, where it is interpolated in the function space.

40 The other part will be treated in an implicit way.

42 Parameters

43 ----------

44 c_nvars : int, optional

45 Spatial resolution, i.e., numbers of degrees of freedom in space.

46 t0 : float, optional

47 Starting time.

48 family : str, optional

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

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

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

52 order : int, optional

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

54 refinements : int, optional

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

56 nu : float, optional

57 Diffusion coefficient :math:\nu.

58 c: float, optional

59 Constant for the Dirichlet boundary condition :math: c

61 Attributes

62 ----------

63 V : FunctionSpace

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

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

66 Denotes the expression :math:\int_\Omega u_t v\,dx.

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

68 Denotes the expression :math:- \nu \int_\Omega \nabla u \nabla v\,dx.

69 g : Expression

70 The forcing term :math:f in the heat equation.

71 bc : DirichletBC

72 Denotes the Dirichlet boundary conditions.

74 References

75 ----------

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

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

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

79 Wells and others. Springer (2012).

80 """

82 dtype_u = fenics_mesh

83 dtype_f = rhs_fenics_mesh

85 def __init__(self, c_nvars=128, t0=0.0, family='CG', order=4, refinements=1, nu=0.1, c=0.0):

86 """Initialization routine"""

88 # define the Dirichlet boundary

89 def Boundary(x, on_boundary):

90 return on_boundary

92 # set logger level for FFC and dolfin

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

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

96 # set solver and form parameters

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

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

99 df.parameters['allow_extrapolation'] = True

101 # set mesh and refinement (for multilevel)

102 mesh = df.UnitIntervalMesh(c_nvars)

103 for _ in range(refinements):

104 mesh = df.refine(mesh)

106 # define function space for future reference

107 self.V = df.FunctionSpace(mesh, family, order)

108 tmp = df.Function(self.V)

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

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

112 super(fenics_heat, self).__init__(self.V)

113 self._makeAttributeAndRegister(

114 'c_nvars', 't0', 'family', 'order', 'refinements', 'nu', 'c', localVars=locals(), readOnly=True

115 )

117 # Stiffness term (Laplace)

118 u = df.TrialFunction(self.V)

119 v = df.TestFunction(self.V)

122 # Mass term

123 a_M = u * v * df.dx

125 self.M = df.assemble(a_M)

126 self.K = df.assemble(a_K)

128 # set boundary values

129 self.bc = df.DirichletBC(self.V, df.Constant(c), Boundary)

130 self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)

132 # set forcing term as expression

133 self.g = df.Expression(

134 '-sin(a*x[0]) * (sin(t) - b*a*a*cos(t))',

135 a=np.pi,

136 b=self.nu,

137 t=self.t0,

138 degree=self.order,

139 )

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

142 r"""

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

145 Parameters

146 ----------

147 rhs : dtype_f

148 Right-hand side for the nonlinear system.

149 factor : float

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

151 u0 : dtype_u

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

153 t : float

154 Current time.

156 Returns

157 -------

158 u : dtype_u

159 Solution.

160 """

162 b = self.apply_mass_matrix(rhs)

164 u = self.dtype_u(u0)

165 T = self.M - factor * self.K

166 self.bc.apply(T, b.values.vector())

167 df.solve(T, u.values.vector(), b.values.vector())

169 return u

171 def __eval_fexpl(self, u, t):

172 """

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

175 Parameters

176 ----------

177 u : dtype_u

178 Current values of the numerical solution (not used here).

179 t : float

180 Current time at which the numerical solution is computed.

182 Returns

183 -------

184 fexpl : dtype_u

185 Explicit part of the right-hand side.

186 """

188 self.g.t = t

189 fexpl = self.dtype_u(df.interpolate(self.g, self.V))

191 return fexpl

193 def __eval_fimpl(self, u, t):

194 """

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

197 Parameters

198 ----------

199 u : dtype_u

200 Current values of the numerical solution.

201 t : float

202 Current time at which the numerical solution is computed.

204 Returns

205 -------

206 fimpl : dtype_u

207 Explicit part of the right-hand side.

208 """

210 tmp = self.dtype_u(self.V)

211 self.K.mult(u.values.vector(), tmp.values.vector())

212 fimpl = self.__invert_mass_matrix(tmp)

214 return fimpl

216 def eval_f(self, u, t):

217 """

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

220 Parameters

221 ----------

222 u : dtype_u

223 Current values of the numerical solution.

224 t : float

225 Current time at which the numerical solution is computed.

227 Returns

228 -------

229 f : dtype_f

230 The right-hand side divided into two parts.

231 """

233 f = self.dtype_f(self.V)

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

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

236 return f

238 def apply_mass_matrix(self, u):

239 r"""

240 Routine to apply mass matrix.

242 Parameters

243 ----------

244 u : dtype_u

245 Current values of the numerical solution.

247 Returns

248 -------

249 me : dtype_u

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

251 """

253 me = self.dtype_u(self.V)

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

256 return me

258 def __invert_mass_matrix(self, u):

259 r"""

260 Helper routine to invert mass matrix.

262 Parameters

263 ----------

264 u : dtype_u

265 Current values of the numerical solution.

267 Returns

268 -------

269 me : dtype_u

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

271 """

273 me = self.dtype_u(self.V)

275 b = self.dtype_u(u)

276 M = self.M

277 self.bc_hom.apply(M, b.values.vector())

279 df.solve(M, me.values.vector(), b.values.vector())

280 return me

282 def u_exact(self, t):

283 r"""

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

286 Parameters

287 ----------

288 t : float

289 Time of the exact solution.

291 Returns

292 -------

293 me : dtype_u

294 Exact solution.

295 """

297 u0 = df.Expression('sin(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t, degree=self.order)

298 me = self.dtype_u(df.interpolate(u0, self.V), val=self.V)

300 return me

303# noinspection PyUnusedLocal

304class fenics_heat_mass(fenics_heat):

305 r"""

306 Example implementing the forced one-dimensional heat equation with Dirichlet boundary conditions

308 .. math::

309 \frac{d u}{d t} = \nu \frac{d^2 u}{d x^2} + f

311 for :math:x \in \Omega:=[0,1], where the forcing term :math:f is defined by

313 .. math::

314 f(x, t) = -\sin(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).

316 For initial conditions with constant c and

318 .. math::

319 u(x, 0) = \sin(\pi x) + c

321 the exact solution of the problem is given by

323 .. math::

324 u(x, t) = \sin(\pi x)\cos(t) + c.

326 In this class the problem is implemented in the way that the spatial part is solved using FEniCS [1]_. Hence, the problem

327 is reformulated to the *weak formulation*

329 .. math:

330 \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.

332 The forcing term is treated explicitly, and is expressed via the mass matrix resulting from the left-hand side term

333 :math:\int_\Omega u_t v\,dx, and the other part will be treated in an implicit way.

335 Parameters

336 ----------

337 c_nvars : int, optional

338 Spatial resolution, i.e., numbers of degrees of freedom in space.

339 t0 : float, optional

340 Starting time.

341 family : str, optional

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

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

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

345 order : int, optional

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

347 refinements : int, optional

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

349 nu : float, optional

350 Diffusion coefficient :math:\nu.

352 Attributes

353 ----------

354 V : FunctionSpace

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

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

357 Denotes the expression :math:\int_\Omega u_t v\,dx.

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

359 Denotes the expression :math:- \nu \int_\Omega \nabla u \nabla v\,dx.

360 g : Expression

361 The forcing term :math:f in the heat equation.

362 bc : DirichletBC

363 Denotes the Dirichlet boundary conditions.

364 bc_hom : DirichletBC

365 Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual

366 fix_bc_for_residual: boolean

367 flag to indicate that the residual requires special treatment due to boundary conditions

369 References

370 ----------

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

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

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

374 Wells and others. Springer (2012).

375 """

377 def __init__(self, c_nvars=128, t0=0.0, family='CG', order=4, refinements=1, nu=0.1, c=0.0):

378 """Initialization routine"""

380 super().__init__(c_nvars, t0, family, order, refinements, nu, c)

382 self.fix_bc_for_residual = True

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

385 r"""

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

388 Parameters

389 ----------

390 rhs : dtype_f

391 Right-hand side for the nonlinear system.

392 factor : float

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

394 u0 : dtype_u

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

396 t : float

397 Current time.

399 Returns

400 -------

401 u : dtype_u

402 Solution.

403 """

405 u = self.dtype_u(u0)

406 T = self.M - factor * self.K

407 b = self.dtype_u(rhs)

409 self.bc.apply(T, b.values.vector())

411 df.solve(T, u.values.vector(), b.values.vector())

413 return u

415 def eval_f(self, u, t):

416 """

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

419 Parameters

420 ----------

421 u : dtype_u

422 Current values of the numerical solution.

423 t : float

424 Current time at which the numerical solution is computed.

426 Returns

427 -------

428 f : dtype_f

429 The right-hand side divided into two parts.

430 """

432 f = self.dtype_f(self.V)

434 self.K.mult(u.values.vector(), f.impl.values.vector())

436 self.g.t = t

437 f.expl = self.dtype_u(df.interpolate(self.g, self.V))

438 f.expl = self.apply_mass_matrix(f.expl)

440 return f

442 def fix_residual(self, res):

443 """

444 Applies homogeneous Dirichlet boundary conditions to the residual

446 Parameters

447 ----------

448 res : dtype_u

449 Residual

450 """

451 self.bc_hom.apply(res.values.vector())

452 return None

455# noinspection PyUnusedLocal

456class fenics_heat_mass_timebc(fenics_heat_mass):

457 r"""

458 Example implementing the forced one-dimensional heat equation with time-dependent Dirichlet boundary conditions

460 .. math::

461 \frac{d u}{d t} = \nu \frac{d^2 u}{d x^2} + f

463 for :math:x \in \Omega:=[0,1], where the forcing term :math:f is defined by

465 .. math::

466 f(x, t) = -\cos(\pi x) (\sin(t) - \nu \pi^2 \cos(t)).

468 and the boundary conditions are given by

470 .. math::

471 u(x, t) = \cos(\pi x)\cos(t).

473 The exact solution of the problem is given by

475 .. math::

476 u(x, t) = \cos(\pi x)\cos(t) + c.

478 In this class the problem is implemented in the way that the spatial part is solved using FEniCS [1]_. Hence, the problem

479 is reformulated to the *weak formulation*

481 .. math:

482 \int_\Omega u_t v\,dx = - \nu \int_\Omega \nabla u \nabla v\,dx + \int_\Omega f v\,dx.

484 The forcing term is treated explicitly, and is expressed via the mass matrix resulting from the left-hand side term

485 :math:\int_\Omega u_t v\,dx, and the other part will be treated in an implicit way.

487 Parameters

488 ----------

489 c_nvars : int, optional

490 Spatial resolution, i.e., numbers of degrees of freedom in space.

491 t0 : float, optional

492 Starting time.

493 family : str, optional

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

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

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

497 order : int, optional

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

499 refinements : int, optional

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

501 nu : float, optional

502 Diffusion coefficient :math:\nu.

504 Attributes

505 ----------

506 V : FunctionSpace

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

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

509 Denotes the expression :math:\int_\Omega u_t v\,dx.

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

511 Denotes the expression :math:- \nu \int_\Omega \nabla u \nabla v\,dx.

512 g : Expression

513 The forcing term :math:f in the heat equation.

514 bc : DirichletBC

515 Denotes the time-dependent Dirichlet boundary conditions.

516 bc_hom : DirichletBC

517 Denotes the homogeneous Dirichlet boundary conditions, potentially required for fixing the residual

518 fix_bc_for_residual: boolean

519 flag to indicate that the residual requires special treatment due to boundary conditions

521 References

522 ----------

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

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

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

526 Wells and others. Springer (2012).

527 """

529 def __init__(self, c_nvars=128, t0=0.0, family='CG', order=4, refinements=1, nu=0.1, c=0.0):

530 """Initialization routine"""

532 # define the Dirichlet boundary

533 def Boundary(x, on_boundary):

534 return on_boundary

536 super().__init__(c_nvars, t0, family, order, refinements, nu, c)

538 self.u_D = df.Expression('cos(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t0, degree=self.order)

539 self.bc = df.DirichletBC(self.V, self.u_D, Boundary)

540 self.bc_hom = df.DirichletBC(self.V, df.Constant(0), Boundary)

542 # set forcing term as expression

543 self.g = df.Expression(

544 '-cos(a*x[0]) * (sin(t) - b*a*a*cos(t))',

545 a=np.pi,

546 b=self.nu,

547 t=self.t0,

548 degree=self.order,

549 )

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

552 r"""

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

555 Parameters

556 ----------

557 rhs : dtype_f

558 Right-hand side for the nonlinear system.

559 factor : float

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

561 u0 : dtype_u

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

563 t : float

564 Current time.

566 Returns

567 -------

568 u : dtype_u

569 Solution.

570 """

572 u = self.dtype_u(u0)

573 T = self.M - factor * self.K

574 b = self.dtype_u(rhs)

576 self.u_D.t = t

578 self.bc.apply(T, b.values.vector())

579 self.bc.apply(b.values.vector())

581 df.solve(T, u.values.vector(), b.values.vector())

583 return u

585 def u_exact(self, t):

586 r"""

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

589 Parameters

590 ----------

591 t : float

592 Time of the exact solution.

594 Returns

595 -------

596 me : dtype_u

597 Exact solution.

598 """

600 u0 = df.Expression('cos(a*x[0]) * cos(t) + c', c=self.c, a=np.pi, t=t, degree=self.order)

601 me = self.dtype_u(df.interpolate(u0, self.V), val=self.V)

603 return me