CG: distribute periodicity and hanging node constraints

[richardschu]

fixes #650
For the structure module, similar problems might appear. Maybe no one solved structure problems with hanging nodes or periodic constraints yet? There, it would actually impact the performance of the nonlinear solver I believe, as we distribute constraints on the linearization vector. I did not test the changes in the structure module yet, but I doubt any of these consider hanging nodes anyways. Will test manually and then update here. Until then, this is on do not merge.

[richardschu]

I checked now the parts for the linear and nonlinear structure solvers. The same issue was also appearing there. With the lines I added, these are fixed (some pictures below ...). Also, for the linear solver, the iteration count went from 28.3 to 12.8 on average over 200 time steps. For the nonlinear solver, I could not really see a difference (probably the test case was not really challenging, only 1.1 Newton iterations on average over 10 load steps).
unsteady linear solver before fix:

unsteady linear solver after fix:

quasistatic solver before fix:

quasistatic solver after fix:

[nfehn]

Can you explain why iteration counts are changing, because issue #650 mostly talks about the impact in terms of postprocessing/visualization?

[nfehn]

I think we should discuss again the design. When looking at the code changes, it might be that we have forgotten it somewhere. I think this originates from the fact that the new line of code is added at many places in the SpatialOperator but for example not in OperatorBase (or in a derived class of OperatorBase). It is not so long ago that we introduced two separate affine constraints objects and maybe we would have done it differently if we had been aware of this. If the spatial operator needs the new lines of code for correct evaluation of the PDE operator, why not let this operator care about imposing the constraints (single responsibility principle)? Similarly, the postprocessor might want to care about this so that we do not per se distribute whenever a solution has been computed. I see this command somewhat similar to zeroing ghost values for example. These lines of code appear in my opinion also closer to the matrix-free operator or closer to the postprocessor.

[nfehn]

Because of your comment "Until then, this is on do not merge.", I converted this PR to a draft.

[nfehn]

Maybe no one solved structure problems with hanging nodes or periodic constraints yet?

Yes, this is not unlikely. I wonder right now what periodic BCs would mean for structure problems? Wouldn't the system matrix be singular in this case so that the problem is not solvable? (edit: in a single direction it would work I guess, if the remaining BCs are consistent with the periodicity constraints)

[kronbichler]

I wonder right now what periodic BCs would mean for structure problems? Wouldn't the system matrix be singular in this case so that the problem is not solvable?

Isn't it the same as with the other singular systems we deal with (e.g. Poisson problem with pure periodic/Neumann conditions)? That means we have to make sure the system/rhs are compatible with the null space and in case of multigrid preserve the null space during the whole hierarchy and coarse solver?

[nfehn]

Most changes are made after solving. This change is made before solving. Why?

[nfehn]

Did you see the TODO a few lines below? This code that is currently commented, but code that we might want to take care of as well.

[richardschu]

Most changes are made after solving. This change is made before solving. Why?

see #651 (comment)
Here, the sol is then used in the Newton solver to update the linearization vector
https://github.com/richardschu/exadg/blob/c5f489f1085b78f951d53cc62f510028c74c38be/include/exadg/solvers_and_preconditioners/newton/newton_solver.h#L79
and evaluate the residual
https://github.com/richardschu/exadg/blob/c5f489f1085b78f951d53cc62f510028c74c38be/include/exadg/solvers_and_preconditioners/newton/newton_solver.h#L104

Otherwise we have the zero hanging node values in the linearization vector.
I was surprised to not see an impact on the nonlinear iteration count, but maybe the solution we are searching for (with zero values at hanging nodes) is still close enough to the "real" one with correct hanging node values.

[nfehn]

I agree that we need this vector for the evaluation of the residual/linearized operator. And it is also clear to me that we need to take these constraints into account correctly when evaluating the integrals. However, from a software design perspective, the problem is in my opinion that responsibilities are beginning to blur if we add these lines here.

[richardschu]

The Operator has all affine_constraints objects, so I think it should also enforce the constraints.

I would not put it in the postprocessor, since the idea (up to now) for the postprocessor was, that it is entirely optional (see throughput applications, where it is not even set up). If we move the lines, I think they have to go together with the
elasticity_operator_nonlinear.set_inhomogeneous_boundary_values(sol);
otherwise we split things up unintuitively(?).

