Re: [firedrake] BCs on velocity for mixed gravity wave system solver
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/...
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/
Hi Andrew, yes, this should be sufficient in my case, my velocity space is always HDiv(U2xV0) + HDiv(U1xV1), where the second part is the vertical velocity. At lowest order I use so far U1 = RT_0, U2 = DG_0, V0 = CG_1, V1 = DG_0, but I also want to work with the next-to-lowest-order case U1 = BDFM_1, U2=DG_1, V0 = CG_2, V1 = DG_1. I will try out what you suggest. Thanks a lot for your help, 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/
On 27 Dec 2014, at 15:33, Andrew McRae <a.mcrae12@imperial.ac.uk> wrote:
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 <mailto: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/... <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 <mailto:e.mueller@bath.ac.uk> http://people.bath.ac.uk/em459/ <http://people.bath.ac.uk/em459/>
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
I've been wondering if we should make a specific boundary condition type that is only permitted with H(div) spaces, and sets the normal component from a scalar. We discussed it over coffee with David and Lawrence but I can't remember what the conclusion was. Cjc On 28 Dec 2014 12:21, "Eike Mueller" <e.mueller@bath.ac.uk> wrote:
Hi Andrew,
yes, this should be sufficient in my case, my velocity space is always HDiv(U2xV0) + HDiv(U1xV1), where the second part is the vertical velocity. At lowest order I use so far
U1 = RT_0, U2 = DG_0, V0 = CG_1, V1 = DG_0,
but I also want to work with the next-to-lowest-order case
U1 = BDFM_1, U2=DG_1, V0 = CG_2, V1 = DG_1.
I will try out what you suggest. Thanks a lot for your help,
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/
On 27 Dec 2014, at 15:33, Andrew McRae <a.mcrae12@imperial.ac.uk> wrote:
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/...
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
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
participants (3)
-
Andrew McRae
-
Colin Cotter
-
Eike Mueller