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

1import logging 

2 

3import dolfin as df 

4import numpy as np 

5 

6from pySDC.core.problem import Problem 

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

8 

9 

10# noinspection PyUnusedLocal 

11class fenics_heat(Problem): 

12 r""" 

13 Forced 1D heat equation with FEniCS and Dirichlet BCs, IMEX, with the mass matrix inverted in the right-hand side. 

14 

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

16 

17 .. math:: 

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

19 

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

21 

22 .. math:: 

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

24 

25 For initial conditions with constant c and 

26 

27 .. math:: 

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

29 

30 the exact solution of the problem is given by 

31 

32 .. math:: 

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

34 

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

36 is reformulated to the *weak formulation* 

37 

38 .. math: 

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

40 

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

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

43 

44 Parameters 

45 ---------- 

46 c_nvars : int, optional 

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

48 t0 : float, optional 

49 Starting time. 

50 family : str, optional 

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

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

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

54 order : int, optional 

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

56 refinements : int, optional 

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

58 nu : float, optional 

59 Diffusion coefficient :math:`\nu`. 

60 c : float, optional 

61 Constant for the Dirichlet boundary condition :math:`c`. 

62 

63 Attributes 

64 ---------- 

65 V : FunctionSpace 

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

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

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

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

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

71 g : Expression 

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

73 bc : DirichletBC 

74 Denotes the Dirichlet boundary conditions. 

75 

76 References 

77 ---------- 

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

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

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

81 Wells and others. Springer (2012). 

82 """ 

83 

84 dtype_u = fenics_mesh 

85 dtype_f = rhs_fenics_mesh 

86 

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

88 """Initialization routine""" 

89 

90 # define the Dirichlet boundary 

91 def Boundary(x, on_boundary): 

92 return on_boundary 

93 

94 # set logger level for FFC and dolfin 

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

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

97 

98 # set solver and form parameters 

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

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

101 df.parameters['allow_extrapolation'] = True 

102 

103 # set mesh and refinement (for multilevel) 

104 mesh = df.UnitIntervalMesh(c_nvars) 

105 for _ in range(refinements): 

106 mesh = df.refine(mesh) 

107 

108 # define function space for future reference 

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

110 tmp = df.Function(self.V) 

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

112 

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

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

115 self._makeAttributeAndRegister( 

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

117 ) 

118 

119 # Stiffness term (Laplace) 

120 u = df.TrialFunction(self.V) 

121 v = df.TestFunction(self.V) 

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

123 

124 # Mass term 

125 a_M = u * v * df.dx 

126 

127 self.M = df.assemble(a_M) 

128 self.K = df.assemble(a_K) 

129 

130 # set boundary values 

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

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

133 

134 # set forcing term as expression 

135 self.g = df.Expression( 

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

137 a=np.pi, 

138 b=self.nu, 

139 t=self.t0, 

140 degree=self.order, 

141 ) 

142 

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

144 r""" 

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

146 

147 Parameters 

148 ---------- 

149 rhs : dtype_f 

150 Right-hand side for the nonlinear system. 

151 factor : float 

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

153 u0 : dtype_u 

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

155 t : float 

156 Current time. 

157 

158 Returns 

159 ------- 

160 u : dtype_u 

161 Solution. 

162 """ 

163 

164 b = self.apply_mass_matrix(rhs) 

165 

166 u = self.dtype_u(u0) 

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

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

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

170 

171 return u 

172 

173 def eval_f_increment(self, base, delta, t): 

174 """ 

175 Evaluate the right-hand side increment. 

176 

177 Parameters 

178 ---------- 

179 base : dtype_u 

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

181 delta : dtype_u 

182 The correction. 

183 t : float 

184 Physical time, accepted for interface compatibility. 

185 

186 Returns 

187 ------- 

188 dtype_f 

189 The increment, with a zero explicit part. 

190 """ 

191 increment = self.eval_f(delta, t) 

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

193 return increment 

194 

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

196 r""" 

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

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

199 

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

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

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

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

204 

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

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

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

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

209 

210 Parameters 

211 ---------- 

212 r : dtype_u 

213 Right-hand side of the correction equation. 

214 factor : float 

215 Implicit prefactor assembled by the sweeper. 

216 base : dtype_u 

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

218 f_base : dtype_f 

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

220 t : float 

221 Physical time, accepted for interface compatibility. 

222 

223 Returns 

224 ------- 

225 dtype_u 

226 The correction. 

227 """ 

