[WIP] Fastscape c++ integration with adapter

[Minerallo]

Description

This PR is a follow-up to #5323.
Since development has moved to a new main branch, this new PR continues the work with a rebased and updated version of the FastScape coupling plugin.

The key update in this version is the use of the FastScapeAdapter interface by @benbovy and @MFraters , which provides a consistent deal.II mesh-based surface coupling for both Box and SphericalShell geometries.

The current objective is to restore functionality for the Box geometry, now using the adapter interface. SphericalShell support will be integrated next.

---

Checklist

Before your first pull request:

• [ ] I have read the guidelines in our CONTRIBUTING.md document.

For all pull requests:

• [ ] I have followed the instructions for indenting my code.

For new features/models or changes of existing features:

• [ ] I have tested my new feature locally to ensure it is correct.
• [ ] I have created a testcase for the new feature/benchmark in the tests/ directory.
• [ ] I have added a changelog entry in the doc/modules/changes directory that will inform other users of my change.

[bangerth]

@Minerallo I've been reading through some parts of the code and found that there are a number of member variables that are either not used, or only ever set. Specifically, these are:

• [ ] fs_seed: never used
• [ ] p: I could not find where this variable is read from the input file. The variable name is too short to tell whether it's used anywhere :-)
• [ ] node_tolerance: I don't think this is used anywhere
• [ ] precision: This variable is set, but I don't think it's ever used.
• [ ] restart: Never used.
• [ ] end_time: This variable is read from the input file, but then never used. If you plan on using it, my preference would still be to get the end time from the SimulatorAccess object, rather than duplicate it in this class.

It's possible that you have plans for some of these variables. But if you don't, would you mind pulling the newest version of the branch (including my commits) and remove these variables? Feel free to check off the checkboxes above once you've done that.

[Minerallo]

@Minerallo I've been reading through some parts of the code and found that there are a number of member variables that are either not used, or only ever set. Specifically, these are:

• [x] fs_seed: never used
• [x] p: I could not find where this variable is read from the input file. The variable name is too short to tell whether it's used anywhere :-)
• [x] node_tolerance: I don't think this is used anywhere
• [x] precision: This variable is set, but I don't think it's ever used.
• [x] restart: Never used.
• [x] end_time: This variable is read from the input file, but then never used. If you plan on using it, my preference would still be to get the end time from the SimulatorAccess object, rather than duplicate it in this class. okay

It's possible that you have plans for some of these variables. But if you don't, would you mind pulling the newest version of the branch (including my commits) and remove these variables? Feel free to check off the checkboxes above once you've done that.

@bangerth I have removed all the variables for now, and will add them back if we need them later. p is the slop exponent for multidirectonal flow, this is computed internally in fastscapelib c++ now, right @benbovy ?

[bangerth]

There is an issue with n_grid_nodes, as just discussed.

For the box model, this variable isn't set at all. For the spherical model, you do

          n_grid_nodes = spherical_vertex_index_map.size();

but then (for both cases), in initialize(), you do

      n_grid_nodes = surface_mesh.n_active_cells();

which seems wrong simply based on the words: The left hand side indicates a number of nodes, the right hand side a number of cells.

Perhaps this kind of patch would seem right? I don't know what the rest of the code in this file does yet, so I'll let you figure out whether this makes sense:

--- a/source/mesh_deformation/fastscapecc.cc
+++ b/source/mesh_deformation/fastscapecc.cc
@@ -59,8 +59,6 @@ namespace aspect
       this->get_pcout() << "Geometry model type: " << typeid(geom_model).name() << std::endl;
 
       init_surface_mesh(geom_model);
-
-      n_grid_nodes = surface_mesh.n_active_cells();
     }
   
     template <int dim>
@@ -136,11 +134,13 @@ namespace aspect
                   }
               }
 
-          n_grid_nodes = spherical_vertex_index_map.size();
+          AssertDimension (spherical_vertex_index_map.size(), surface_mesh.n_used_vertices());
         }
       else
         AssertThrow(false, ExcMessage("FastScapecc plugin only supports Box or Spherical Shell geometries."));
 
+      n_grid_nodes = surface_mesh.n_used_vertices();
+
 
       // Having create the mesh, now set up the DoFHandlers and constraints objects
       surface_mesh_dof_handler.reinit(surface_mesh);

If that does make sense, please apply the patch for you locally.

[Minerallo]

yes I would think that n_used_vertices() make more sense to use to set the vectors using the patch.

Quantity	Value
surface_mesh.n_active_cells()	100
surface_mesh.n_used_vertices()	120

[benbovy]

which seems wrong simply based on the words: The left hand side indicates a number of nodes, the right hand side a number of cells.

FWIW, a Fastscapelib grid has no concept of a cell. Fastscapelib's algorithms for flow routing and erosion only need a grid (or a graph) of nodes connected to their direct neighbors. What those algorithms also need is an area surrounding each node, which may correspond to the area of a cell around the node (e.g., cell center) although the Fastscapelib doesn't need to know the exact geometry of the cell (which is expensive to store). There's some explanation here too: https://fastscapelib.readthedocs.io/en/latest/guide_grids.html#grid-representation

Something that we discussed during our last ASPECT/Fastscapelib coupling sprint in Potsdam: a Fastscapelib grid node may correspond to the center (centroid) of an ASPECT (deal.II) mesh cell, a node's neighbors may correspond to all adjacent cells and a node's "area" may correspond to the cell area. I'm not sure to remember how data fields in ASPECT are mapped to cell vertices (i.e., nodes) vs. cells in a deal.II mesh, but this is important to take into account when mapping ASPECT data onto a Fastscapelib grid.

[bangerth]

I see. I'll let the two of you figure out what is the appropriate approach. The current state does not seem useful since the variable is first set for the spherical shell grid in init_surface_mesh() (but not for the box geometry), and then immediately overwritten at the end of initialize(). I don't understand the code well enough to say how to fix this and will gladly let you figure this out :-)

[benbovy]

In #6204 the n_grid_nodes variable was set only once in initialize() after calling init_surface_mesh() for the appropriate geometry.

https://github.com/benbovy/aspect/blob/0a1a7f951bd878667d67b7fd0ed9996a91a58bde/source/mesh_deformation/fastscapecc.cc#L51-L101

I'm not sure why it has been changed since then.

I can try helping to figure out how to fix any issue related to n_grid_nodes (as it was implemented in #6204). However I might not be familiar enough with ASPECT to figure out what exactly is going on and what is needed, unfortunately.

[bangerth]

OK, that suggests that

n_grid_nodes = surface_mesh.n_active_cells();

is indeed the right piece of code, and the presence of

n_grid_nodes = spherical_vertex_index_map.size();

is wrong? @Minerallo I think I'll leave this to you then!

[Minerallo]

If we set n_grid_nodes = surface_mesh.n_active_cells(); during initialization, it should work for both Box and SphericalShell geometries. That means I might be able to remove the line:

n_grid_nodes = spherical_vertex_index_map.size();

However, I think getting the Box geometry work, as simpler case will help me to better understand the logic.

[benbovy]

p is the slop exponent for multidirectonal flow, this is computed internally in fastscapelib c++ now, right @benbovy ?

@Minerallo in Fastscapelib C++ this is a parameter of the multi_flow_router flow operator. You can find some examples of usage here: https://fastscapelib.readthedocs.io/en/latest/guide_flow.html#flow-routing-strategies.

[bangerth]

@Minerallo FYI, I implemented the use of deal.II's send/receive functions. This should make it substantially easier to switch to the map.