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

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_f_increment(self, base, delta, t):

172 """

173 Evaluate the right-hand side increment.

174

175 Parameters

176 ----------

177 base : dtype_u

178 The base state, unused: the implicit part is linear.

179 delta : dtype_u

180 The correction.

181 t : float

182 Physical time, accepted for interface compatibility.

183

184 Returns

185 -------

186 dtype_f

187 The increment, with a zero explicit part.

188 """

189 increment = self.eval_f(delta, t)

190 increment.expl = self.dtype_u(self.V, val=0.0)

191 return increment

192

193 def solve_system_delta(self, r, factor, base, f_base, t):

194 r"""

195 Solve :math:`\delta - factor\,[f(w+\delta) - f(w)] = r`, i.e.

196 :math:`(M - factor\,K)\,\delta = M r` with **zero** boundary data.

197

198 This is the piece ``linear_implicit=True`` cannot supply on this backend. That shortcut

199 reuses the stock ``solve_system``, which applies *inhomogeneous* Dirichlet data to whatever

200 right-hand side it is handed, and a correction must carry zero boundary data. Applying

201 ``bc_hom`` instead is the whole difference.

202

203 Without this the sweeper falls back to the substitution :math:`y = w + \delta`, which is

204 exact but reads the level's :math:`\mathcal{O}(1)` state. That is merely no benefit while

205 the level is at backend precision, and is a wrong answer once it is not -- measured here as

206 1.4e-05 with the coarse level at ``float32``.

207

208 Parameters

209 ----------

210 r : dtype_u

211 Right-hand side of the correction equation.

212 factor : float

213 Implicit prefactor assembled by the sweeper.

214 base : dtype_u

215 Base state, unused: the implicit operator is linear.

216 f_base : dtype_f

217 ``f`` evaluated at ``base``, unused for the same reason.

218 t : float

219 Physical time, accepted for interface compatibility.

220

221 Returns

222 -------

223 dtype_u

224 The correction.

225 """

226 b = self.apply_mass_matrix(r)

227 delta = self.dtype_u(self.V, val=0.0)

228 system = self.M - factor * self.K

229 self.bc_hom.apply(system, b.values.vector())

230 df.solve(system, delta.values.vector(), b.values.vector())

231 return delta

232

233 def __eval_fexpl(self, u, t):

234 """

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

236

237 Parameters

238 ----------

239 u : dtype_u

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

241 t : float

242 Current time at which the numerical solution is computed.

243

244 Returns

245 -------

246 fexpl : dtype_u

247 Explicit part of the right-hand side.

248 """

249

250 self.g.t = t

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

252

253 return fexpl

254

255 def __eval_fimpl(self, u, t):

256 """

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

258

259 Parameters

260 ----------

261 u : dtype_u

262 Current values of the numerical solution.

263 t : float

264 Current time at which the numerical solution is computed.

265

266 Returns

267 -------

268 fimpl : dtype_u

269 Explicit part of the right-hand side.

270 """

271

272 tmp = self.dtype_u(self.V)

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

274 fimpl = self.__invert_mass_matrix(tmp)

275

276 return fimpl

277

278 def eval_f(self, u, t):

279 """

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

281

282 Parameters

283 ----------

284 u : dtype_u

285 Current values of the numerical solution.

286 t : float

287 Current time at which the numerical solution is computed.

288

289 Returns

290 -------

291 f : dtype_f

292 The right-hand side divided into two parts.

293 """

294

295 f = self.dtype_f(self.V)

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

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

298 return f

299

300 def apply_mass_matrix(self, u):

301 r"""

302 Routine to apply mass matrix.

303

304 Parameters

305 ----------

306 u : dtype_u

307 Current values of the numerical solution.

308

309 Returns

310 -------

311 me : dtype_u

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

313 """

314

315 me = self.dtype_u(self.V)

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

317

318 return me

319

320 def __invert_mass_matrix(self, u):

321 r"""

322 Helper routine to invert mass matrix.

323

324 Parameters

325 ----------

326 u : dtype_u

327 Current values of the numerical solution.

328

329 Returns

330 -------

331 me : dtype_u

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

333 """

334

335 me = self.dtype_u(self.V)

336

337 b = self.dtype_u(u)

338 M = self.M

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

340

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

342 return me

343

344 def u_exact(self, t):

345 r"""

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

347

348 Parameters

349 ----------

350 t : float

351 Time of the exact solution.

352

353 Returns

354 -------

355 me : dtype_u

356 Exact solution.

357 """

358

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

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

361

362 return me

363

364

365# noinspection PyUnusedLocal

366class fenics_heat_mass(fenics_heat):