228 b = self.apply_mass_matrix(r) 

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

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

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

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

233 return delta 

234 

235 def __eval_fexpl(self, u, t): 

236 """ 

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

238 

239 Parameters 

240 ---------- 

241 u : dtype_u 

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

243 t : float 

244 Current time at which the numerical solution is computed. 

245 

246 Returns 

247 ------- 

248 fexpl : dtype_u 

249 Explicit part of the right-hand side. 

250 """ 

251 

252 self.g.t = t 

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

254 

255 return fexpl 

256 

257 def __eval_fimpl(self, u, t): 

258 """ 

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

260 

261 Parameters 

262 ---------- 

263 u : dtype_u 

264 Current values of the numerical solution. 

265 t : float 

266 Current time at which the numerical solution is computed. 

267 

268 Returns 

269 ------- 

270 fimpl : dtype_u 

271 Explicit part of the right-hand side. 

272 """ 

273 

274 tmp = self.dtype_u(self.V) 

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

276 fimpl = self.__invert_mass_matrix(tmp) 

277 

278 return fimpl 

279 

280 def eval_f(self, u, t): 

281 """ 

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

283 

284 Parameters 

285 ---------- 

286 u : dtype_u 

287 Current values of the numerical solution. 

288 t : float 

289 Current time at which the numerical solution is computed. 

290 

291 Returns 

292 ------- 

293 f : dtype_f 

294 The right-hand side divided into two parts. 

295 """ 

296 

297 f = self.dtype_f(self.V) 

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

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

300 return f 

301 

302 def apply_mass_matrix(self, u): 

303 r""" 

304 Routine to apply mass matrix. 

305 

306 Parameters 

307 ---------- 

308 u : dtype_u 

309 Current values of the numerical solution. 

310 

311 Returns 

312 ------- 

313 me : dtype_u 

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

315 """ 

316 

317 me = self.dtype_u(self.V) 

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

319 

320 return me 

321 

322 def __invert_mass_matrix(self, u): 

323 r""" 

324 Helper routine to invert mass matrix. 

325 

326 Parameters 

327 ---------- 

328 u : dtype_u 

329 Current values of the numerical solution. 

330 

331 Returns 

332 ------- 

333 me : dtype_u 

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

335 """ 

336 

337 me = self.dtype_u(self.V) 

338 

339 b = self.dtype_u(u) 

340 M = self.M 

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

342 

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

344 return me 

345 

346 def u_exact(self, t): 

347 r""" 

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

349 

350 Parameters 

351 ---------- 

352 t : float 

353 Time of the exact solution. 

354 

355 Returns 

356 ------- 

357 me : dtype_u 

358 Exact solution. 

359 """ 

360 

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

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

363 

364 return me 

365 

366 

367# noinspection PyUnusedLocal 

368class fenics_heat_mass(fenics_heat): 

369 r""" 

370 Forced 1D heat equation with FEniCS and Dirichlet BCs, IMEX, with the mass matrix applied instead of inverted. 

371 

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

373 

374 .. math:: 

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

376 

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

378 

379 .. math:: 

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

381 

382 For initial conditions with constant c and 

383 

384 .. math:: 

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

386 

387 the exact solution of the problem is given by 

388 

389 .. math:: 

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

391 

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

393 is reformulated to the *weak formulation* 

394 

395 .. math: 

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

397 

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

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

400 

401 Parameters 

402 ---------- 

403 c_nvars : int, optional 

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

405 t0 : float, optional 

406 Starting time. 

407 family : str, optional 

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

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

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

411 order : int, optional 

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

413 refinements : int, optional 

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

415 nu : float, optional 

416 Diffusion coefficient :math:`\nu`. 

417 c : float, optional 

418 Constant for the Dirichlet boundary condition :math:`c`. 

419 

420 Attributes 

421 ---------- 

422 V : FunctionSpace 

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

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

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

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

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

428 g : Expression 

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

430 bc : DirichletBC 

431 Denotes the Dirichlet boundary conditions. 

432 bc_hom : DirichletBC 

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

434 fix_bc_for_residual: boolean 

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

436 

437 References 

438 ---------- 

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

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

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

442 Wells and others. Springer (2012). 

443 """ 

