[WIP] Add topograhpy infrastructure

[MFraters]

An alternative implementation to #814 for the topography interface (issue #780). The difference is that this implementation can have different topography models per feature. It is still very much a work in progress, since not all features and geometries are implemented, but is as already workable, as can be seen in the example below. @Djneu would this also work for your use-cases? @alanjyu where you also interested in this, and if so, does this work for you as well?

I am thining of turning the topography input from a single value into a surface, like the min and max depht. That would allow to seamlessly use Litho1.0 data. (as a side-note, I probably want to rework the current implementation of surfaces to be more flexible and in line with the rest, but that is for a different pull request.)

This is made with the following world builder file:

{
  "version":"1.1",
  "coordinate system":{"model":"spherical", "depth method":"begin at end segment"},
    "surface temperature":293.15,
  "force surface temperature":true,
  "features":
  [
    {"model":"continental plate", "name":"Continent plate A", "coordinates":[[-1,10],[41,10],[41,15],[-1,15]], "min depth":-1e2,
     "topography models":[{"model":"uniform","topography":50000}],
     "temperature models":[{"model":"uniform", "min depth":-1e2, "max depth":95e3, "temperature":500}],
     "composition models":[{"model":"uniform", "compositions":[0], "min depth":-1e2}]},

    {"model":"continental plate", "name":"Continent plate B", "coordinates":[[20,12],[31,12],[31,15],[20,15]], "min depth":-1e2,
     "topography models":[{"model":"uniform","topography":70000}],
     "temperature models":[{"model":"uniform", "min depth":-1e2, "max depth":95e3, "temperature":500}],
     "composition models":[{"model":"uniform", "compositions":[1], "min depth":-1e2}]},

    {"model":"continental plate", "name":"Continent plate C", "coordinates":[[10,12],[20,12],[20,15],[10,15]], "min depth":-1e2,
     "topography models":[{"model":"uniform","topography":-10000}],
     "temperature models":[{"model":"uniform", "min depth":-1e2, "max depth":95e3, "temperature":500}],
     "composition models":[{"model":"uniform", "compositions":[2], "min depth":-1e2}]},

      // no topography and everywhere in the continent. Should not change any topography
    {"model":"continental plate", "name":"Continent plate D", "coordinates":[[-1,10],[41,10],[41,15],[-1,15]], "min depth":-1e2,
     "composition models":[{"model":"uniform", "compositions":[3], "min depth":-1e2, "max depth":25e3}]},

    {"model":"oceanic plate", "name":"oceanic plate A", "coordinates":[[-1,0],[41,0],[41,10],[-1,10]], "min depth":-1e2,
     "temperature models":[{"model":"half space model", "max depth":95e3, "min depth":-1e2, "spreading velocity":0.005, "ridge coordinates":[[[10,0], [10,5]],[[20,5],[20,10]]]}],
     "composition models":[{"model":"uniform", "compositions":[4], "min depth":-1e2}]},

    {"model":"oceanic plate", "name":"oceanic plate B", "coordinates":[[12.5,2.5],[17.5,2.5],[17.5,5],[12.5,5]], "min depth":-1e2,
     "topography models":[{"model":"uniform", "topography":10000}],
     "composition models":[{"model":"uniform", "compositions":[5], "min depth":-1e2}]}
  ]
}

[Djneu]

Sorry I didn't get a chance to look at the older pull request, I ended up getting sick last week! I think this sounds like a good way to do this, I think in the initial pull request we were just focused on getting a simple interface implemented so it was quite basic.

One of the main points of discussion if I recall was defining a surface type so that we could tell ASPECT whether it is adding topography by either thickening at the top, or in the case of isostasy, thickening at the bottom so that the entire lithosphere is uplifted. I didn't look through all the code, but was this in here somewhere?

[MFraters]

Sorry I didn't get a chance to look at the older pull request, I ended up getting sick last week!
Oh, I am sorry to hear that. I hope you feel better now.

One of the main points of discussion if I recall was defining a surface type so that we could tell ASPECT whether it is adding topography by either thickening at the top, or in the case of isostasy, thickening at the bottom so that the entire lithosphere is uplifted. I didn't look through all the code, but was this in here somewhere?

hmm, I may need to think a bit more about this, but what the world builder currently does is that is just works from depth 0, which is independent of the topography. With that, in gwb-grid I can compute the a depth from the surface and a depth from a reference, that happens here (notice the last two line):

WorldBuilder/source/gwb-grid/main.cc

Lines 945 to 967 in 4a5fa5d

	grid_x[counter] = x_min + (static_cast<double>(i) - 1.0) * dlong;
	const double longitude = grid_x[i];
	domain_height [counter]= outer_radius - inner_radius;
	cell_height[counter] = domain_height[counter] / static_cast<double>(n_cell_z);
	grid_z[counter] = inner_radius + (static_cast<double>(j) - 1.0) * cell_height[counter];
	grid_depth_wrt_surface[counter] = domain_height[counter] - (static_cast<double>(j) - 1.0) * cell_height[counter];
	const double radius = grid_z[counter];
	const double x = radius * std::cos(longitude);
	const double z = radius * std::sin(longitude);
	std::vector<std::array<unsigned ,3>> properties;
	properties.push_back({{6,0,0}}); // topography
	const std::array<double,2> coords = {{x,z}};
	const double topography = world->properties(coords, grid_depth_wrt_surface[counter],properties)[0];
	grid_x[counter] = x_min + (static_cast<double>(i) - 1.0) * dlong;
	domain_height [counter]= outer_radius + topography - inner_radius;
	cell_height[counter] = domain_height[counter] / static_cast<double>(n_cell_z);
	grid_z[counter] = inner_radius + (static_cast<double>(j) - 1.0) * cell_height[counter];
	grid_depth_wrt_surface[counter] = domain_height[counter] - (static_cast<double>(j) - 1.0) * cell_height[counter];
	grid_depth_wrt_reference[counter] = domain_height[counter] - (static_cast<double>(j) - 1.0) * cell_height[counter] - topography;
	counter++;

And later I pass in the depth with respect to the surface:

WorldBuilder/source/gwb-grid/main.cc

Line 1657 in 4a5fa5d

std::vector<double> output = world->properties(coords, grid_depth_wrt_surface[i],properties);

So how it is currently setup, the caller (gwb-grid or aspect) is responsible to put in the correct depth. This means that we could make a switch in aspect which depth to put into the world builder. But this would probably require to know the topography in aspect. I don't know whether this is desirable or feasible in aspect since it has been a while since I have looked at the topography and world builder code in aspect. Do you think this would work or not?

[Djneu]

I think this sounds good for now, and we can update it later if it is necessary when isostasy is implemented. From the third point on option 2 in #780, I think the idea was that we pass the topography separately and then ASPECT will decide how to apply the topography, but this was something that we still need to implement within ASPECT for it to work.