Creating a FESystem with bulk and surface FE spaces

[manuel.q...@gmail.com]

Hello, 

I need to solve a system of many PDEs (currently 18). The actual number changes since the PDEs are defined by a truncated infinite sum and some recursive equations. Anyway, the first two equations are defined on a 2D or 3D domain (call it Omega) and the rest is only defined on the boundary of Omega. My first idea was to create a FESystem object and attach to it 2 FE spaces defined on Omega and 16 FE spaces defined on the boundary. That is, I tried these options 

• FESystem<1, 2> fe(FE_Q<2, 2>(degree), 2, FE_Q<1, 2>(degree), 16)
• FESystem<2, 2> fe(FE_Q<2, 2>(degree), 2, FE_Q<1, 2>(degree), 16)

However, that produced compilation errors. What I understand is that the FESystem must be created by either FE spaces defined in Omega or FE spaces defined in the manifold, but not a combination. As a way around I tried (a very ugly hack) where I defined all FE spaces in Omega and then add constraints to all DoFs associated with the internal nodes for the variables that are defined on the manifold. This sounds like a tremendous waste of resources since most of the equations in my system live in the manifold, so I end up constraining most of the DoFs. And even worse, it produced a matrix that seems to be sometimes singular and sometimes non-singular. That is, I can solve the system and get accurate results (I have the exact solution) for some grids but in other grids the matrix is singular so I get NaNs. I suspect this problem is a consequence of constraining most of my DoFs. FYI, the solver that I am using is the MUMPS implementation in PETSc. 

I wonder if there is a proper or standard way to construct systems that are defined by some equations in a domain and some in a manifold of the domain. I definitely don't think constraining most of my DoFs is the way to go. 

Thanks and all the best, 

Manuel Quezada 

[Wolfgang Bangerth]

Manuel,

> I need to solve a system of many PDEs (currently 18). The actual number
> changes since the PDEs are defined by a truncated infinite sum and some
> recursive equations. Anyway, the first two equations are defined on a 2D or 3D
> domain (call it Omega) and the rest is only defined on the boundary of Omega.
> My first idea was to create a FESystem object and attach to it 2 FE spaces
> defined on Omega and 16 FE spaces defined on the boundary. That is, I tried
> these options
>

> * FESystem<1, 2> fe(FE_Q<2, 2>(degree), 2, FE_Q<1, 2>(degree), 16)
> * FESystem<2, 2> fe(FE_Q<2, 2>(degree), 2, FE_Q<1, 2>(degree), 16)

>
> However, that produced compilation errors. What I understand is that the
> FESystem must be created by either FE spaces defined in Omega or FE spaces
> defined in the manifold, but not a combination.

Yes, that is correct -- because FESystem<dim,spacedim> is derived from
FiniteElement<dim,spacedim>, the elements that form the system must have a
unique 'dim' parameter.

> As a way around I tried (a
> very ugly hack) where I defined all FE spaces in Omega and then add
> constraints to all DoFs associated with the internal nodes for the variables
> that are defined on the manifold. This sounds like a tremendous waste of
> resources since most of the equations in my system live in the manifold, so I
> end up constraining most of the DoFs. And even worse, it produced a matrix
> that seems to be sometimes singular and sometimes non-singular. That is, I can
> solve the system and get accurate results (I have the exact solution) for some
> grids but in other grids the matrix is singular so I get NaNs. I suspect this
> problem is a consequence of constraining most of my DoFs. FYI, the solver that
> I am using is the MUMPS implementation in PETSc.

I don't know about the origins of the illposedness, but this is a valid
approach. ASPECT, for example, does this kind of thing when implementing a
free boundary. You could take a look here
https://github.com/geodynamics/aspect/blob/master/include/aspect/mesh_deformation/free_surface.h
https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc

In particular, ASPECT just doesn't do anything at all with the interior
degrees of freedom in the solve phase:
https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc#L345-L352
In other words, we happily solve a linear system with a matrix with a lot of
zero rows and columns. I've got questions about this here:
https://github.com/geodynamics/aspect/issues/3110
But it seems to work.

> I wonder if there is a proper or standard way to construct systems that are
> defined by some equations in a domain and some in a manifold of the domain. I
> definitely don't think constraining most of my DoFs is the way to go.