444 

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

446 """Initialization routine""" 

447 

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

449 

450 self.fix_bc_for_residual = True 

451 

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

453 r""" 

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

455 

456 Parameters 

457 ---------- 

458 rhs : dtype_f 

459 Right-hand side for the nonlinear system. 

460 factor : float 

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

462 u0 : dtype_u 

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

464 t : float 

465 Current time. 

466 

467 Returns 

468 ------- 

469 u : dtype_u 

470 Solution. 

471 """ 

472 

473 u = self.dtype_u(u0) 

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

475 b = self.dtype_u(rhs) 

476 

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

478 

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

480 

481 return u 

482 

483 def eval_f(self, u, t): 

484 """ 

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

486 

487 Parameters 

488 ---------- 

489 u : dtype_u 

490 Current values of the numerical solution. 

491 t : float 

492 Current time at which the numerical solution is computed. 

493 

494 Returns 

495 ------- 

496 f : dtype_f 

497 The right-hand side divided into two parts. 

498 """ 

499 

500 f = self.dtype_f(self.V) 

501 

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

503 

504 self.g.t = t 

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

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

507 

508 return f 

509 

510 def fix_residual(self, res): 

511 """ 

512 Applies homogeneous Dirichlet boundary conditions to the residual 

513 

514 Parameters 

515 ---------- 

516 res : dtype_u 

517 Residual 

518 """ 

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

520 return None 

521 

522 

523# noinspection PyUnusedLocal 

524class fenics_heat_mass_timebc(fenics_heat_mass): 

525 r""" 

526 Forced 1D heat equation with FEniCS, IMEX with the mass matrix applied, and time-dependent Dirichlet BCs. 

527 

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

529 

530 .. math:: 

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

532 

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

534 

535 .. math:: 

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

537 

538 and the boundary conditions are given by 

539 

540 .. math:: 

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

542 

543 The exact solution of the problem is given by 

544 

545 .. math:: 

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

547 

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

549 is reformulated to the *weak formulation* 

550 

551 .. math: 

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

553 

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

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

556 

557 Parameters 

558 ---------- 

559 c_nvars : int, optional 

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

561 t0 : float, optional 

562 Starting time. 

563 family : str, optional 

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

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

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

567 order : int, optional 

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

569 refinements : int, optional 

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

571 nu : float, optional 

572 Diffusion coefficient :math:`\nu`. 

573 c : float, optional 

574 Constant :math:`c` added to the exact solution and hence to the time-dependent Dirichlet boundary condition. 

575 

576 Attributes 

577 ---------- 

578 V : FunctionSpace 

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

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

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

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

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

584 g : Expression 

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

586 bc : DirichletBC 

587 Denotes the time-dependent Dirichlet boundary conditions. 

588 bc_hom : DirichletBC 

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

590 fix_bc_for_residual: boolean 

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

592 

593 References 

594 ---------- 

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

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

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

598 Wells and others. Springer (2012). 

599 """ 

600 

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

602 """Initialization routine""" 

603 

604 # define the Dirichlet boundary 

605 def Boundary(x, on_boundary): 

606 return on_boundary 

607 

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

609 

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

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

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

613 

614 # set forcing term as expression 

615 self.g = df.Expression( 

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

617 a=np.pi, 

618 b=self.nu, 

619 t=self.t0, 

620 degree=self.order, 

621 ) 

622 

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

624 r""" 

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

626 

627 Parameters 

628 ---------- 

629 rhs : dtype_f 

630 Right-hand side for the nonlinear system. 

631 factor : float 

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

633 u0 : dtype_u 

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

635 t : float 

636 Current time. 

637 

638 Returns 

639 ------- 

640 u : dtype_u 

641 Solution. 

642 """ 

643 

644 u = self.dtype_u(u0) 

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

646 b = self.dtype_u(rhs) 

647 

648 self.u_D.t = t 

649 

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

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

652 

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

654 

655 return u 

656 

657 def u_exact(self, t): 

658 r""" 

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

660 

661 Parameters 

662 ---------- 

663 t : float 

664 Time of the exact solution. 

665 

666 Returns 

667 ------- 

668 me : dtype_u 

669 Exact solution. 

670 """ 

671 

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

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

674 

675 return me