Coverage for pySDC/implementations/problem_classes/HeatEquation_1D_FEniCS_matrix

1import logging

3import dolfin as df

4import numpy as np

6from pySDC.core.Problem import ptype

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

10# noinspection PyUnusedLocal

11class fenics_heat(ptype):

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

100

101 # set mesh and refinement (for multilevel)

102 mesh = df.UnitIntervalMesh(c_nvars)

103 for _ in range(refinements):

104 mesh = df.refine(mesh)

105

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()[:]))

110

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 )

116

117 # Stiffness term (Laplace)

118 u = df.TrialFunction(self.V)

119 v = df.TestFunction(self.V)

120 a_K = -1.0 * df.inner(df.nabla_grad(u), self.nu * df.nabla_grad(v)) * df.dx

121

122 # Mass term

123 a_M = u * v * df.dx

124

125 self.M = df.assemble(a_M)

126 self.K = df.assemble(a_K)

127

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)

131

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 )

140

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}`.

144

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.

155

156 Returns

157 -------

158 u : dtype_u

159 Solution.

160 """

161

162 b = self.apply_mass_matrix(rhs)

163

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

168

169 return u

170

171 def __eval_fexpl(self, u, t):

172 """

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

174

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.

181

182 Returns

183 -------

184 fexpl : dtype_u

185 Explicit part of the right-hand side.

186 """

187

188 self.g.t = t

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

190

191 return fexpl

192

193 def __eval_fimpl(self, u, t):

194 """

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

196

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.

203

204 Returns

205 -------

206 fimpl : dtype_u

207 Explicit part of the right-hand side.

208 """

209

210 tmp = self.dtype_u(self.V)

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

212 fimpl = self.__invert_mass_matrix(tmp)

213

214 return fimpl

215

216 def eval_f(self, u, t):

217 """

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

219

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.

226

227 Returns

228 -------

229 f : dtype_f

230 The right-hand side divided into two parts.

231 """

232

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

237

238 def apply_mass_matrix(self, u):

239 r"""

240 Routine to apply mass matrix.

241

242 Parameters

243 ----------

244 u : dtype_u

245 Current values of the numerical solution.

246

247 Returns

248 -------

249 me : dtype_u

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

251 """

252

253 me = self.dtype_u(self.V)

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

255

256 return me

257

258 def __invert_mass_matrix(self, u):

259 r"""

260 Helper routine to invert mass matrix.

261

262 Parameters

263 ----------

264 u : dtype_u

265 Current values of the numerical solution.

266

267 Returns

268 -------

269 me : dtype_u

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

271 """

272

273 me = self.dtype_u(self.V)

274

275 b = self.dtype_u(u)

276 M = self.M

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

278

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

280 return me

281

282 def u_exact(self, t):

283 r"""

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

285

286 Parameters

287 ----------

288 t : float

289 Time of the exact solution.

290

291 Returns

292 -------

293 me : dtype_u

294 Exact solution.

295 """

296

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)

299

300 return me

301

302

303# noinspection PyUnusedLocal

304class fenics_heat_mass(fenics_heat):

305 r"""

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

307

308 .. math::

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

310

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

312

313 .. math::

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

315

316 For initial conditions with constant c and

317

318 .. math::

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

320

321 the exact solution of the problem is given by

322

323 .. math::

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

325

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*

328

329 .. math:

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

331

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.

334

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`.

351

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

368

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 """

376

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"""

379

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

381

382 self.fix_bc_for_residual = True

383

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}`.

387

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.

398

399 Returns

400 -------

401 u : dtype_u

402 Solution.

403 """

404

405 u = self.dtype_u(u0)

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

407 b = self.dtype_u(rhs)

408

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

410

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

412

413 return u

414

415 def eval_f(self, u, t):

416 """

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

418

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.

425

426 Returns

427 -------

428 f : dtype_f

429 The right-hand side divided into two parts.

430 """

431

432 f = self.dtype_f(self.V)

433

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

435

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)

439

440 return f

441

442 def fix_residual(self, res):

443 """

444 Applies homogeneous Dirichlet boundary conditions to the residual

445

446 Parameters

447 ----------

448 res : dtype_u

449 Residual

450 """

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

452 return None

453

454

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

459

460 .. math::

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

462

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

464

465 .. math::

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

467

468 and the boundary conditions are given by

469

470 .. math::

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

472

473 The exact solution of the problem is given by

474

475 .. math::

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

477

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*

480

481 .. math:

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

483

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.

486

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`.

503

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

520

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 """

528

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"""

531

532 # define the Dirichlet boundary

533 def Boundary(x, on_boundary):

534 return on_boundary

535

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

537

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)

541

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 )

550

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}`.

554

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.

565

566 Returns

567 -------

568 u : dtype_u

569 Solution.

570 """

571

572 u = self.dtype_u(u0)

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

574 b = self.dtype_u(rhs)

575

576 self.u_D.t = t

577

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

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

580

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

582

583 return u

584

585 def u_exact(self, t):

586 r"""

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

588

589 Parameters

590 ----------

591 t : float

592 Time of the exact solution.

593

594 Returns

595 -------

596 me : dtype_u

597 Exact solution.

598 """

599

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)

602

603 return me

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

112 statements