Where is in your opinion the right place to enforce constraints on a vector?
The Operator has the affine_constraints, so it should do it, right?
The sol could already have the right BCs enforced when calling solve_nonlinear(sol).
So should we add a operator->enforce_constraints(), which is called before solve_nonlinear and solve_linear(sol) from the driver(s)? But that is also not really intuitive from the driver side: How would I expect that I have to update the vectors constraints due to internals of how the operator works? For the postprocessor, we could think about it, but then it would be different for linear and nonlinear solvers, which is also kinda confusing. So I would rather say, maybe introduce a operator->enforce_constraints() to bundle the constraints, but when calling operator->solve_nonlinear(sol) or operator->solve_linear(sol), this should be done in the operator side automatically where needed.

[richardschu]

Did you see the TODO a few lines below? This code that is currently commented, but code that we might want to take care of as well.

I think we should remove it, when coming up with something new here.

[richardschu]

Ok so I repeated the tests with the linear elasticity solver (due to #652 I cannot use the steady solver currently).
I changed my preconditioner settings to whatever I had yesterday, but adding line
https://github.com/richardschu/exadg/blob/c5f489f1085b78f951d53cc62f510028c74c38be/include/exadg/structure/spatial_discretization/operator.cpp#L1007
changes the iteration count from 23 on average to 9.35.
So since I could repeat this experiment, I think this is actually a thing. I would explain it via:

1. we solve some system and get a solution, the hanging node values are zero (rather than the interpolated value on the edge/face). Inhom Dirichlet BCs however, are enforced(!).
2. we hand that solution over to the time integrator. (For the the steady solver it might thus be a postprocessing issue similar to the Poisson solver!)
3. the time integrator keeps those zero values at the hanging nodes over time.
4. The displacement, velocity and acceleration values that enter the rhs have the zeros at the hanging nodes at all times, BUT the solution itself id "rougher" than the one without hanging nodes, since we have those zero "spikes" at the hanging nodes. I think the multigrid does not like those rough solutions (and righthandsides), because they cannot be resolved properly on coarser grids, right? So thats the point where I think the solver performance is influenced.
5. distributing the hanging node constraints gets rid of the zeroes, but the physics are actually wrong, aren't they? We enforced zero displacement/velocity/acceleration at the hanging nodes?
   ...that is at least my understanding for the linear solver only.

[richardschu]

Isn't it the same as with the other singular systems we deal with (e.g. Poisson problem with pure periodic/Neumann conditions)? That means we have to make sure the system/rhs are compatible with the null space and in case of multigrid preserve the null space during the whole hierarchy and coarse solver?

I would agree on that, but I have not tested it myself!

[richardschu]

I think we should discuss again the design. When looking at the code changes, it might be that we have forgotten it somewhere. I think this originates from the fact that the new line of code is added at many places in the SpatialOperator but for example not in OperatorBase (or in a derived class of OperatorBase). It is not so long ago that we introduced two separate affine constraints objects and maybe we would have done it differently if we had been aware of this. If the spatial operator needs the new lines of code for correct evaluation of the PDE operator, why not let this operator care about imposing the constraints (single responsibility principle)? Similarly, the postprocessor might want to care about this so that we do not per se distribute whenever a solution has been computed. I see this command somewhat similar to zeroing ghost values for example. These lines of code appear in my opinion also closer to the matrix-free operator or closer to the postprocessor.

My logic was to add it at places, where we call set_inhomogeneous_boundary_values(), so maybe we should allow that function to also distribute the hanging node constraints? However, from my understanding the hanging node constraints might be treated differently in some cases, since they are x1 = 0.5*(x2+x3) instead of x1 = 3, so how we treat the constraints in and around MatrixFree plays a role. I did not fully dig into this yet.

[nfehn]

@kronbichler agree. I should make my statement more precise. I do not want to argue that one could not deal with this technically, but wondered whether we could expect the problem for such a setup to be solvable with the present implementation apart from constraints.distribute() (I am not aware that we did anything related to "rhs compatible with the nullspace" for the structure module). So I was wondering how a reasonable problem setup with periodic BCs would look like for the Structure module. Have you something in mind, @kronbichler ? which problem did you think of or does one of your applications needs this, @richardschu?

[richardschu]

which problem did you think of or does one of your applications needs this, @richardschu?

No I have no application in mind that needs that. It is simply because that is the name of the object affine_constraints_periodicity_and_hanging_nodes 😄
In fact, I just played around with the hanging nodes to create tests on my end and stumbled upon this. I do not even use hanging nodes rn even though they are extremely cool 💯

[nfehn]

@richardschu thanks for your explanations. What I do not understand is why the operator evaluation is affected by this. During initialization, I thought that MatrixFree receives both AffineConstraints objects and then evaluates the constraints correctly in matrixfree->loop(). It might be that I had a wrong picture in mind of what is happening internally when implementing these things. So maybe @kronbichler could explain what we do wrong such that the operator evaluations are done incorrectly (or at least differently than what we expected before). Do I remember correctly that we decided for read_dof_values() instead of read_dof_values_plain()?

[kronbichler]

What I do not understand is why the operator evaluation is affected by this. During initialization, I thought that MatrixFree receives both AffineConstraints objects and then evaluates the constraints correctly in matrixfree->loop(). It might be that I had a wrong picture in mind of what is happening internally when implementing these things.

In general, the design internally is such that FEEvaluation::read_dof_values will resolve the constraints that are handed in to MatrixFree::reinit. So in case all data access is via this function (and does not mix read_dof_values_plain or dealii::DoFCellAccessor that accesses the data from the plain solution entries via cell->get_dof_indices, all data should be consistent with the given constraints. If we see different iteration counts, this indicates that the solution does not get appropriately interpolated to the given hanging node/periodicity constraints somewhere in the solver process. It could be the residual, the computation of the right hand side/body forces, the boundary conditions, or similar issues.

I made a small search

$ grep -n read_dof_values_plain ~/Work/dg/exadg/include/ -r
/home/martin/Work/dg/exadg/include/exadg/postprocessor/normal_flux_calculation.cpp:81:      integrator.read_dof_values_plain(solution);
/home/martin/Work/dg/exadg/include/exadg/structure/spatial_discretization/operator.h:430:  // While dealii::FEEvaluation::read_dof_values_plain() would take into account inhomogeneous
/home/martin/Work/dg/exadg/include/exadg/structure/spatial_discretization/operators/body_force_operator.cpp:83:      integrator.read_dof_values_plain(src);
/home/martin/Work/dg/exadg/include/exadg/fluid_structure_interaction/precice/quad_coupling.h:275:    phi.read_dof_values_plain(data_vector);

which shows all appearances of the read_dof_values_plain function. At least the body force operator looks suspicious, but I did not look more closely into it yet. @richardschu is this function used? Similarly, I searched for get_dof_indices to see where we used it, and it appears as if most uses apart from

exadg/include/exadg/operators/operator_base.cpp

Lines 2174 to 2180 in 0075db2

	if(is_mg)
	cell_v->get_mg_dof_indices(dof_indices);
	else
	cell_v->get_dof_indices(dof_indices);
	for(auto const & i : dof_indices)
	weights[i] += 1.;

are either used in conjunction with AffineConstraints object of the right type or in the postprocessor.

If the spatial operator needs the new lines of code for correct evaluation of the PDE operator, why not let this operator care about imposing the constraints (single responsibility principle)? Similarly, the postprocessor might want to care about this so that we do not per se distribute whenever a solution has been computed. I see this command somewhat similar to zeroing ghost values for example. These lines of code appear in my opinion also closer to the matrix-free operator or closer to the postprocessor.

The question I have here is at which level we would like the interpolation of constraints to happen. As I elaborated previously, one concern with AffineConstraints::distribute is that it could be expensive, so doing it near the solver might be problematic and defeating the ambition to do high-performance computing. At the same time, we already have the appropriate infrastructure in MatrixFree to do it more efficiently, so something like matrix_free->cell_loop() that internally does FEEvaluation::read_dof_values() (resolving constraints) and then calls FEEvaluation::set_dof_values_plain could be an appropriate design. It might also help us with some long-term plans I've had, namely eliminating those degrees of freedom from the vectors in linear/nonlinear solvers that are always constrained (to zero or other DoFs), thus reducing the size of the vectors. This vectors would only be compatible within their use in MatrixFree, but then we would need a way to convert them to the regular vectors. This might be in line with that goal, albeit I admit it is a bit ambitious.

[richardschu]

is this function used?
No, I did not use a rhs for these tests @kronbichler .

I think one solution would be to introduce

1. virtual void OperatorBase::enforce_hanging_nodes_and_periodicity_constraints(sol);(similar toset_inhomogeneous_boundary_values) where we do a matrix_free->cell_loop() using the right DoF index, where we have the affine_constraints_hanging_nodes_and_periodicity`, see https://github.com/exadg/exadg/blob/ff65a5fc9815b5aa9cd1d8768dc0de404d88581b/include/exadg/structure/spatial_discretization/operator.cpp#L144C1-L153C95
2. Overwrite in all CG solvers
   exadg/poisson/spatial_discretization/laplace_operator.h
exadg/structure/spatial_discretization/operators/elasticity_operator_base.h
   'exadg/structure/spatial_discretization/operators/mass_operator.h'
3. replace where we have now affine_constraints_hanging_nodes_and_periodicity.distribute(sol) with that new
   operator->enforce_hanging_nodes_and_periodicity_constraints(sol) to avoid constraints.distrubute(sol).

What do you think @nfehn?