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
On 28 Sep 2016, at 17:11, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
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.
One thing I can imagine might be happening is that the iterative solvers only solve the problem to (by default) a relative tolerance of 1e-8. In your home-baked code do you use some kind of direct method which solves to a tighter tolerance? Lawrence
Dear Lawrence, Yes, I think I solve it directly. I just use numpy.linalg.solve. Could you indicate how to increase the tolerance? Thanks, Tomasz ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 28 September 2016 17:22:40 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] FSI coupling disagreement
On 28 Sep 2016, at 17:11, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
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.
One thing I can imagine might be happening is that the iterative solvers only solve the problem to (by default) a relative tolerance of 1e-8. In your home-baked code do you use some kind of direct method which solves to a tighter tolerance? Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 28 Sep 2016, at 17:38, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Dear Lawrence,
Yes, I think I solve it directly. I just use numpy.linalg.solve. Could you indicate how to increase the tolerance?
Please see http://firedrakeproject.org/solving-interface.html#specifying-solution-metho... I think you want the ksp_rtol solver parameter. Lawrence
Thank you Lawrence, that solved it. Tomasz ________________________________ From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk> Sent: 28 September 2016 18:03:29 To: firedrake@imperial.ac.uk Subject: Re: [firedrake] FSI coupling disagreement
On 28 Sep 2016, at 17:38, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
Dear Lawrence,
Yes, I think I solve it directly. I just use numpy.linalg.solve. Could you indicate how to increase the tolerance?
Please see http://firedrakeproject.org/solving-interface.html#specifying-solution-metho... I think you want the ksp_rtol solver parameter. Lawrence _______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (2)
-
Lawrence Mitchell
-
Tomasz Salwa [RPG]