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