Solutions in de Rham complex with boundary conditions.

[ethrelfall]

I wanted to ask how to find FEEC solutions in (what I believe to be) the de Rham complex with boundary conditions.

I am able to find the harmonic mode of the k=1 Hodge Laplacian on e.g. a 2D annulus with the natural boundary condition $u \cdot n = 0$ (see https://github.com/ethrelfall/Finite-element-exterior-calculus/blob/main/eigenvalue_problem_single_vortex.py for my script). This gives the usual vortex $(-y/(x^2+y^2), x/(x^2+y^2))$.

How can I convert the script to find the harmonics subject to the alternate boundary condition $u \times n = 0$? I am looking for the solution $(x/(x^2+y^2), y/(x^2+y^2))$ which is div-free and curl-free. I think the governing equations are the same and the only difference lies in the BCs.

I think that these two sets of solutions are related by Poincar'e duality (i.e. swap the components and change their relative sign) but I wonder if the $u \times n = 0$ one can be found directly in Firedrake.

Thanks,

Ed.

[colinjcotter]

Hi Ed, The alternate boundary condition on the H(curl) space is enforced strongly rather than weakly, i.e. we restrict the finite element space to the subspace satisfying the boundary condition. This is done via DirichletBC just like for the scalar case. The way this works is that you provide a boundary value that has 2D vector shape, and then the nodal variables associated with the boundary are applied to it to set the conditions. So it looks superficially like you are setting the full vector value, but actually only the tangential component is used. The way this is done practically is that we keep the matrix for the unrestricted finite element space, and replace boundary rows with zeros everywhere except for a 1 on the diagonal. This means that there will be a bunch of 1 eigenvalues corresponding to the bcs, which should be ignored in a SLEPc analysis (you are looking for the zero modes, anyway). all the best --cjc

[ethrelfall]

Thanks, I've tried adding e.g.

bc = DirichletBC(V.sub(0), as_vector([0.0,0.0]), {21,22}) #21 ,22 are the inner and outer boundaries

petsc_a = assemble(a).M.handle
petsc_m = assemble(m, bcs=bc).M.handle

(no other changes)

but I get the same vortex solution out of this ... guess I'm doing it wrong?

[colinjcotter]

I think you probably need to apply the bcs in a as well.

[ethrelfall]

Cheers Colin. Tried that but I am getting a boundary mode (eigenvalue = 1) as the lowest mode.

I'll keep experimenting with it, thanks again for the support.

[wence-]

Those are probably the eigenvalues from the identity on the dirichlet nodes. You can shift them away, or extract the sub matrix that doesn't include the bcs and solve the eigenvalue problem for that.

[ScottMacLachlan]

@wence- : what's the recommended way to do either of these? Is there an easy way to extract the submatrix that doesn't include the bcs?

[wence-]

@wence- : what's the recommended way to do either of these? Is there an easy way to extract the submatrix that doesn't include the bcs?

"Easy" is relative. Having assembled the operator you can use:

A = assemble(a).petscmat
Asub = A.createSubMatrix(isrow, iscol)

And go from there, how do you compute isrow and iscol? Something like this:

V = FunctionSpace(...)
bc = DirichletBC(V, ...)
A = assemble(a, bcs=bc).petscmat
localrows = np.arange(V.dof_dset.size, dtype=IntType)
nonbcs = np.setdiff1d(localrows, bc.nodes)
globalnonbc = V.dof_dset.lgmap.apply(nonbcs)
isrows = PETSc.IS().createGeneral(globalnonbc)
Asub = A.createSubMatrix(isrows, isrows)

(Fully untested).

[colinjcotter]

This is probably not it, but you have eps_gen_hermitian but your problem is not Hermitian. On 23 Jul 2022, at 07:26, Cotter, Colin J ***@***.***> wrote:  The point is that Ed is looking for the discrete harmonic 1-forms which span the nullspace. There should be an n-k=2-1=1 dimensional subspace for the operator with zero trace bcs, spammed by a function that looks like a rotating flow going around the torus. But he’s getting the minimum eigenvalue as 1. On 23 Jul 2022, at 00:36, Lawrence Mitchell ***@***.***> wrote:  Those are probably the eigenvalues from the identity on the dirichlet nodes. You can shift them away, or extract the sub matrix that doesn't include the bcs and solve the eigenvalue problem for that. — Reply to this email directly, view it on GitHub<#2501 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4S2YLLOML4GFXNBSLTVVMWATANCNFSM54LVEELA>. You are receiving this because you commented.Message ID: ***@***.***>

[wence-]

The point is that Ed is looking for the discrete harmonic 1-forms which span the nullspace. There should be an n-k=2-1=1 dimensional subspace for the operator with zero trace bcs, spammed by a function that looks like a rotating flow going around the torus. But he’s getting the minimum eigenvalue as 1.

Ah, not paying enough attention, sorry.

[ethrelfall]

Thanks guys.

Removing that "Hermitian" line didn't change what I get.

To be clear, the circulating solution comes out fine, the one I am now trying to find is a purely radial flow (like draining bathtub) and the problem is the u \times n = 0 boundary condition.

"spammed" by a function - :)

[colinjcotter]

Can you try solving a Helmholtz version of the equations (add a mass term to the vector variable and a source term) and check if the bcs are being correctly enforced? Is this a quadrilateral mesh? On 23 Jul 2022, at 12:20, ethrelfall ***@***.***> wrote:  Thanks guys. Removing that "Hermitian" line didn't change what I get. To be clear, the circulating solution comes out fine, the one I am now trying to find is a purely radial flow (like draining bathtub) and the problem is the u \times n = 0 boundary condition. "spammed" by a function - :) — Reply to this email directly, view it on GitHub<#2501 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4SYJ3WXZTMUCDOQPSTVVPIQXANCNFSM54LVEELA>. You are receiving this because you commented.Message ID: ***@***.***>

[ethrelfall]

I can try that but it might not be today (sorry don't think I can match your awesome level of commitment).

It's a fairly coarse triangular mesh (mesh is on my github).

Thanks again.

[ethrelfall]

I've tried writing a simple test for the BC - it looks like the vector Dirichlet condition is applied across both cpts? (NB this is a simpler example with no harmonics and not an eigensolve ...)

https://github.com/ethrelfall/Finite-element-exterior-calculus/blob/main/test_material/square_BC_test_submit.py

Thanks.

[colinjcotter]

The point is that Ed is looking for the discrete harmonic 1-forms which span the nullspace. There should be an n-k=2-1=1 dimensional subspace for the operator with zero trace bcs, spammed by a function that looks like a rotating flow going around the torus. But he’s getting the minimum eigenvalue as 1. On 23 Jul 2022, at 00:36, Lawrence Mitchell ***@***.***> wrote:  Those are probably the eigenvalues from the identity on the dirichlet nodes. You can shift them away, or extract the sub matrix that doesn't include the bcs and solve the eigenvalue problem for that. — Reply to this email directly, view it on GitHub<#2501 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4S2YLLOML4GFXNBSLTVVMWATANCNFSM54LVEELA>. You are receiving this because you commented.Message ID: ***@***.***>

[colinjcotter]

Looks like it, doesn't it! I think that the smoking gun would be a failing convergence test for some test problem with tangential boundary conditions. Easiest way to do it is to cook up a solution, apply vector Helmholtz to it (analytically), and use that as the RHS. all the best --cjc

…

________________________________ From: ethrelfall ***@***.***> Sent: 25 July 2022 13:31 To: firedrakeproject/firedrake ***@***.***> Cc: Cotter, Colin J ***@***.***>; Comment ***@***.***> Subject: Re: [firedrakeproject/firedrake] Solutions in de Rham complex with boundary conditions. (Discussion #2501) I've tried writing a simple test for the BC - it looks like the vector Dirichlet condition is applied across both cpts? (NB this is a simpler example with no harmonics and not an eigensolve ...) https://github.com/ethrelfall/Finite-element-exterior-calculus/blob/main/test_material/square_BC_test_submit.py Thanks. — Reply to this email directly, view it on GitHub<#2501 (reply in thread)>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ABOSV4RRXKYQ7PXJDO7WS4LVV2CKDANCNFSM54LVEELA>. You are receiving this because you commented.Message ID: ***@***.***>