Hi Tomasz,

This is a subtlety of UFL. *All* arguments and coefficients in an interior facet integral have to be "restricted". This even applies to ones which are continuous and for which the "+" and "-" restrictions are therefore the same. You can just arbitrarily decide to restrict the continuous things to the plus side, or take avg() of the lack of symmetry offends your sense of aesthetics.

Regards,

David

On Tue, 26 Apr 2016 at 15:46 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:

Hi,

Previous problem continued.

I have a working solver for separate subsystems in mixed function space now. I'd like to introduce the coupling terms at the interface, but they cause a UFL exception "Discontinuous type Argument must be restricted."

Full code is here: https://bitbucket.org/mmtjs/lin_coupled_fd_2d

The relevant part from lib/coupled_system.py is:

...

# define the normal vector at the interface

        self.n = fd.FacetNormal(self.mesh)
        self.n_int = (self.B.I("+")*self.W.I("-") - self.B.I("-")*self.W.I("+"))*self.n("+")  # btw how can write this to file to verify?

...

# the second line causes exception (and later similar terms too)

        a_U = fd.dot( self.trial_s * self.B.I, self.v_s * self.B.I) * fd.dx \
                    + fd.dot( self.v_s, self.n_int ) * self.trial_f * fd.dS

Regards,

Tomasz

---

From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk>
Sent: 21 April 2016 11:04

To: firedrake
Subject: Re: [firedrake] Laplace equation on subdomain

 

Cool, works now, dynamic case too. I'll do the same for the structural part.

Cheers,

Tomasz

---

From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of David Ham <David.Ham@imperial.ac.uk>

Sent: 21 April 2016 09:38

To: firedrake
Subject: Re: [firedrake] Laplace equation on subdomain

Hi Tomasz,

See the new pull request. In essence, the projection is not safe in the indicator function case.

Regards,

David

On Wed, 20 Apr 2016 at 15:53 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:

Hi David,

Thank you, now it converges properly.

However, I noticed that the error depends on the mesh resolution beyond the fluid subdomain, which I think should not be the case (at least when I refine it in x direction). Perhaps the solution is ok, but the source is in the way I impose the analytical one, or compute the error. Can you diagnose this?

Bitbucket's repo is updated.

Regards,

Tomasz

---

From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of David Ham <David.Ham@imperial.ac.uk>

Sent: 20 April 2016 10:05

To: firedrake
Subject: Re: [firedrake] Laplace equation on subdomain

Hi Tomasz,

I have just put in a pull request fixing your problem. In essence, you were only applying half of the indicator function fix. In order to correctly use indicator functions to solve problems over a subdomain, it is necessary to *both* use the custom boundary condition to force the solution in the rest of the domain to zero *and* multiply the integrands in your forms by the indicator function. You were only doing the former and not the latter. Making this fix causes the error to drop to 10^-5. I haven't done a convergence test to check that the error converges at the right rate, but you should really do so.

Regards,

David

On Tue, 19 Apr 2016 at 11:52 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:

Hi David,

Here's the code:

https://bitbucket.org/mmtjs/laplace_extended_domain

It's roughly the same as before, just with a gmsh-generated mesh.

Regards,

Tomasz

---

From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of David Ham <David.Ham@imperial.ac.uk>
Sent: 19 April 2016 09:54
To: firedrake

Subject: Re: [firedrake] Laplace equation on subdomain

Hi Tomasz,

Can you please post the code you are using for this (preferably by pointing us to a git repo we can clone). I can't diagnose what is going wrong just by looking at the visualisation, I'm afraid.

Regards,

David

On Mon, 18 Apr 2016 at 12:12 Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:

Hi,

I plotted the difference between analytical and numerical solutions (see attached). It doesn't converge to the right solution with mesh refinement. Is it possible that the Neumann BC is not applied at the interface during the computation, although the extended Dirichlet BCs enforce resulting zero solution beyond the fluid subdomain?

Cheers,
Tomasz

________________________________________
From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk>
Sent: 15 April 2016 16:00
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] Laplace equation on subdomain

> On 15 Apr 2016, at 15:54, Tomasz Salwa [RPG] <mmtjs@leeds.ac.uk> wrote:
>
> Hi David,
>
> It works partially, i.e. the nodes beyond the fluid subdomain are "turned off" now, but it looks like the interface is "smeared". I suspect that the assembled matrix's entries involve also integration over elements one layer beyond the subdomain boundary. Is it possible to exclude them?

If you have neumann conditions on your subdomain boundary you've only constrained the gradient of the solution.  Therefore on the boundary nodes you may have non-zero solution values.

Because the domain you visualise in paraview extends outside the subdomain you'll have (in 1D) something like:


........
        \............
       ^ boundary node

Note how the cell at the boundary has a non-zero coefficient value are one end, but a zero value at the other (assuming piecewise liners). Paraview will linearly interpolate this basis function when visualising, which may be the artefact you're seeing.

Can you construct an analytic solution and determine that you converge at the correct rate as you refine the mesh?  This is a safer test than eyeballing the solution.

Lawrence