Dear all,

I think I fixed this now, and the gravity wave system is solved with the Dirichlet BC n.u=0 at the top- and bottom boundary of the spherical shell. Since PETSc uses left-preconditionining by default, all I have to do is enforce the boundary condition on u in the preconditioner for the mixed system:

Constructor of GravityWave solver class:

[…]
self._bcs = [DirichletBC(self._W2, 0.0, "bottom"),DirichletBC(self._W2, 0.0, "top”)]
[…]

solve() method of mixed preconditioner class:

[…]
# Pressure solve
p.assign(0.0)
self._pressure_solver.solve(self._rtilde_p,p)
# Backsubstitution for velocity 
assemble(self._dt_half * div(self._utest) * p*self._dx,tensor=self._tmp_u)
self._tmp_u += self._rtilde_u
self._mutilde.divide(self._tmp_u,u)
for bc in self._bcs:
    bc.apply(u)
[…]

Since the solution is in the Krylov space span{r',(P^{-1}A).r',…,(P^{-1}A)^k.r'} with r' = P^{-1}r, it does automatically fulfill the BC.

I checked that the corresponding nodes are indeed set to zero by looking them up via bc.nodes.

Thanks a lot again,

Eike



On 27 Dec 2014, at 15: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/

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake