geodynamics / aspect

A parallel, extensible finite element code to simulate convection in both 2D and 3D models.
https://aspect.geodynamics.org/
Other
227 stars 237 forks source link

How to call constraints.distribute() in solve_stokes(). #5412

Closed bangerth closed 1 year ago

bangerth commented 1 year ago

While looking at #5398 (see also #5405), I came to realize that solve_stokes() has a design problem that I don't quite know how to work around.

At its core, this is the situation: In solve_stokes(), we create a distributed_stokes_solution vector that has only one or two blocks (and so represents a subset of the blocks and degrees of freedom of a complete solution vector), but then we call

  current_constraints.distribute(distributed_stokes_solution);

where current_constraints is an AffineConstraints object that represents constraints on all constraints on DoFs. (current_constraints is initialized in compute_current_constraints(), which fills it with a set of constraints for all variables, not just a subset, such as the Stokes variables.)

The fact that this used to work is a bit of a miracle before because one could imagine that there are constraints also for other variables and I have no idea what AffineConstraints::distribute() is supposed to do with constraints for DoFs that simply are not in distributed_stokes_solution because the latter is a truncated vector. This looks like a pre-existing bug to me, even though it wasn't caught with an assertion.

The question is how to deal with this in solve_stokes() now that

  current_constraints.distribute(distributed_stokes_solution);

is simply no longer legal. I can see a number of options:

Opinions on these matters? I must admit that I'm a bit stuck :-)

tjhei commented 1 year ago

I vote for option 1 and I think you can use add_selected_constraints to do what you want: https://www.dealii.org/developer/doxygen/deal.II/classAffineConstraints.html#a1a2dae8cbd98611a6cfbbe7c0281abcb

bangerth commented 1 year ago

Thanks for pointing out add_selected_constraints. It turns out that that function has serious problems, so I deprecated it but took the majority of the implementation for what I do in https://github.com/dealii/dealii/pull/16099.