Depth of geometry models with free surface not consistent after first time step

[MFraters]

The way the depth is currently calculated in all the geometry models is after the first time step only correct if the free surface is not used. When the free surface is used, the topography is not accounted for, resulting in a incorrect depth. This may be relevant for for example depth dependent material models and depth dependent mesh refinement.

The challenge in fixing this problem is that the change of the geometry has to be calculated and interpolated. There are several solutions possible, and I am wondering what your thoughts on this matter are.

The first solution would be to place al responsibility with the geometry models. It is possible that they like the topography post-processor loop over all the surface elements and interpolate the topography between those points to calculate the correct depth.

A second solution can be to share the responsibility between the geometry model and the free surface code. The free surface code would calculate and interpolate it (or do something fancy with a manifold as mentioned in https://github.com/geodynamics/aspect/issues/369) and the geometry model only has to ask what the topography component of the depth is and account for it.

Personally I think the second option is the preferred one since it prevents a lot of code duplication, safes a lot of computation if the free surface code computes the topography once everytime the free surface is updated and keeps the geometry models clean. Or is there a better/simpler way?

Here is a picture of the problem (to make this I needed to comment out the asserts in the box geometry model, because when I ask for depth, it can't deal with a negative depth).

[ian-r-rose]

This is tricky. For a significantly deformed mesh, it is not even immediately clear what depth means, or whether it has much physical significance. Depending upon the application, I could see it meaning one of several things:

• The distance in the direction of gravity below the surface (I think what you are intending here)
• The distance in the direction of gravity below the reference surface (perfectly fine for small deformations, and what is implicitly done already)
• A "hydrostatic" depth, given by converting between the pressure solution to a depth given a reference gravity and density: z = P/(rho_0 g_0)

Part of the reason that the suggestions you are making are complicated is that (as far as I can see) you would need to construct some sort of function that gives depth, which is dependent on both the current (deformed) shape of the mesh and the gravity model. Something like depth = f (mesh, gravity, position) The first argument of f is a problem: not only is the calculation of depth kind of tricky (you would need to shoot in the direction of gravity until you hit the surface), but it is also a nonlocal operation. Since the mesh is a distributed object, the processor which knows the position of the surface is in general not the same as the processor which wants to know the depth at a point.

These considerations are not really a solution to your issue, but they do illustrate the challenges. Is there a particular depth-dependent module that you are unable to use at the moment?

[bangerth]

As an aside, the "hydrostatic" depth, I believe, is related to the geodesic related to following a curve from the point to the surface along which the gravity vector is always tangent. Because the gravity vector is, in general, different everywhere, this curve need not be a straight line from the given point to the surface. Of course, in practice, this is going to be very difficult to compute.

I, too, think that computing such a depth may be interesting to do, but will be incredibly difficult to implement efficiently in a way that allows for its use in all of the places where we currently use it for the depth.

[MFraters]

@ian-r-rose The first option is indeed what for me would be the logical thing what I would call depth. The second option is indeed the approximation we are now using, but with 'large' deformations (whatever that means) this is not valid any more. The third option is if I understand it correctly indeed related to the first option, and should result in the ideal case in exactly the same way. Calculating this is not trivial though, because we only calculate the total pressure, and not the adiabatic/hydrostatic pressure (although there is a function is aspect which says it can differentiate between the adiabatic and nonadiabatic pressure or something like that, but I am not completely sure what it exactly does.)

I do understand it is not a trivial exercise, but I do not understand why it should be so difficult and inefficient. What in my eyes should happen is the following:

1. the following points should only happen if the free surface module is used, otherwise use the old thing. This makes sure no efficiently is lost when the free surface is not used.
2. for every update of the mesh, loop over all the surface points and store the deviation from the original surface. This should be trivial to generate since the topography post processor does already something similar.
3. Then every time the depth is needed, we only need to interpolate between the 2/4 points at the surface surrounding this. There are two parts which are tricky:
   1. It will be tricky to do the interpolation efficient.
   2. We need to know the location of a point at the surface to be able to interpolate. Although we might officially need to use the gravity vector direction, I think that the geometry models implicitly force people to use certain directions for the geometries. For example, for the box model, people will only use vertical, and for spherical cases they will use something spherical. It makes in my eyes no sense to use a spherical gravity for a box model. The magnitude may of course change to their preference, but the direction will most likely stay the same. So in the case of the box, I know where my point on the surface is (except for depth), and also for the chunk I also know where the point on the surface with with the help of the manifold (again except for the depth).

If we do the interpolation within a function in the free surface code (where the topography is also stored), then the only thing the geometry model then needs to do is call a function freesurface::get_topography(position) to get the topography, since the free surface code has already direct access to the manifold if needed. This topography can then be used to correct the depth.

Am I missing something?

[ian-r-rose]

One big problem is still that the topography, in general, lives on a different processor than the point at which you want to query the depth. That is going to require a lot of all-to-all communication, which should be avoided if at all possible.

[MFraters]

Hmm, good one. I am not that experienced with mpi programming, but I guess that we would only need to share the topography to every node, so every node has it's own copy of the topography. I do not expect this to take a lot of memory (since it is only one double value per surface point). In this way the communication would be limited to toe moment the grid is changed.

[jperryhouts]

I think one of the big goals of parallelizing code is supposed to be that the local memory usage is independent of the model size (ie model size / number of cpus or something). If you have to store information about the surface on every node, it might make sense to do it in a fixed-size vector which is its self already an interpolation of the real surface.

That said, I think @bangerth's comment is a good one, that to do this right, finding the "depth" would involve integrating the length of the geodesic tangent everywhere to the gravity field until you reach the surface, which sounds obviously more complicated than almost anyone would care to waste CPU cycles on. Maybe in practice this isn't super important, but I just don't know.

I personally think the "hydrostatic" depth should be fine for most applications. That way code can safely assume that depth is always positive, it has an obvious physical meaning, and it doesn't involve any awkward MPI kludges.

On 10/15/2015 09:13 AM, MFraters wrote:

Hmm, good one. I am not that experienced with mpi programming, but I guess that we would only need to share the topography to every node, so every node has it's own copy of the topography. I do not expect this to take a lot of memory (since it is only one double value per surface point). In this way the communication would be limited to toe moment the grid is changed.

— Reply to this email directly or view it on GitHub https://github.com/geodynamics/aspect/issues/643#issuecomment-148441328.

[bangerth]

I don't think that computing the "real" depth, taking into account free surfaces, is feasible in any kind of efficient way by tracing back to the surface for individual points. And we cannot replicate the topography to all processors either: if, for example, you have 10,000 processors, then O(1000) of those are at the surface (in 3d) and that's too much data to store on every processor.

I think the closest you can get is to not compute the depth for every point, but instead to note that the depth d(z) satisfies the differential equation d'(z)=1, d(0)=0. You can write this as a partial differential equation for points x,y,z in every geometry: \vec g/||\vec g|| \cdot nabla d(x,y,z) = 1, d=0 at the surface. This is a transport equation that we can solve numerically, and it would give us a finite element field d_h(x,y,z) that at every point gives us an approximation to the depth. The equation can of course be solved in parallel as well.

Whether that's worth the effort I don't know. @MFraters -- what do you actually need this for? I'd rather document what the depth() function does than change fundamentally how it operates.

[MFraters]

I guess that that is indeed the most correct way to calculate it. For now implementing this method isn't worth the effort yet for me, so I guess documenting this problem it should be fine. Maybe we should also put an assert in the depth calculation, to prevent it from being used after the initial timestep when the free surface is on, although for small deformations it might not be much of a problem.

The reason I am interested in the depth is that I want to test rheologies which are depth dependent (so not using the full pressure, as we are doing now), because they are said to converge quite a bit better. But I can also imagine that doing the refinement by depth or the boundary conditions might also be useful.

[bangerth]

I see. But as I said, I don't think it is easy enough to compute.

I started updating the documentation, but need to run now. Will get back to it tonight.

[bangerth]

See #649 for improved documentation in this area.