Fixed compositional field (and temperature) boundary conditions on deforming boundaries

[anne-glerum]

By default, fixed temperature boundary conditions are prescribed on both the inflow and outflow parts of the user-selected boundaries. For compositional fields, the default is to only prescribed values on inflow parts. To determine what is an inflow or outflow part of the boundary, the velocity at the boundary is queried in the replace_outflow_boundary_ids() function in helper_functions.cc. This can be done by asking the prescribed velocity plugin for a value, or, in case of a deforming boundary, by extracting the velocity at the location in question from the current_linearization_point.

However, for deforming boundaries, the velocity at the boundary moving up or down does not necessarily correspond to actual outflow or inflow. For example, a free surface moves up and down according to (a projection of) the velocity at the boundary, without any in/outflow. A boundary that is deformed by FastScape, is moved by uplift and by erosion or deposition. In this case, an uplift ("outflow") could coincide with deposition of material, for which a boundary condition would be appropriate. Uplift plus erosion would I think not require a boundary condition to be prescribed. Depending on the relative magnitudes of uplift and erosion, the new surface height could be larger or smaller than the previous one, but the value on the boundary should reflect the material that was eroded above it, not a prescribed value.

In the advection equations, the mesh velocity is subtracted from the material velocity first. @gassmoeller suggested to do a similar thing to determine in/outflow on deforming boundaries. For a free surface using the normal projection, I think this would work (depending on the projection used). It would lead to a zero velocity at the surface, and thus no outflow boundaries.
For other mesh deformation methods, I can't fully wrap my head around it yet, but I think it could work. For example, when an uplifted area is eroded more than it is uplifted, the mesh velocity would be negative and the boundary would be considered an outflow boundary (as I think it should). When uplift and deposition occur, the boundary would be considered an inflow boundary, which I think is also correct. Any thoughts on this solution would be appreciated!

[anne-glerum]

I've given this some more thought and set up a benchmark based on the /benchmarks/diffusion_of_hill/1_sine_zero_flux.prm setup (see attached):

• Gravity is turned off.
• All boundaries are free slip, except for the top which is prescribed zero velocity (otherwise no in/outflow is computed).
• Hence there are no driving forces and Stokes velocity is zero.
• The initial topography describes a hill spanning half of the top boundary.
• The top boundary is deformed by the diffusion mesh deformation plugin.
• There are three compositional fields.
• Only one field (layer_1) occupies part of the hill at t0, other fields are 0 at t0.
• The boundary condition for layer_1 is 1 at t0 and 0 afterwards. Ideally, the BC would be the initial composition.
• The boundary condition for layer_3 is 0 at t0 and 1 afterwards.

Due to diffusion, I expect the hill to be eroded away, and the rest of the domain to fill up. The eroding hill should stay layer_1, the deposited material should be layer_3.

I made 2 movies for both main and for my branch adjust_v_for_outflow_BC showing the evolution of layer_1 and layer_3. In the latter branch, I remove the mesh velocity from the material velocity when computing whether outflow occurs along the top boundary. The behaviour of layer_1 and layer_3 is then as I would expect.

Note that I had to set the whole boundary to outflow at t0 and t1, because:

1. The compositional boundary conditions are computed at the beginning of the timestep before mesh deformation.
2. Mesh deformation is first computed at t1, not t0.
3. Stokes and mesh velocity are therefore 0 during t0 and t1, and the whole top boundary would become a prescribed field boundary. That means that layer_1 and layer_3 both become 1 along the top boundary otherwise.

Movies: https://nextcloud.gfz-potsdam.de/s/8eHNbSxCcs7Z5e5
prm file attached: 1_sine_zero_flux_ci_main.txt

On a side note (I should make this a separate issue): if the mesh has adaptive mesh refinement (coarser cells at the bottom of the domain), the mesh velocity along the top boundary shows anomalies at the locations of the edges of the largest cells. They disappear when a uniform mesh is used.