T-Model Configuration

[danieldouglas92]

This a first pass at testing out the "T-Model" configuration, with a 2D ASPECT model and a 2D LandLab model. The test model that is included in this PR does not average the LandLab topography along the z-axis, though this could be easily achieved within the update_until() function in the python file. This test shows that when the resolution of the ASPECT/LandLab meshes are the same on the surface, the topography in LandLab and ASPECT exactly align where the two models intersect (see figure below). The magenta line is the topography output by LandLab, and the topography exactly matches the displacement of the ASPECT mesh.

The only change that I had to make to ASPECT code to get this to work was in the definition of the surface points in landlab.cc. In 2D, the surface points (for any given y-value on the LandLab mesh) was always set equal to the x-value of the LandLab mesh, and the surface value of the ASPECT mesh. When determining the mapping in set_evalution_points, this meant that points at y=100 km and y= 0 km on the LandLab mesh were considered equally close to the ASPECT surface.

As with the other PR I opened, I'm sure there is a more robust way to do this and it'll be interesting to see how this works for unstructured grids, but I got around this issue by just setting the y-value of the surface points to be infinity if they were far enough away from the ASPECT surface. Therefore, if the LandLab mesh extends from -Y to + Y, and the ASPECT surface is at y=0, then only the points on the LandLab mesh near y=0 will be in direct communication with ASPECT. However, because this requires a tolerance, I doubt that this will be as robust for unstructured grids. Still, below is what the model looks like without the changes in this PR:

Looking forward to hearing what you think @bangerth @tjhei!

[tjhei]

As with the other PR I opened, I'm sure there is a more robust way to do this and it'll be interesting to see how this works for unstructured grids, but I got around this issue by just setting the y-value of the surface points to be infinity if they were far enough away from the ASPECT surface.

Interesting. And good to know that this way things are already close to working! :-)

I am not sure about this approach though: You now have a bunch of invalid evaluation points "high in the sky" that don't contribute anything (we are unable to evaluate them, because they are outside the mesh and we can't interpolate velocities to it) and that are basically unused. I am surprised this doesn't crash.

My approach would be the following:

1. In the python script, take the 2d grid points from Landlab and create a set evaluation points at the y=0 plane (one for each row). Send those to ASPECT.
2. When returning elevation/velocities from Landlab, average and return the value in the corresponding evaluation point.
3. When ASPECT gives a velocity, duplicate in that row and give that information to Landlab.

[danieldouglas92]

@tjhei Yes that definitely makes more sense! I've quickly put something together that yields the same model result, but instead passes a 1D line of evaluation points between landlab/aspect instead of a 2D plane when the ASPECT dim==2

[danieldouglas92]

@tjhei I updated this PR to remove everything except for the example model since all of the code required for this was already merged in other PR's.

[tjhei]

Can you fix indentation? Thanks.