Allow plugins to modify sparsity pattern

[maxrudolph]

This allows for constraints that can potentially modify the sparsity pattern to be evaluated prior to the creation of the system matrices. This allows for prescribed pressure, for instance, as required for Pull Request #1864 .

[tjhei]

I think we should be calling compute_current_constraints() and use current_constraints when setting up the sparsity pattern. Then the signal would be triggered automatically already and your change is not necessary.
It is a lucky accident that the current system works, because we are missing the information about constraints that change over time when setting up the matrix. It works, because all diagonal blocks (except for the pressure, which is why you got the error) are marked as "coupled" so the diagonal is correctly in the pattern. This means we have additional entries in the sparsity pattern that are unnecessary.
@bangerth do you agree? /run-tests

[maxrudolph]

@tjhei would it work to just call compute_current_constraints() after constraints.close() and then constraints.merge(current_constraints)?

[bangerth]

I don't think this patch can be correct.

We call compute_current_constraints() at the top of solve_time_step(). But you are already making use of current_constraints in the place where we build the sparsity pattern -- i.e., in the function setup_system_matrix() (which is where all of your changes are), which is called from setup_dofs() which is called at the beginning of run(). In other words, at the point where you are calling setup_dofs(), current_constraints is not yet initialized.

I think a better way is to find a place to hook into the mechanism that builds the constraints. At the point where current_constraints is initialized, the sparsity pattern is already built. (One could argue that that is a bug -- current_constraints can only differ from constraints by constraints that do not change the sparsity pattern of the matrix; that appears to be the case in all uses we've had so far.)

[maxrudolph]

@bangerth would the approach suggested above work:

call compute_current_constraints() after constraints.close() and then constraints.merge(current_constraints)?

[bangerth]

I don't think it's that easy. compute_current_constraints() is (and must be) called once per time step, but the sparsity pattern is built only once before the time loop even starts. So current_constraints is just the wrong place to put something in that you need for the sparsity pattern. I'd just create a new signal that's called right after we create the constraints right now but before we create the sparsity pattern. A listener to the signal can then add new constraints and they will propagate into the sparsity with some luck. (Though I'm not sure that's the case for identity constraints DoFs put into the ConstraintMatrix -- but that still seems like the right direction to look.)

[tjhei]

So current_constraints is just the wrong place to put something in that you need for the sparsity pattern.

I don't understand. Of course current_constraints might be changing and we haven't specified if it is allowed to change the type of constraint over time or only the values. But I assume that current_constraints computed at t=0 is more correct to use than constraints for make_sparsity_pattern.

This patch computes current_constraints before using them. I think that this change not only allows us to add additional constraints (like the pressure dofs), but will also improve memory usage and performance, because we carry around plenty of unused sparsity pattern entries for boundary conditions that are only in current_constraints.

[maxrudolph]

@tjhei I did notice that in some of the failed tests are showing fewer nonzero entries. Interesting.
Did you look at the failed tests? Some of them look like they differ numerically only at the level of roundoff (e.g. composition_active_with_melt) but others looks like they differ more substantially (e.g. inner_core).

[tjhei]

I did notice that in some of the failed tests are showing fewer nonzero entries. Interesting.

This is what I predicted. ;-)

but others looks like they differ more substantially (e.g. inner_core).

I wouldn't worry about that. I think inner_core is a test where we use a very under-resolved mesh so that it runs fast enough.

[bangerth]

I don't understand. Of course |current_constraints| might be changing and we haven't specified if it is allowed to change the type of constraint over time or only the values. But I assume that |current_constraints| computed at t=0 is more correct to use than |constraints| for |make_sparsity_pattern|.

Yes, we agree.

This patch computes |current_constraints| before using them.

But from what? The patch calls it from setup_dofs(), which is significantly before we even start a time step. It makes me uncomfortable to call compute_current_constraints() before we know that everything is set up for the first time step. I'm not saying that something bad actually happens -- it just seems logically wrong to do this. I'm going to add that I think the current place where we call this is also questionable: We call compute_current_constraints() *before* we call the update() function on all of the plugins. I think it should actually be after.

I think that this change not only allows us to add additional constraints (like the pressure dofs), but will also improve memory usage and performance, because we carry around plenty of unused sparsity pattern entries for boundary conditions that are only in |current_constraints|.

That, however, is an interesting point. Have you measured what difference it makes?

[tjhei]

But from what? The patch calls it from setup_dofs(), which is significantly
before we even start a time step.

I agree, this is a problem. We could delay constructing the sparsity patterns until we need them the first time. Then this setup would actually happen at t=0 and everything is set up correctly by then. I imagine this could be done by inspecting a flag right before we do the assembly. This would allow us to change the type of constraints during a computation too if we wanted to do that in the future. I can prepare a patch if you agree that this is acceptable.

[tjhei]

but will also improve memory usage and performance
That, however, is an interesting point. Have you measured what difference it
makes?

Of course not. ;-)
I don't think it makes a big difference, especially because the number of boundary dofs is quite small in a large computation. For small problems it is signifcant (not that performance matters here), see /tests/matrix_nonzeros_3 for example:

-Total system matrix nnz: 1,773
+Total system matrix nnz: 1,613
 system matrix nnz by block: 
          288         112           0           0           0
          112           0           0           0           0
-           0           0         289           0           0
+           0           0         129           0           0
            0           0           0         972           0
            0           0           0           0           0

[maxrudolph]

@tjhei @bangerth Do I understand correctly from this discussion that Timo will resolve this issue through a new pull request? If so, I could close this one.

[tjhei]

Do I understand correctly from this discussion that Timo will resolve this issue through a new pull request?

Yes, I will prepare something but leave this open for now. This way you guys can remind me in case I don't get it done soon.

[tjhei]

addressed in #1875