The "proper" way to deal with this would probably be to have two DoFHandler
objects with dim and dim-1, and then to build a coupled linear system. I
believe that we have most of the building blocks for this already in place
(take, for example, a look at step-60 and step-70) but there is no concrete
example of coupling bulk and surface problems. It has been on my TODO list for
quite a while already to write a tutorial program for this and see what it
would take to make that work properly.

Out of curiosity, what is the problem you are trying to solve?

Best
W>

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Jean-Paul Pelteret]

Hi Manuel,

As a small addition to what Wolfgang has already suggested, you can also considering using the hp framwwork in conjunction with FE_Nothing elements in order to, effectively, disable elements in the interior of the volumetric domain instead of constraining them directly. This is the equivalent of saying that you simply wouldn’t assign DoFs for those elements, but rather only for the elements that are positioned along the boundary of the domain. You’d still have to constrain the DoFs that are off of the surface itself, but at least there are a lot fewer DoFs that require this.

Best,

Jean-Paul

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
--- You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/6bb61a54-9428-2174-48ae-83b649c1e3a7%40colostate.edu.

[manuel.q...@gmail.com]

Hi Professor, 

Thanks for the quick reply. It is nice to talk to you again. It has been a while since we last spoke. I think it was during an AGU conference in NOLA. 

I don't know about the origins of the illposedness, but this is a valid
approach. ASPECT, for example, does this kind of thing when implementing a
free boundary. You could take a look here
https://github.com/geodynamics/aspect/blob/master/include/aspect/mesh_deformation/free_surface.h
https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc

In particular, ASPECT just doesn't do anything at all with the interior
degrees of freedom in the solve phase:
https://github.com/geodynamics/aspect/blob/master/source/mesh_deformation/free_surface.cc#L345-L352

I took a look at that code and it seems to me that the DoFs that are constrained are only in a boundary and not in the bulk. If this is true, the number of constrained DoFs is significantly smaller than in my problem. 

 

In other words, we happily solve a linear system with a matrix with a lot of
zero rows and columns. I've got questions about this here:
https://github.com/geodynamics/aspect/issues/3110
But it seems to work.

That is an interesting discussion. So, is adding a constraint via add_line setting the whole line with zeros? I thought that it was adding zeros for all columns in the i-th row and a 1 in the diagonal. I should have read the documentation more carefully. Of course, I agree with the discussion in issue #3110. It is surprising that this can be done and that it works well in ASPECT. I strongly suspect this is the source of my problem. The matrix that I get behaves very weird. For a given mesh, I can solve the system and the solution is accurate, but if I consider a slightly different mesh the system becomes singular. The whole matrix (by zeroing out entire rows) is technically singular, so I can't blame the solver for not being able to solve it properly for all cases. And even more when the vast majority of the rows are filled with zeros. If I am right, why not simply add a 1 in the diagonal via add_entry? I am sure there is a strong reason for not doing this; otherwise, you would be doing it in ASPECT.  

> I wonder if there is a proper or standard way to construct systems that are
> defined by some equations in a domain and some in a manifold of the domain. I
> definitely don't think constraining most of my DoFs is the way to go.

The "proper" way to deal with this would probably be to have two DoFHandler
objects with dim and dim-1, and then to build a coupled linear system. I
believe that we have most of the building blocks for this already in place
(take, for example, a look at step-60 and step-70) but there is no concrete
example of coupling bulk and surface problems.

Yes, I thought about this but I was hoping to find an easier way. I will for sure consider this as an option. Thanks for the suggestion. I also thought on decoupling the system by having two DoFHandler objects and do some type of fixed point iteration where I use known information from the boundary (or an initial guess) when solving the equations in the bulk and similarly with the equations in the boundary. Of course the problem is that I am not sure if that iterative process would converge. I am also concerned that this might add some instability similar to the "added mass" effect in FSI problems. 

 

It has been on my TODO list for
quite a while already to write a tutorial program for this and see what it
would take to make that work properly. 

 

Out of curiosity, what is the problem you are trying to solve?

I am solving the Helmholtz equation to model scattering. The interesting part of this particular problem is that it is possible to apply arbitrarily accurate absorbing boundary conditions by including some boundary terms that come from a Karp's expansion of the exact solution in the exterior of the domain (the far field solution). This expansion makes the system of PDEs considerably more complicated (see the attached notes.pdf), but since the absorbing boundary condition is arbitrarily accurate (if I include enough terms in the expansion), I can use high-order FEs and observe the expected rates of convergence. This is true for artificial external boundaries that are (potentially) very close to the scatter. This has been done with IGA (see the attached PDF). At the moment I am collaborating with one of the authors of that paper to do a "standard" FE implementation that can reproduce some of their results and then we plan to move to a more complicated situation. 

