Re: [firedrake] nonlinear diffusion with space-time Robin condition
Hi Onno, As Stephan correctly (as usual) points out, your primary speed issue is that you are using solve instead of LinearVariationalProblem and LinearVariationalSolver. This means you unnecessarily duplicate a whole load of setup every timestep. In this case, the simplest thing to do is to delete your special case for theta = 0 and just have the theta = 0 case run through the nonlinear code path. On my machine this gets you two time units in about two seconds. If you really want to, you can set the solver options in the theta = 0 case to switch off the outer nonlinear solve. However my observation is that this makes negligible additional difference because PETSc only does one solver iteration anyway (it does an extra residual evaluation, but this seems not to be a significant cost). One would not expect Firedrake to be significantly faster than Matlab in this case because the problem is so tiny. The reason that you are observing failure in the theta > 0 case is this term: sqrt(max_value(2.0*h/3.0,0.0)) When this gets differentiated, you get a 1/sqrt(h) term, which has a pole at h=0 (which is not a huge surprise since the expression is actually not differentiable at h=0). This inserts NaNs into the matrix and the solver dies. Firedrake is doing the right thing, it’s just that your problem is not well defined there. You have two options: 1. Use an initial guess which is not zero, and hope that the newton solve never actually hits h=0 at the boundary (which it probably won’t). 2. Change the form to something differentiable (or at least with finite formal derivative) at h=0. Regards, David PS. In many circumstances, solver options are also a significant speed factor. However in this case (a) the tiny problem largely makes this irrelevant and (b) because you have a 1D problem, the default ilu(0) preconditioner acts like a direct solver, which is almost certainly the best option for tiny problems like this. From: <firedrake-bounces@imperial.ac.uk> on behalf of Onno Bokhove <O.Bokhove@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Sunday, 13 May 2018 at 13:15 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] nonlinear diffusion with space-time Robin condition Dear Firedrakers, (i) Will & I tried to build a theta-time-step solver including Crank-Nicolson for theta=1/2 but only get it to work for theta=0 (Forward Euler); see the attached code at "HERE2" which matches (32) in the attached pdf. [This is really a test set-up atm for a larger question on dynamic boundary conditions for an elliptic problem, but the system may be amusing by itself as a tentative FD-example.) (ii) The working explicit Forward Euler code for theta=0 is very slow, at least 10x slower than the same algorithm in a FEM matlab code (-). See "HERE1" in the attached code; why is the code so slow? [I have also attached the basic matlab code.] This example may be amusing because of the space-time Robin condition (14) applied to (11), leading to the Firedrake weak form (32) --when time is discretised as well. The question can also be solved, slightly differently wrt the implementation of the Robin boundary condition, using finite differences, which makes for an amusing comparison. I am not even sure whether this space-time bc and nonlinear diffusion eqn is well-posed but the finite difference and FEM numerics seems to be consistent. Thank you and best wishes, Onno
HI David and Stephan, Thanks. Per the below: (i) The expression is a^1.5 but was written as a*sqrt(a); now that is corrected the nonlinear solver seems to work fine; apparently that rewrite matters for the Newton iteration. (ii) The explicit code is now faster (still a factor 2 slower than matlab) but the implicit code works too for larger time step and different order polynomials. So even in 1D FD is helpful. Cheers, Onno ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Ham, David A <david.ham@imperial.ac.uk> Sent: Sunday, May 13, 2018 4:38 PM To: firedrake Subject: Re: [firedrake] nonlinear diffusion with space-time Robin condition Hi Onno, As Stephan correctly (as usual) points out, your primary speed issue is that you are using solve instead of LinearVariationalProblem and LinearVariationalSolver. This means you unnecessarily duplicate a whole load of setup every timestep. In this case, the simplest thing to do is to delete your special case for theta = 0 and just have the theta = 0 case run through the nonlinear code path. On my machine this gets you two time units in about two seconds. If you really want to, you can set the solver options in the theta = 0 case to switch off the outer nonlinear solve. However my observation is that this makes negligible additional difference because PETSc only does one solver iteration anyway (it does an extra residual evaluation, but this seems not to be a significant cost). One would not expect Firedrake to be significantly faster than Matlab in this case because the problem is so tiny. The reason that you are observing failure in the theta > 0 case is this term: sqrt(max_value(2.0*h/3.0,0.0)) When this gets differentiated, you get a 1/sqrt(h) term, which has a pole at h=0 (which is not a huge surprise since the expression is actually not differentiable at h=0). This inserts NaNs into the matrix and the solver dies. Firedrake is doing the right thing, it’s just that your problem is not well defined there. You have two options: 1. Use an initial guess which is not zero, and hope that the newton solve never actually hits h=0 at the boundary (which it probably won’t). 2. Change the form to something differentiable (or at least with finite formal derivative) at h=0. Regards, David PS. In many circumstances, solver options are also a significant speed factor. However in this case (a) the tiny problem largely makes this irrelevant and (b) because you have a 1D problem, the default ilu(0) preconditioner acts like a direct solver, which is almost certainly the best option for tiny problems like this. From: <firedrake-bounces@imperial.ac.uk> on behalf of Onno Bokhove <O.Bokhove@leeds.ac.uk> Reply-To: firedrake <firedrake@imperial.ac.uk> Date: Sunday, 13 May 2018 at 13:15 To: firedrake <firedrake@imperial.ac.uk> Subject: [firedrake] nonlinear diffusion with space-time Robin condition Dear Firedrakers, (i) Will & I tried to build a theta-time-step solver including Crank-Nicolson for theta=1/2 but only get it to work for theta=0 (Forward Euler); see the attached code at "HERE2" which matches (32) in the attached pdf. [This is really a test set-up atm for a larger question on dynamic boundary conditions for an elliptic problem, but the system may be amusing by itself as a tentative FD-example.) (ii) The working explicit Forward Euler code for theta=0 is very slow, at least 10x slower than the same algorithm in a FEM matlab code (-). See "HERE1" in the attached code; why is the code so slow? [I have also attached the basic matlab code.] This example may be amusing because of the space-time Robin condition (14) applied to (11), leading to the Firedrake weak form (32) --when time is discretised as well. The question can also be solved, slightly differently wrt the implementation of the Robin boundary condition, using finite differences, which makes for an amusing comparison. I am not even sure whether this space-time bc and nonlinear diffusion eqn is well-posed but the finite difference and FEM numerics seems to be consistent. Thank you and best wishes, Onno
On 15 May 2018, at 11:49, Onno Bokhove <O.Bokhove@leeds.ac.uk> wrote:
HI David and Stephan,
Thanks.
Per the below: (i) The expression is a^1.5 but was written as a*sqrt(a); now that is corrected the nonlinear solver seems to work fine; apparently that rewrite matters for the Newton iteration.
The Jacobian is computed by automatic differentiation, so although mathematically a^1.5 and a*sqrt(a) are the same, they are differentiated differently to: 1.5 a^0.5 or 0.5 a / sqrt(a) + sqrt(a) The latter has a numerical singularity at 0. In floating point, if a is zero, you compute 0.0 / 0.0 => NaN Lawrence
Yes realised that now of course! ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: Tuesday, May 15, 2018 12:01 PM To: firedrake@imperial.ac.uk Subject: Re: [firedrake] nonlinear diffusion with space-time Robin condition
On 15 May 2018, at 11:49, Onno Bokhove <O.Bokhove@leeds.ac.uk> wrote:
HI David and Stephan,
Thanks.
Per the below: (i) The expression is a^1.5 but was written as a*sqrt(a); now that is corrected the nonlinear solver seems to work fine; apparently that rewrite matters for the Newton iteration.
The Jacobian is computed by automatic differentiation, so although mathematically a^1.5 and a*sqrt(a) are the same, they are differentiated differently to: 1.5 a^0.5 or 0.5 a / sqrt(a) + sqrt(a) The latter has a numerical singularity at 0. In floating point, if a is zero, you compute 0.0 / 0.0 => NaN Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (3)
-
Ham, David A
-
Lawrence Mitchell
-
Onno Bokhove