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