Start up your Matlab(R) from within the icetools-x.x directory and run
>> icetools_benchmark
This will solve a slab of ice 1000m long and 100m thick with a tilt of 5 deg using the analytical solution to this problem as inflow and outflow boundary conditions. The surface is described as a free surface and a no-slip condition is implemented at the bed boundary. The coordinate system is defined with positive x-axis downstream the ice flow an positive z-axis normal to the bedrock upward. The current icetools default setup will iterate 25 times to get a good viscosity estimate in this test case.
The mesh used in the icetools_benchmark example with 278 triangular elements.
The results of the each iteration are stored in the local directory Icetools is installed in and after the last iteration the result is displayed within Matlab(R). Each iteration is stored in ".mat" files called run_001_X.mat, where X is the number of iterations taken to create the data. There is also an initial run which is called run_001_i.mat.
Icetools provides a function called plot_result.m which can be used to compare calculated velocity fields with the analytical solution. Run:
>> plot_result('001_1')
to get the flow field and the comparison with the analytical solution for the first timestep. One can easily see that there is quite a difference still between the numerics and the analytical solution due to the fact that only one viscosity iteration was taken.
The flowfield and error estimate after one iteration.
By running
>> plot_result('001_25')
you can display the result after 25 iterations, which is now quite better then the initial result. The figure below demonstrates nicely how gridsize influences the result of the numerical model. Play with the gridsize (see below) to see even more improvements to the numerical solution to the given problem.
The flowfield and error estimate after 25 iterations.
The model parameters, like initial viscosity of the glacier, A and n from Glen's flow law, number of iterations and other parameters can be modified within the icetools_benchmark.m file. I assume with a decent knowledge of Matlab(R) it should be easy to understand the source code. Also the documentation of the Getfem++ package might help.
The mesh is based on the test.geo file, which is the input file for Gmsh. Gmsh than creates a file called test.msh, which is the Gmsh version of the computational mesh. Getfem++ had problems with importing Gmsh mesh files in 2D cases but never problems with mesh files created with GiD, another mesh generator. To avoid problems in 2D cases a Python script with the name gmsh2gid comes along with Icetools to convert a Gmsh mesh into a GiD mesh. The GiD mesh can than be used in Icetools as an input. To mesh the test.geo file and convert the mesh to GiD format run:
$ gmsh -2 -o test.msh -format msh1 test.geo $ ./gmsh2gid test.msh test.gid
The first command uses gmsh to create the mesh and saves it in gmsh mesh format version 1 into the test.msh file. Thereafter gmsh2gid is used to convert the mesh into the GiD format. At the moment gmsh2gid only supports gmsh mesh format version 1, so make sure you always save meshes in that format. If you want to play with gridsizes, there is a command line flag called -clscale, which changes the characteristic length scale of the mesh. A command like the one below would change the gridsize to 1/2 of the original grid size of 20 m.
$ gmsh -2 -o test.msh -format msh1 -clscale 0.5 test.geo
The general syntax of gmsh2gid is:
$ ./gmsh2gid inputfile outputfile
Now you can go and play with parameters in the Icetools code and change grid sizes of the mesh. Gmsh has an extensive documentation, so it should be easy to learn how to make complex geometries for Icetools.