On 02/05/16 09:44, Fryderyk Wilczynski wrote:
Hello all,
I'm solving a system of pdes that are linear in (n+1) terms (in time disretised system).
When I write the equations in Firedrake as a bilinear form and solve using:
...
solve(a == L, out_u, bcs=[bc0, bc1, bc2, bc3], solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'})
I get a stable (energy dissipative) scheme, everything seems to be working fine.
If I reformulate the equations so that I am solving the residual form:
...
solve(F == 0, u, bcs=[bc0, bc1, bc2, bc3],solver_parameters={'pc_type': 'lu', 'pc_factor_mat_solver_package': 'mumps'} )
all other things remainimg equal; the scheme breaks down after a number of timesteps, and throws an exception:
If your formulation is actually linear, then these two (especially with an LU solver) should be equivalent. The fact that they're not suggests that maybe your problem as formulated has some nonlinearity that it shouldn't (I think your timestepping scheme results in a linear formulation) and as the system evolves this gets larger. Although eyeballing the formulation it doesn't look like this is the case. That they are not indicates that, possibly, they are not. Can you try some of the advice on debugging convergence failures we have here: http://firedrakeproject.org/solving-interface.html#debugging-convergence-fai... Moreover, something funky is going on. I don't know when your residual formulation breaks down, but I ran it for a while here (up to t=0.5) and had no problems. Additionally, if I refactor your bilinear form code such that I use a LinearVariationalProblem/LinearVariationalSolver pair. The only difference in this case being that you don't rebuild the operator and preconditioner on every solve, then the bilinear formulation blows up at t=0.245. So not sure! Cheers, Lawrence