Dear all (especially David and Lawrence), I work with a coupled fluid-structure interaction problem. I have a plain python code for verification. For a simple test setup (see attached file init_cond_py.pdf), I expect the same result from Firedrake as from my code. The result is very close, but slightly different. Grey plots show results from my code, coloured from Firedrake. The energy looks basically the same. However, in the inspection of convergence of energy with dt (expected second order convergence), Firedrake code results in a slight energy drift, which is unexpected. The uncoupled subsystems behave properly, so I suspect that the source of this discrepancy lies in the coupling terms. The relevant part is ... # given the definitions: V_f = FunctionSpace(mesh, "CG", 1) V_s = VectorFunctionSpace(mesh, "CG", 1) mixed_V = V_f * V_s trial_f, trial_s = TrialFunctions(mixed_V) v_f, v_s = TestFunctions(mixed_V) phi = Function(V_f) U = Function(V_s) n_normal = FacetNormal(mesh) # n - normal vector at the interface #I_s, I_f - indicator functions, 1(0) in solid, 0(1) in fluid, respectively n = (I_s("+") * I_f("-") - I_s("-") * I_f("+")) * n_normal("-") .... # the coupling terms are: ... + dot( avg(v_s), n ) * avg(trial_f) * dS ... + dot( avg(v_s), n ) * avg(phi) * dS ... - dot( n, avg(trial_s) ) * avg(v_f) * dS ... - dt * dot( n, avg(U) ) * avg(v_f) * dS Can it be that the energy drift is caused by the averages above or by a method Firedrake assembles those terms? I would welcome any suggestions, as I have run out of ideas what I may be doing wrong. The full code is here: https://bitbucket.org/mmtjs/lin_coupled_fd_2d with relevant part in lib/coupled_system.py Thank you, Tomasz