How to allow incompatible velocity boundary conditions in incompressible models?

[gassmoeller]

This is a problem that came up in one of @johnjarmitage 's models. ASPECT currently applies a pressure right-hand side compatibility modification for compressible models, but not for incompressible models (source/simulator/core.cc:628). This is based on the reasoning that the right hand side should be compatible in incompressible models anyway, because there should be no velocity convergence.

Now what about the case where we have incompatible boundary conditions that we can not change (e.g. imprecisely interpolated velocities from another computation)? In this case incompressible models currently fail (depending on the magnitude of the overall velocity divergence imposed by the velocity boundary conditions and the solver tolerance). This can be seen in the attached slightly modified platelike-boundary cookbook. I impose a small y- velocity at the top boundary, all other boundaries are closed. The remaining Stokes residual linearly correlates with the magnitude of the influx, and already small influx velocities (1e-3 of the average velocities of the model) prevent convergence. Enabling the pressure right hand side compatibility modification immediately leads to convergence.

Should we add a parameter to enforce a pressure right hand side compatibility modification for incompressible models? Or should we somehow automatically detect this case and then enable the modification? I do not think the current way of simply throwing a non-convergence error is an acceptable solution, because while there should not be any velocity divergence in incompressible models there are models that need this kind of setup.

failing_example.prm.txt

[naliboff]

My vote would be to add a parameter that enforces a tight right hand side compatibility rather than do it automatically. I think by default the user should be aware if the boundary mass fluxes do not balance in an incompressible model, which is currently the case through the model not converging.

A scenario to keep in mind: applying the pressure modification to an incompressible model with a high velocity divergence could produce odd dynamics if a pressure-dependent rheology is used. While it is up to the user to make sure this is not the case, perhaps it is worth putting in a cut-off for the maximum allowable pressure modification (or velocity divergence value)?

[tjhei]

Physically, enabling it means that we add a constant sink (or source) of material over the whole domain. That is not something I would think is sensible to have on by default. In most cases it is probably more appropriate to leave one of the boundaries open.

Should we add a parameter to enforce a pressure right hand side compatibility modification for incompressible models?

That is okay with me.

Or should we somehow automatically detect this case and then enable the modification?

You can compute the correction factor and then decide what to do based on the value (error / warning / apply correction). I am not exactly sure how to make this comparison robust though.

[bangerth]

My preference is to do this:

• In the place where we currently interpolate the boundary values for the velocity, if all boundaries are Dirichlet, then also compute the net flux into/out of the domain.
• If this value is nonzero to some non-negligible degree (maybe tied to the solver tolerance somehow?) then error out and tell the user that that's not going to work.
• Possible, we could introduce a flag that a user could then set (but that is off by default) that would compute some correction. That could be done on the incompressibility, but it could also be done on the boundary values. I might think that the latter might actually be more desirable.

[gassmoeller]

In the place where we currently interpolate the boundary values for the velocity, if all boundaries are Dirichlet, then also compute the net flux into/out of the domain.

Sounds reasonable. Modifying the boundary condition is more accurate, and also relatively simple for incompressible models, because we only need to worry about the boundary-integrated velocity, and not rho*u. I will try to come up with something over the next weeks.

[bangerth]

I think the question is how exactly you want to modify the boundary velocity. The first step ought to be to integrate the net flux through the boundary. This needs to use the interpolated boundary velocity that we also use to compute constraints, i.e., F = \int_{\partial \Omega} (I_h g) \cdot n where I_h is the interpolation operator. Now that you know the total flux F, you somehow need to replace I_h g by I_h g - F / |\partial\Omega| n -- not the factor of n in there. It may be that this just involves substituting the rhs of the constraints by values that are slightly different (i.e., not having to call the function that computes the constraints again), but it's not entirely clear to me how and where to do this.

[tjhei]

This reminds me of http://dx.doi.org/10.1016/j.jmaa.2016.01.075 where we correct boundary conditions (in our case) caused by the boundary interpolation (not incompatible boundary conditions though).