Once I have a working code that reproduces some of the results in the attached paper, I can share the code with you. If you think the problem and code are interesting for a tutorial then we can work together in such a tutorial. It might be an interesting problem for a tutorial since we solve a system of any number of PDEs (by only changing the number of terms in the expansion) some of which are in a domain and others in the boundary. 

Bests, 

Manuel 

[manuel.q...@gmail.com]

Hi Jean-Paul, 

I have taken a look at FE_Nothing and hp::DoFHandler. If I understand correctly, the idea is that all the interior cells will use an FE_Nothing space and the boundary cells will use a standard FE_Q space, right? My question is what happens to the DoFs that are shared by a cell with a FE_Nothing space and another cell with a FE_Q space? This is of course going to happen at all cells in the boundary since some DoFs of those cells are indeed in the boundary while others are still internal DoFs. Does this make sense? I suppose I can apply a constraint to those internal DoFs, similar to my previous "hack". The advantage following your suggestion is that the number of DoFs that I need to constrain is much smaller. 

Thanks for your suggestion. I think this is a very nice approach, 

Manuel 

[Andrew McBride]

Hi Manuel,

I’m not sure if this helps, but some time ago we developed a code for surface elasticity in deal.ii. In this problem the surface of a bulk has its own energetic structure. One effectively constructs system matrices for volume cells and surface cells and these contribute to the same global dofs.

For more details see https://github.com/mac-a/surface_elasticity and https://arxiv.org/abs/1506.01361 . This is 6 years old so I can’t promise it will work but it might be of some use.

Best,

Andrew

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/b09a005d-8849-4c70-a851-1997aa398ec5n%40googlegroups.com.

[manuel.q...@gmail.com]

Hi Andrew, 

This definitely helps! Although I still have to check the code in detail, I found and scanned through your paper a few days ago (before posting my question here). I went to the link in the paper to find and look at the code, but the link is broken so I decided to ask for help here. It is nice that now I have access to your code! Thanks a lot. 

Bests, 

Manuel 

[Jean-Paul Pelteret]

HI Manuel,

I’ve got a brief answer to your followup questions, but let me first say that Andrew’s suggestion, which seems to be the most elegant/optimal approach to the problem, is one that is possible to for you to realise. It’s also made me think that perhaps some of the tools that have been added somewhat recently to help tackle immersed problems could be of help? I have zero experience with these tools, so I’m just going to put it out there for you to investigate and draw your own conclusions from:

https://dealii.org/current/doxygen/deal.II/namespaceNonMatching.html

https://dealii.org/current/doxygen/deal.II/step_60.html

I’m sure that if you have any questions on those, then the other developers would be willing to help out to answer them!

My question is what happens to the DoFs that are shared by a cell with a FE_Nothing space and another cell with a FE_Q space? This is of course going to happen at all cells in the boundary since some DoFs of those cells are indeed in the boundary while others are still internal DoFs. Does this make sense? I suppose I can apply a constraint to those internal DoFs, similar to my previous "hack”.

Yes, that’s correct. You’d be doing exactly the same thing as you’d detailed before, but just for fewer elements / DoFs.

The advantage following your suggestion is that the number of DoFs that I need to constrain is much smaller. 

Indeed. So not only is your global linear system smaller, but the AffineConstraints object also has less entries. Why this might be important is due to a comment that can be found in the AffineConstraints introductory notes:

The class is meant to deal with a limited number of constraints relative to the total number of degrees of freedom, for example a few per cent up to maybe 30 per cent; and with a linear combination of M other degrees of freedom where M is also relatively small (no larger than at most around the average number of entries per row of a linear system). It is not meant to describe full rank linear systems.

Its hard to say whether this might have a performance effect wrt. what you’re trying to do, but nevertheless using Fe_Nothing to limit the number of DoFs that you have to constrain would likely mean that you could consider finer discretisation before you have to start considering what the impact of the number of constraints have on performance (and anything else).

Best,

Jean-Paul

To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/b09a005d-8849-4c70-a851-1997aa398ec5n%40googlegroups.com.