Add the finishing touches on importing global topography.

[bangerth]

This isn't quite read yet, but once done it should put the finishing touches on for #5421. I'm opening this now already because this way I can get to see whether the tests succeed.

[bangerth]

I still have to make sure the test output is what the CI check expects, but other than that this is complete.

I do have a data problem, though, that I do not know whether it's in the software or in the input file. If I plot the output of the new surface_elevation_spherical_shell_ascii_data_global_3d in VTU format and with some more refinement steps, I get these elevation maps:

This is clearly not correct: West Africa is not at an elevation of 3km and higher, and much of Europe is not (yet) under water. The high point of Earth is also not at 12170m, though the low point given as -11480m is pretty accurate. @alarshi any thoughts whether that must be a problem in my implementation or in the input file?

[bangerth]

Looking around some more, the abnormally high places seem to be around the equator, the abnormally low ones further north.

[bangerth]

That can't be it. The data file actually looks correct:

Help welcome!

[bangerth]

This is going to require a bit more work. The data file has values between approximately -6700 and +5600 meters of elevation. When I put this code into spherical_shell.cc,

      template <int dim>
      double
      SphericalManifoldWithTopography<dim>::
      topography_for_point(const Point<dim> &x_y_z) const
      {
        if (dynamic_cast<const InitialTopographyModel::ZeroTopography<dim>*>(topo) != nullptr)
          return 0;
        else
          {
            Assert (dim==3, ExcNotImplemented());

            // The natural coordinate system of the sphere geometry is r/phi/theta.
            // This is what we need to query the topography with. So start by
            // converting into this coordinate system
            const std::array<double, dim> r_phi_theta = Utilities::Coordinates::cartesian_to_spherical_coordinates(x_y_z);

            // Grab lon,lat coordinates
            Point<dim-1> surface_point;
            for (unsigned int d=0; d<dim-1; ++d)
              surface_point[d] = r_phi_theta[d+1];
            const double elevation = topo->value(surface_point);
            Assert (elevation < 10000, ExcInternalError());

            static double max=0, min=0;
            if (elevation > max)
              {
                std::cout << "New max elevation: " << elevation << std::endl;
                max = elevation;
              }
            if (elevation < min)
              {
                std::cout << "New min elevation: " << elevation << std::endl;
                min = elevation;
              }

            return elevation;
          }
      }

it reports min and max elevations as follows:

New min elevation: -6523.63
New max elevation: 5262.77

This is consistent. But when I put corresponding code into the surface_elevation postprocessor, I get these values:

New min elevation: -11050
New max elevation: 11136.6

That's clearly not right.

I would be ok if we merged this as-is so that at least @alarshi and @jdannberg can set up their models and try them out, even though we know that there's a bug somewhere. I will tackle the question next week then.

[gassmoeller]

I agree it is good to have the tests in and the functionality enabled, even though we know there is a bug somewhere in the already merged code. It is certainly not this PR that is causing the issue. Thanks for looking into it, maybe @alarshi can take a look as well. Just as a wild guess: Since the amplitudes are approximately twice of what we would expect, maybe the surface deformation is applied twice somewhere in the manifold?