external mesh deformation base class

[tjhei]

FYI @bangerth

[bangerth]

we need compute_initial_deformation_on_boundary(). we should get the initial deformation from landlab

[bangerth]

need to refactor code in base class where we build a vector of initial deformations, so that we query by node index, not by node point.

[tjhei]

@bangerth Thank you. This is now ready for a first review.

We will need to switch to 9.6 for the tests to pass, see geodynamics#6638

Additionally, I have not done the refactor for the initial conditions so far, but I think this could be merged without.

[tjhei]

@bangerth ping

[tjhei]

@bangerth Ping.

The tests now pass.

[bangerth]

Here's the PointDataOut class. I'm going to have to go teach, but will be back at this for the rest when back.

[bangerth]

Do we eventually want to make this so that in parallel, we don't re-use ids?

[tjhei]

That would make sense, but the class has currently no knowledge about the MPI communicator to use.

[bangerth]

Yes, perhaps an issue for the future.

[bangerth]

Here's the rest. I apologize for taking so long to actually get to this! :-(

[bangerth]

Yes. We use this vector to evaluate at all boundary DoFs of the current process because this is what DoFTools::extract_boundary_dofs does:

  /**
   * Extract all degrees of freedom which are at the boundary and belong to
   * specified components of the solution. The function returns its results in
   * the form of an IndexSet that contains those entries that correspond to
   * these selected degrees of freedom, i.e., which are at the boundary and
   * belong to one of the selected components.
   *
   * By specifying the @p boundary_ids variable, you can select which boundary
   * indicators the faces have to have on which the degrees of freedom are
   * located that shall be extracted. If it is an empty list (the default), then
   * all boundary indicators are accepted.
   *
   * This function is used in step-11 and step-15, for example.
   *
   * @note If the DoFHandler object is defined on a
   * parallel Triangulation object, then the computed index set
   * will contain only those degrees of freedom on the boundary that belong to
   * the locally relevant set (see                                                    **********************************
   * @ref GlossLocallyRelevantDof "locally relevant DoFs"), i.e., the function
   * only considers faces of locally owned and ghost cells, but not of
   * artificial cells.
   * [...]
   */
  template <int dim, int spacedim>
  IndexSet
  extract_boundary_dofs(const DoFHandler<dim, spacedim>    &dof_handler,
                        const ComponentMask                &component_mask = {},
                        const std::set<types::boundary_id> &boundary_ids = {});

[bangerth]

You may want to adjust this part of the comment. We are already consistent on the locally relevant DoFs. The call to make_consistent() is probably not wrong, but I think also not necessary.

[bangerth]

Your choice :-)

[tjhei]

Thank you @bangerth. I replied to some of your comments and will address everything else in the next few days.

[bangerth]

Tell me if you want me to take another look!

[tjhei]

I am stuck with one problem prompted by your remarks:

You may want to adjust this part of the comment. We are already consistent on the locally relevant DoFs.

I don't think this is true: On each rank we receive data in all evaluation points that are contained within any of our locally owned cells. We then proceed to set each surface DoF of the ASPECT mesh to the value of the closest point and store this in v_interpolated. On each rank we therefore touch all locally active DoFs (all DoFs in our own cells) and write the velocity. This means on a processor boundary the two ranks likely disagree and there is no way to communicate which value is correct (the closer evaluation point should win!).

This means:

1. It is currently not consistent
2. If we make it consistent with make_consistent_in_parallel, the owner likely wins, which is often incorrect.

I don't see an easy way to do this correctly, because it would involve communicating values and distance to evaluation point for each DoF on a processor interface and then picking the best value (probably by the DoF owner followed by a broadcast to the others).

Sigh.

[tjhei]

That said, I would also be okay to merge what we have here and incrementally improve. The tests are passing...

[bangerth]

I am stuck with one problem prompted by your remarks:

You may want to adjust this part of the comment. We are already consistent on the locally relevant DoFs.

I don't think this is true: On each rank we receive data in all evaluation points that are contained within any of our locally owned cells. We then proceed to set each surface DoF of the ASPECT mesh to the value of the closest point and store this in v_interpolated. On each rank we therefore touch all locally active DoFs (all DoFs in our own cells) and write the velocity. This means on a processor boundary the two ranks likely disagree and there is no way to communicate which value is correct (the closer evaluation point should win!).

Oh yes, I get it now. You're computing this also going for ghost nodes, which I thought was enough to be consistent. But you may be computing the wrong thing. So yes, you need to call the make_consistent() function.

