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:

Best wishes,

Fryderyk