Coverage for pySDC/projects/StroemungsRaum/problem_classes/newton_step.py: 100%

1import dolfin as df 

2 

3 

4class NewtonStep(df.NonlinearProblem): 

5 r""" 

6 Newton problem for a single SDC node-to-node step :math:`M w + \Delta t_{QI} N(w) = rhs`. 

7 

8 The right-hand side handed over by the sweeper is an assembled vector, and it is subtracted 

9 from the residual here, at the algebraic level. The alternative -- writing it into the 

10 variational form as :math:`\int_\Omega rhs \cdot v\,dx` so that dolfin's high level 

11 ``solve(F == 0, ...)`` interface can be used -- applies the mass matrix to it, which then 

12 has to be undone by a mass matrix solve beforehand. That round trip is exact in exact 

13 arithmetic, so it buys nothing, while costing a solve per node per sweep and capping the 

14 attainable accuracy at the tolerance of that solve. 

15 

16 ``rhs`` and ``bcs`` are set per solve rather than at construction, so that the form and its 

17 Jacobian can be compiled once and reused with the step size carried by a ``Constant``. 

18 

19 Parameters 

20 ---------- 

21 F : UFL form 

22 Residual form of the step, *without* the right-hand side term. 

23 J : UFL form 

24 Jacobian of ``F``. 

25 

26 Attributes 

27 ---------- 

28 rhs : GenericVector 

29 Right-hand side vector for the current solve, subtracted from the residual. 

30 bcs : list of DirichletBC 

31 Boundary conditions for the current solve, applied in residual form. 

32 """ 

33 

34 def __init__(self, F, J): 

35 super().__init__() 

36 self.F_form = F 

37 self.J_form = J 

38 self.rhs = None 

39 self.bcs = [] 

40 

41 def F(self, b, x): 

42 df.assemble(self.F_form, tensor=b) 

43 b.axpy(-1.0, self.rhs) 

44 for bc in self.bcs: 

45 bc.apply(b, x) 

46 

47 def J(self, A, x): 

48 df.assemble(self.J_form, tensor=A) 

49 for bc in self.bcs: 

50 bc.apply(A)