[WIP] Fastscape c++ integration with adapter #6389
Conversation
Fastscapelib and Xtensor are treated as external dependencies (using find_package in the CMake configuration).
| = | ||
| // fastscapelib::make_spl_eroder(flow_graph, kff, n, m, 1e-5) | ||
| fastscapelib::make_spl_eroder(flow_graph, kff, n, m, 1e-5); | ||
| fastscapelib::make_spl_eroder(flow_graph, 1e-6, 0.4, 1, 1e-5); |
There was a problem hiding this comment.
The values make more sense than the variables, m (area exponent) and n (slope exponent) were flipped. The right order is:
fastscapelib::make_spl_eroder(flow_graph, kff, m, n, 1e-5);(side note: this is also why I prefer the names area_exp / slope_exp over m / n to avoid confusion and one character variable names).
| std::cout << "KFF : " << kff<< std::endl; | ||
|
|
||
| //Only for raster grid | ||
| //To propose if we add back the raster grid option later |
There was a problem hiding this comment.
Fastscapelib's diffusion is indeed only for the fixed-resolution raster grid since this is a requirement of the ADI numerical scheme used to efficiently solve diffusion.
If I understood correctly, there's an aspect diffusion plugin in the works? Or maybe already working for box geometry? I guess that it will be easier to reuse that plugin on the deal.II grid used by Fastscape instead of introducing a special case here?
There was a problem hiding this comment.
ASPECT has a surface diffuser in a different mesh deformation plugin. I'm not entirely sure at the moment whether one can have more than one surface deformation plugin, though.
There was a problem hiding this comment.
Separately, though: ASPECT's diffuser works on the ASPECT mesh, but that is for most simulations going to be substantially coarser than FastScape++'s mesh. As a consequence, I think at least in the long run it would be useful to have surface diffusion work within FastScape++.
| const double surface_cell_area = surface_mesh.begin_active()->measure(); | ||
|
|
||
| // 1. Create the FastScape grid adapter | ||
| // TODO: I haven't looked at the GridAdapter yet, but we should see whether we can't get rid of the const_cast. |
There was a problem hiding this comment.
I think it should work without the const_cast? The Fastscapelib grid adapter class internally keeps a non-const reference to the underlying mesh (although the adapter class shouldn't change any state of the underlying mesh).
| // auto diffusion_eroder = fastscapelib::make_diffusion_adi_eroder(grid, kdd); | ||
|
|
||
| // uplift rate in meters per year for Fastscape | ||
| std::vector<double> uplift_rate_in_m_year(vz.size()); |
There was a problem hiding this comment.
| std::vector<double> uplift_rate_in_m_year(vz.size()); | |
| xt::xarray<double> uplift_rate_in_m_year({ vz.size() }); |
Create directly an xtensor container instead of creating an std::vector and later use xt::adapt to make it xtensor-compatible. Unless this variable is later used by some function that specifically requires an std::vector (note that if those functions require STL-compatible they should probably work with xtensor containers as well, maybe with little adjustement). See, e.g., https://xtensor.readthedocs.io/en/latest/quickref/basic.html.
This could probably be done with other variables as well.
| } | ||
| // auto elevation = xt::adapt(h, shape); | ||
| auto elevation_old = xt::adapt(h_old, shape); | ||
| auto uplift_rate = xt::adapt(uplift_rate_in_m_year, shape); |
There was a problem hiding this comment.
another variable as xt::adapt(...) may not be needed (see https://github.com/geodynamics/aspect/pull/6389/files#r2154251150).
| auto row_bounds = xt::view(uplift_rate, xt::keep(0, -1), xt::all()); | ||
| row_bounds = 0.0; | ||
| auto col_bounds = xt::view(uplift_rate, xt::all(), xt::keep(0, -1)); | ||
| col_bounds = 0.0; |
There was a problem hiding this comment.
This is likely wrong. It tries to create views of uplift_rate as a 2-dimensional array, but uplift_rate is actually a 1-dimensional array with size = n_grid_nodes (i.e., the grid adapter's size).
I guess that the goal is to set fixed "boundary" conditions? This can be done in a grid-agnostic way by iterating over the flow graph's base level nodes, e.g.,
for (auto& bi : flow_graph.base_levels())
{
uplift_rate[bi] = 0.0;
}(I haven't run and checked the code above)
| // auto diffusion_eroder = fastscapelib::make_diffusion_adi_eroder(grid, kdd); | ||
|
|
||
| // Create an xarray with correct shape_vec | ||
| xt::xarray<double> uplift_rate_vec = xt::zeros<double>({vz.size()}); |
There was a problem hiding this comment.
Would be safe to have an assertion check that vz.size() equals grid.size() (I guess it should be equal).
| uplift_rate_vec[i] = vz[i] / year_in_seconds; | ||
|
|
||
| // Reshape uplift_rate to match grid shape (same as elevation) | ||
| auto uplift_rate = xt::adapt(uplift_rate_vec.data(), shape); |
There was a problem hiding this comment.
Unless I miss something, this variable uplift_rate is not needed. uplift_rate_vec should already have the right shape (1-dimensional array with size equal to grid size) and may be reused directly below.
Same for other variables auto var_name = xt::adapt(...) (like elevation and elevation_old below): they are likely not needed either. Instead you can create xt::xarray containers directly and assign the data elements to them like if they were std::vector containers.
| // TODO: You are using these three arrays in the following way: | ||
| // temp[0][i] = topography value at some location | ||
| // temp[1][i] = DoF index within surface_dof_handler for this location | ||
| // temp[2][i] = surface_uplift_rate for this location | ||
| // I think this would be more transparent if you had the following: | ||
| // struct ValuesAtSurfaceVertex | ||
| // { | ||
| // double topography; | ||
| // double surface_uplift_rate; | ||
| // }; | ||
| // and then you store everything in a | ||
| // std::multimap<types::global_dof_index, ValuesAtSurfaceVertex> surface_vertex_to_surface_values; | ||
| // Note that here I chose a multimap because you visit vertices more than once from all of | ||
| // the adjacent cells, at least in the current code you enter the same vertex into | ||
| //the temp vectors multiple times. I *think* that you should get the same values every | ||
| //time you visit a vertex, so perhaps a std::map is a better choice. | ||
| std::vector<std::vector<double>> temporary_variables(3, std::vector<double>()); |
There was a problem hiding this comment.
I agree using a struct / map would be clearer.
I'm struggling to follow the logic here. I guess that the goal of using temporary variables is for gathering / merging together data from the MPI processes, right? In the case of a single process, the data is simply copied through the temporary variables... Is that correct?
For better clarity if would be nice if this logic could be refactored into a separate function that returns new (xtensor) array containers for topography and surface uplift rate. If it is too expensive to create and re-allocate new array containers each time, the function could instead take topography and surface uplift rate arrays as input, maybe resize them and update their values (same for the temporary structures if needed). IIUC we could even try to optimize the case of a single MPI process (i.e., bypass temporary variables).
There was a problem hiding this comment.
@benbovy I think you're interpreting the reasons behind this code right (but then, of course, I didn't write it and so don't know). When you do plain MPI, working with regular arrays such as those here makes sense. But in ASPECT, we can use deal.II functions such as Utilities::MPI::isend/ireceive that can deal with arbitrary data structures and so sending a map is no more work than sending a vector. I think that will make the code easier to read.
In the long run, I want to separate out the code that evaluates the ASPECT solution at a set of points into a separate function, because we also need that when interfacing with other surface codes. But that's for a later step.
|
@Minerallo FYI, I implemented the use of deal.II's send/receive functions. This should make it substantially easier to switch to the map. |
…copy-patch1 Ensure we use the vertical direction, not the z-direction.
Minerallo
changed the title
4 participants
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
FastScapeAdapterinterface 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:
For all pull requests:
For new features/models or changes of existing features: