What do the AffineConstraints::distribute_local_to_global(...) functions do on a matrix/vector level?

[Kyle Schwiebert]

Hello all,

Thank you for taking the time to look at my question. First, I'll ask a couple of basic questions about the built-in functions, and then I'll give a few details of why I ask.

Does the inhomogeneity-handling call distribute_local_to_global(...,false) do the following:

Say that we have a constraint x_i = a. I would guess that the function takes the following steps: (matrix is M, RHS is right hand side)

1) Set M(i,i) = L and set RHS(i) = 0, where L is somehow "nice" for the matrix's spectrum.

2) Set the rest of the entries of the ith row and column of M = 0. Perhaps technically we don't bother to zero the column since we plan to call AffineConstraints::set_zero.

3) For each DoF j which is coupled to DoF i and is not itself constrained, we do RHS(j) -= a*M(j,i).

Now considering the same constraint with instead a call to distribute_local_to_global(local_rhs,dof_inds,global_rhs) I would guess that we just set RHS(i) = 0 and do nothing else special.

Most importantly, is the above correct?

If so, I am having trouble understanding where I might be going wrong. I'm solving a time dependent coupled problem with a decoupled solver. By decoupled solver I mean that information from the second domain only appears in the right hand side. To make matters one step more complex, one term in my matrix is time-dependent and so I do not recompute the whole matrix. My procedure is this:

1) Assemble the static part of the matrix with a zero right hand side. This matrix will never be touched again. The vector, RHS_bc, which will also never be touched again contains no nonzero entries except for those created by the call to distribute_local_to_global(...,false).

2) I then assemble the time-dependent terms into a separate matrix and vector, RHS_td, using distribute_local_to_global(...,false) and add the RHS_bc to RHS_td.

3) I then assemble the coupled terms from the other solver into RHS_c with distribute_local_to_global(local_rhs,dof_inds,RHS_c).

4) I then add RHS_c to RHS_td and solve (being careful to call constraints.set_zero(solution) before the solve and constraints.distribute(solution) after the solve). I am of course using a custom class similar to the linear operator classes to bundle the vmult() outputs of the static and time-dependent matrices.

However, I am seeing blowup in my solution and I'm seeing some indication that the inhomogeneous boundary conditions are to blame. Note that I also tried calling constraints.set_zero(RHS_c) in step 4 above and as I expected that had no change (since presumably constrained entries are already zero).

Thank you again for taking the time to engage with my question and in advance for any tips you are able to offer.

-Kyle

[Luca Heltai]

Kyle,

what you are asking is a very difficult question. If you take a look at the (new) step-86 branch (https://github.com/dealii/dealii/pull/15540), you’ll find a similar approach to what you describe.

Everything you said is correct, with possibly one problem (that I have encountered very recently, and may be related to what you are observing).

Boundary conditions and hanging nodes have a very delicate interplay. In particular, when you *build* the constraints, if one of the hanging nodes (dependent) has one of its dependencies on the boundary, then when you add non-homogeneous boundary conditions to the constraints, the hanging node constraints are also modified to take into account the values of the boundary conditions. This computation is done *when you compute the constraints for BC*, and it holds the same values for the rest of the life of the affine constraints. When you `distribute`, the constraints are computed with the BC that you had at the beginning.

This information is retained even if you change later your (time dependent) boundary conditions, and may result in the hanging nodes being compute in a wrong manner (in your solution, you may see discontinuous solutions instead of continuous ones where hanging nodes interact with boundary nodes).

This affects the final call to “distribute”, not the “distribute_local_to_global” (with false as last argument), since that call condensates away the constraints. So your process is correct, but `distribute` does not do the right job, because the constraints on the hanging nodes have not been updated to take into account the new boundary conditions.

To verify that this is the case, try doing the following when you update the boundary conditions: clear the constraints, and recreate them from scratch, and see if this solves your problem. The `distribute_local_to_global` will behave exactly in the same way, but `distribute` should now do the right thing.

L.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/0589f784-e1ed-4699-bf30-a7358b3acee1n%40googlegroups.com.

[Wolfgang Bangerth]

On 7/21/23 04:48, Luca Heltai wrote:
> Boundary conditions and hanging nodes have a very delicate interplay.

As just another data point, step-26 re-builds the matrix in every time step
just so that we can apply boundary conditions and hanging node constraints
from scratch.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Kyle Schwiebert]

Luca,

Thank you very much for the reply. That is all very very helpful.

I see what you are saying about the hanging node constraints. I am using a mesh from gmsh directly (without any sort of refinement), which I am pretty sure is conforming. My code does include a call to make_hanging_node_constraints just to be safe, but commenting that out does not seem to cause any difference in behavior.

It is very helpful that you mention time-dependent constraints, which I will certainly need to implement sooner rather than later (but the BCs in the current problem are actually constant). I actually believe the scheme I listed above does not totally work in that case. If we have x_i = a(t), then RHS_bc(i) must change from having the contribution -m(j,i)a(t) -> -m(j,i)a(t + dt). Unfortunately, I do not see a way to do that other than to rebuild that static part of the matrix in order to call distribute_local_to_global(...,false) again. Something like multiplying each entry by a factor does not work since some entries in RHS_bc could contain information from more than one constrained node.

Perhaps then I should look into using distribute_local_to_global(...,true) and let the time_dependent assembly handle the boundary conditions.

I assume that the call to distribute_local_to_global(...,true) instead does:

1) Sets M(i,i) = L and RHS(i) = L*a.

2) Zeros out row i except for the diagonal entry.

The extra trick here is that since x_i is not being set to zero prior to the solve step, if you are reusing the preconditioner across time steps, you need to control L so it is always the exact same, but that is not too difficult in this easy situation. It is difficult if one is using periodic or other constraints which constrain one DoF to others (rather than to a constant) there does not seem to be an easy way to get access to which DoFs those are for a given constrained DoF. This would inhibit efficient manual setting of the scaling of the rows. The current code does not have time-dependent or periodic BCs, so these considerations are perhaps for a later date for me to consider.

I found a major issue in an unrelated part of my code, so it is mostly working now, but I still have a feeling there is still at least a minor issue with how I am handling the BCs.

Regards,

Kyle

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---

You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/B5IJB1O-BSw/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/03b57226-f66e-011a-3adb-58cdf995d661%40colostate.edu.