[aspect-devel] Questions about surface process

[Wolfgang Bangerth]

Yimin,

> I am a student of University of Chinese Academy of Sciences and I have been
> using ASPECT to do simulations on crack propagation of small-scaled geological
> bodies in the passing year. I have found that the existing velocity pojection
> strategies in the free-surface module introduces too much error to the result
> (the 'normal' strategy leads to severe distortion of cells, and the 'vertical'
> strategy ignores horizontal displacements), while the sticky-air method will
> double the computational complexity of the problem (my model can 'grow' tall
> so the air layer must be thick enough). Therefore I intended to apply a
> diffusive surface process on the free-surface, but when I looked into this
> problem I found it quite difficult because the partition changes each time we
> execute refinement and coarsening. The finest way seems to be transporting
> faces between processes using functions provided by p4est (I don't know if
> there are any), then it will be unavoidable to read the codes of p4est. Could
> you please tell me whether you have started to do such things or not?  If it
> has been implemented and will soon be released, then I will be happily looking
> forward to the next version of ASPECT; otherwise, I would like to make my
> effort in this part of work.

Let me answer this whole thread at once, or at least those that John didn't 
already answer. I'm quite impressed how far you had already gotten. You 
definitely put the finger where it hurts.

As you already noticed, the difficulty in merging bulk and surface processes 
is the creation of a surface mesh that has the same partitioning as the subset 
of cells that are at the surface. This would work through the 
extract_boundary_mesh(), which of course is not implemented for 
parallel::distributed::Triangulation objects.

The reason for this is that p4est does not allow user control over the 
partitioning. As a consequence, it would be possible to create the *coarse* 
mesh of the surface mesh on every processor, because every processor knows the 
entire coarse mesh of the volume mesh. But that's where things break down: if 
you look at the implementation of the function in 
source/grid/grid_generator.cc, you'll realize that the boundary mesh has to be 
recursively refined based on the volume mesh, but if you can not control the 
partitioning of the boundary mesh, this will not be possible. (I'm sure one 
could come up with a scheme that actually creates this boundary mesh, at the 
cost of a significant amount of communication. But, importantly, you will 
still end up with two meshes where the same processor will not at the same 
time have access to cells and their boundary faces -- because of the 
partitioning issue --, and that will prevent you from *using* this boundary 
mesh in a way that couples bulk and surface processes because you can't 
exchange data efficiently between the two meshes.)

I've thought a bit about alternatives over the past couple of days. I think 
the best -- if maybe not most scalable approach is the following:
* deal.II has a class called parallel::shared::Triangulation. For this
   kind of triangulation, *every* processor has a copy of the entire mesh.
   That limits how many cells you will be able to deal with, but since surface
   meshes typically have far fewer cells than the volume mesh, this should be
   good enough to go to a few hundred processors or so.
* Because every process has the entire boundary triangulation, it should
   be relatively straightforward to write the equivalent of
   GridGenerator::extract_boundary_mesh that takes a p::distributed::Tria
   as argument and returns a p::shared::Tria as result. Internally, it would
   basically work exactly the same as the existing function, with the
   following exception: in each iteration, every process would mark those
   cells for refinement of which it knows, and then there needs to be
   a communication step in which processors exchange their lists
   and generate a combined list of boundary cells to refine. If you
   represent the individual lists as vector<bool>, then the parallel
   combination is essentially just an 'or' operation.
* At the end of this process, every process has a complete copy of the
   parallel::shared::Triangulation, and you can set the subdomain_id
   of all cells in such a way that the partitioning of the boundary mesh
   matches the one of the adjacent volume cells. Some processors will
   simply not own any boundary cells.

With all of this in place, you will be able to define equations that then 
couple volume and surface terms.

If you'd like to follow this direction, please let us know. The function I'm 
outlining above would be independently useful for deal.II, and we'd be very 
happy to walk you through the process of writing patches for deal.II!

Best
  Wolfgang

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