Hi Eike,

In general, we don't actually have proper functionality for that.

However, assuming you're using some H(div) space for u, including an HDiv(U2xV0) + HDiv(U1xV1) sort of thing, you can obtain the BC u.n=0 by doing, for example,

bcs = [DirichletBC(W2, 0.0, "bottom"), DirichletBC(W2, 0.0, "top")]

Why does this work? Well, this causes any degrees of freedom on the 'bottom' and 'top' boundaries to be set to 0. But the degrees of freedom on the boundary in an H(div) space are precisely u.n (technically with some sort of scaling, but scaled 0 is still 0). And any interior DOFs, or DOFs on other facets, have vanishing normal component on the original facet. So indeed, this does the right thing.

If your velocity was in something like VectorFunctionSpace(mesh, "CG", n), unfortunately this approach wouldn't work, and I think you'd have to apply the BC weakly by adding on a surface integral. But in the H(div) case, what I described above happens to work.

Hope this helps!

Best,

Andrew

On 27 December 2014 at 14:14, Eike Mueller <e.mueller@bath.ac.uk> wrote:

Dear firedrakers,

I’ve got a first, lowest-order version of the matrix-free solver for the mixed 3d gravity-wave problem working now (pushed to “vertical" branch of firedrake-helmholtzsolver). I precondition with the matrix-free multigrid in pressure space, and it all converges nicely (although it is probably horribly inefficient in places). The equations for velocity u, pressure p and buoyancy b are

du/dt = \grad(p) + b*\hat{z}

dp/dt = -c^2 \div(u)

db/dt = -N^2 u.\hat{z}

Full details on the solver are in GravityWaves.tex in the directory https://github.com/firedrakeproject/firedrake-helmholtzsolver/tree/vertical/doc/latex

However, I’m struggling to incorporate the boundary condition u.n=0 at the bottom and top of the domain (this is a spherical shell in my case - but should all also work for a 2d annulus) into the solver. I presume it is a Dirichlet BC, but of the form n_0*u_0 + n_1*u_1 + n_2*u_2 = 0, and I do not set the values explicity as in, for example, (u_0,u_1,u_2) = (0,0,0). How can I specify this? Also, where in the solver do I actually have to enforce it?

In the mixed solver I apply the mixed operator, so I think I need to enforce it in the u-equation there.

Also, in the mixed preconditioner I do a back-substitution for u after the pressure solve, so I think I have to do it there as well.

Thanks a lot,

Eike

--

Dr Eike Hermann Mueller
Research Associate (PostDoc)

Department of Mathematical Sciences
University of Bath
Bath BA2 7AY, United Kingdom

+44 1225 38 5803
e.mueller@bath.ac.uk
http://people.bath.ac.uk/em459/