367 r"""

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

369

370 .. math::

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

372

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

374

375 .. math::

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

377

378 For initial conditions with constant c and

379

380 .. math::

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

382

383 the exact solution of the problem is given by

384

385 .. math::

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

387

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

389 is reformulated to the *weak formulation*

390

391 .. math:

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

393

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

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

396

397 Parameters

398 ----------

399 c_nvars : int, optional

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

401 t0 : float, optional

402 Starting time.

403 family : str, optional

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

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

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

407 order : int, optional

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

409 refinements : int, optional

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

411 nu : float, optional

412 Diffusion coefficient :math:`\nu`.

413

414 Attributes

415 ----------

416 V : FunctionSpace

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

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

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

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

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

422 g : Expression

423 The forcing term :math:`f` in the heat equation.

424 bc : DirichletBC

425 Denotes the Dirichlet boundary conditions.

426 bc_hom : DirichletBC

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

428 fix_bc_for_residual: boolean

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

430

431 References

432 ----------

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

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

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

436 Wells and others. Springer (2012).

437 """

438

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

440 """Initialization routine"""

441

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

443

444 self.fix_bc_for_residual = True

445

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

447 r"""

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

449

450 Parameters

451 ----------

452 rhs : dtype_f

453 Right-hand side for the nonlinear system.

454 factor : float

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

456 u0 : dtype_u

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

458 t : float

459 Current time.

460

461 Returns

462 -------

463 u : dtype_u

464 Solution.

465 """

466

467 u = self.dtype_u(u0)

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

469 b = self.dtype_u(rhs)

470

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

472

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

474

475 return u

476

477 def eval_f(self, u, t):

478 """

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

480

481 Parameters

482 ----------

483 u : dtype_u

484 Current values of the numerical solution.

485 t : float

486 Current time at which the numerical solution is computed.

487

488 Returns

489 -------

490 f : dtype_f

491 The right-hand side divided into two parts.

492 """

493

494 f = self.dtype_f(self.V)

495

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

497

498 self.g.t = t

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

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

501

502 return f

503

504 def fix_residual(self, res):

505 """

506 Applies homogeneous Dirichlet boundary conditions to the residual

507

508 Parameters

509 ----------

510 res : dtype_u

511 Residual

512 """

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

514 return None

515

516

517# noinspection PyUnusedLocal

518class fenics_heat_mass_timebc(fenics_heat_mass):

519 r"""

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

521

522 .. math::

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

524

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

526

527 .. math::

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

529

530 and the boundary conditions are given by

531

532 .. math::

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

534

535 The exact solution of the problem is given by

536

537 .. math::

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

539

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

541 is reformulated to the *weak formulation*

542

543 .. math:

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

545

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

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

548

549 Parameters

550 ----------

551 c_nvars : int, optional

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

553 t0 : float, optional

554 Starting time.

555 family : str, optional

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

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

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

559 order : int, optional

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

561 refinements : int, optional

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

563 nu : float, optional

564 Diffusion coefficient :math:`\nu`.

565

566 Attributes

567 ----------

568 V : FunctionSpace

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

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

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

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

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

574 g : Expression

575 The forcing term :math:`f` in the heat equation.

576 bc : DirichletBC

577 Denotes the time-dependent Dirichlet boundary conditions.

578 bc_hom : DirichletBC

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

580 fix_bc_for_residual: boolean

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

582

583 References

584 ----------

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

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

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

588 Wells and others. Springer (2012).

589 """

590

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

592 """Initialization routine"""

593

594 # define the Dirichlet boundary

595 def Boundary(x, on_boundary):

596 return on_boundary

597

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

599

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

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

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

603

604 # set forcing term as expression

605 self.g = df.Expression(

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

607 a=np.pi,

608 b=self.nu,

609 t=self.t0,

610 degree=self.order,

611 )

612

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

614 r"""

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

616

617 Parameters

618 ----------

619 rhs : dtype_f

620 Right-hand side for the nonlinear system.

621 factor : float

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

623 u0 : dtype_u

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

625 t : float

626 Current time.

627

628 Returns

629 -------

630 u : dtype_u

631 Solution.

632 """

633

634 u = self.dtype_u(u0)

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

636 b = self.dtype_u(rhs)

637

638 self.u_D.t = t

639

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

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

642

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

644

645 return u

646

647 def u_exact(self, t):

648 r"""

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

650

651 Parameters

652 ----------

653 t : float

654 Time of the exact solution.

655

656 Returns

657 -------

658 me : dtype_u

659 Exact solution.

660 """

661

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

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

664

665 return me

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

123 statements