This means:

1. It is currently not consistent
2. If we make it consistent with make_consistent_in_parallel, the owner likely wins, which is often incorrect.

I don't see an easy way to do this correctly, because it would involve communicating values and distance to evaluation point for each DoF on a processor interface and then picking the best value (probably by the DoF owner followed by a broadcast to the others).

Hm, good point. And if two processes write constraints into locally active DoFs that are different, then make_consistent() likely trips up (or at least should).

At first I was thinking that you need to make sure that you only ever create constraints for locally owned DoFs. That avoids the conflicts, but you're right that that doesn't mean that these are correct. In particular, you'd end up with a situation where you'd get different results depending on the number of processes.

Do you need to split the loop you have into two? Perhaps something like this (I think this is what you are proposing above):

• Loop over all points and DoFs like you currently have, and for each DoF (locally owned or not) record the closest eval point.
• Exchange this information by taking the minimum for each entry in the vector over all processes. (No idea how to actually achieve that, but let's stick with a conceptual idea.)
• Do another loop like above. If the distance between a point and a DoF is equal to the minimum computed in the second step, then do create a constraint. (Might of course run into problems when there are two eval points at exactly the same distance from a DoF.)
• Call make_consistent().

[bangerth]

That said, I would also be okay to merge what we have here and incrementally improve. The tests are passing...

Let's see if we can come up with a solution over the next couple of days. If not, we'll just go with what you have for now...

[tjhei]

Perhaps something like this (I think this is what you are proposing above)

Yes, this is what I was hinting at. But, as you said, it would involve a custom communication step.

• Might of course run into problems when there are two eval points at exactly the same distance from a DoF.

You make an important point. I think this situation is not only theoretical but very likely to happen considering one use case are structured grids that do not line up.

I think the only reasonable option is to (somehow) precompute a globally injective mapping from evaluation point index to DoF. The owner of the DoF can play the role of the tie breaker for points with equal distance. With this mapping, each processors loops over their list and writes their values into a global vector. After a call to compress(), every rank has consistent information and can create their AffineConstraints.

I would construct the list in the following way:

1. Build a vector of (DoF index, distance, evaluation point index, rank) on each rank (similar to what I have implemented but keep track of the index instead of copying the value)
2. Send each tuple to the owner of that DoF
3. Receive a list of tuples for my own DoFs, merge with my own list (keeping the one with smaller distance)
4. Send the new list of tuples to the respective owner of the evaluation point (this is the rank stored in the tuple)
5. Compile a final list on each rank (own list plus anything received in 4)
6. Create a map from evaluation point index to DoF (some evaluation points might not appear in this map and are therefore not used!)

With this map in place it is as simple as iterating over all evaluation points, looking them up in the map, and setting the values.

This sounds like a custom MPI routine where I can not use any linear algebra or other high-level deal.II functionality. Luckily, this is all point to point communication.

Any suggestions for simplification?

[tjhei]

ping @bangerth

[bangerth]

Yes, tough problem. You're right that that sounds like a function best suited to namespace VectorTools or namespace Utilities::MPI. Perhaps into the latter, with an interface of the form

   std::vector< [rank,node_point index on that rank] >        // one entry for each evaluation point
   find_closest_node_points (const std::vector<Point<spacedim>> &locally_owned_node_points,
                                                   const std::vector<Point<spacedim>> &locally_owned_evaluation_points)

I think that encapsulates the interface you need. The steps you mention seem reasonable to me. Perhaps you need to pass more information in than just the locally owned node points (the locally relevant ones? but then you also have to express ownership).

I did think that at one point you have to perform some kind of min-operations on the elements of a vector of distances (or a vector of (distance, rank) tuples) to figure out who has the closest node. (If you store 1/distance, then it just becomes a max operation, and if you need a tie breaker, that involves the rank part of the pair.) You can implement that as a custom functor with Utilities::MPI::all_reduce(), assuming that you can store all answers on each process.

[bangerth]

Do you think you can break this problem into smaller chunks that you could independently implement and test? I think that would make it much easier to review and convince ourselves that that's the way to go.

[tjhei]

Yeah, sounds good. I think we can implement these things independently of the mesh deformation class (the code in this PR) and unit test it.
So:

1. We can merge this PR now if you agree and update it to use the new functionality later.
2. The function that sets up the communication pattern as discussed can be unit tested and merged separately.
3. They can be combined as a follow-up.

I think we should do all of this in this repo even though 2 could go into deal.II eventually.