{"action":"created","issue":{"url":"https://api.github.com/repos/geodynamics/aspect/issues/643","labels_url":"https://api.github.com/repos/geodynamics/aspect/issues/643/labels{/name}","comments_url":"https://api.github.com/repos/geodynamics/aspect/issues/643/comments","events_url":"https://api.github.com/repos/geodynamics/aspect/issues/643/events","html_url":"https://github.com/geodynamics/aspect/issues/643","id":111396109,"number":643,"title":"Depth of geometry models with free surface not consistent after first time step","user":{"login":"MFraters","id":7631629,"avatar_url":"https://avatars.githubusercontent.com/u/7631629?v=3","gravatar_id":"","url":"https://api.github.com/users/MFraters","html_url":"https://github.com/MFraters","followers_url":"https://api.github.com/users/MFraters/followers","following_url":"https://api.github.com/users/MFraters/following{/other_user}","gists_url":"https://api.github.com/users/MFraters/gists{/gist_id}","starred_url":"https://api.github.com/users/MFraters/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/MFraters/subscriptions","organizations_url":"https://api.github.com/users/MFraters/orgs","repos_url":"https://api.github.com/users/MFraters/repos","events_url":"https://api.github.com/users/MFraters/events{/privacy}","received_events_url":"https://api.github.com/users/MFraters/received_events","type":"User","site_admin":false},"labels":[],"state":"open","locked":false,"assignee":null,"milestone":null,"comments":7,"created_at":"2015-10-14T13:16:39Z","updated_at":"2015-10-21T13:50:41Z","closed_at":null,"body":"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. \r\n\r\nThe 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. \r\n\r\nThe 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.\r\n\r\nA 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.\r\n\r\nPersonally 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?\r\n\r\nHere 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).\r\n\r\n![depth_problem](https://cloud.githubusercontent.com/assets/7631629/10484491/9abc43ea-7285-11e5-833b-f157d22ece58.png)"},"comment":{"url":"https://api.github.com/repos/geodynamics/aspect/issues/comments/149902070","html_url":"https://github.com/geodynamics/aspect/issues/643#issuecomment-149902070","issue_url":"https://api.github.com/repos/geodynamics/aspect/issues/643","id":149902070,"user":{"login":"bangerth","id":7504421,"avatar_url":"https://avatars.githubusercontent.com/u/7504421?v=3","gravatar_id":"","url":"https://api.github.com/users/bangerth","html_url":"https://github.com/bangerth","followers_url":"https://api.github.com/users/bangerth/followers","following_url":"https://api.github.com/users/bangerth/following{/other_user}","gists_url":"https://api.github.com/users/bangerth/gists{/gist_id}","starred_url":"https://api.github.com/users/bangerth/starred{/owner}{/repo}","subscriptions_url":"https://api.github.com/users/bangerth/subscriptions","organizations_url":"https://api.github.com/users/bangerth/orgs","repos_url":"https://api.github.com/users/bangerth/repos","events_url":"https://api.github.com/users/bangerth/events{/privacy}","received_events_url":"https://api.github.com/users/bangerth/received_events","type":"User","site_admin":false},"created_at":"2015-10-21T13:50:40Z","updated_at":"2015-10-21T13:50:40Z","body":"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.\r\n\r\nI 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.\r\n\r\nWhether 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."}}