Analytical benchmark

Explicit time stepping (recommended)

In this example, we want to compute the viscous displacement of the upper mantle resulting from a cylindrical ice load with a radius of 1000 km and a thickness of 1 km. To do so, we first generate the computation domain and the load:

using FastIsostasy, CairoMakie

W, n = 3f6, 7   # square domain with size 2*W and 2^n points in each dimension
domain = RegionalDomain(W, n)

H_ice_0 = zeros(domain)             # Load: 0 at the beginning...
H_ice_1 = 1f3 .* (domain.R .< 1f6)  # cylinder afterwards
fig = plot_load(domain, H_ice_1)

This looks good! plot_load is a function of FastIsostasy that, similar to all functions plot_*, is only loaded if the user is using Makie. These functions are here to help the user visualize quickly the results that were obtained and are summarized in the API reference.

First, we wrap the ice load as a time interpolator so that the boundary doncitions of FastIsostasy can be updated internally (if you want to update the ice thickness externally, have a look at this). Then we define a solid earth structure, a sea-level model (that is inactive by default) and an output object. The simulation is subsequently executed with run!(). The Earth structure deserves particular attention, since it determines key properties of the problem. We here define a single layer boundary at 88 km depth. This means that the lithosphere is 88 km thick and the mantle below it is assumed to be a Maxwell body with a viscosity of 1f21 Pa s. The lithosphere density is set to 0, such that its displacement does not contribute to the pressure term of the isostatic adjustment (this is done to match Bueler et al., 2007).

t_ice = [0, 1, 50f3]
H_ice = [H_ice_0, H_ice_1, H_ice_1]
it = TimeInterpolatedIceThickness(t_ice, H_ice, domain)

bcs = BoundaryConditions(domain, ice_thickness = it)

solidearth = SolidEarth(    # Use the same geometry as in Bueler et al. (2007).
    domain,
    layer_boundaries = [88f3],
    layer_viscosities = [1f21],
    rho_litho = 0f0,
)

sealevel = SeaLevel()               # Will be inactive if not specified.

nout = NativeOutput(vars = [:u],    # only store viscous displacement.
    t = [100, 500, 1500, 5000, 10_000, 50_000f0])

sim = Simulation(domain, bcs, sealevel, solidearth; nout = nout)
run!(sim)
fig = plot_()

The viscous displacement field that we obtained is quite intuitive: its shape is largely determined by the forcing and its maximal amplitude is a fraction of the ice thickness is about \dfrac{ \rho\mathrm{ice} }{ \rho\mathrm{mantle} }. This simple example presents the advantage to have an analytical solution, which can be plotted along our solution:

fig = plot_transect(sim, [:u], analytic_cylinder_solution = true)

Our numerical solution yields very small error!

FastIsostasy was designed to be computationally efficient. We are therefore particularly interested in the time needed for the computation, which is stored in the nout object of the simulation.

println("Computation time using explicit time stepping: $(sim.nout.computation_time)")

Computation time using explicit time stepping: 0.9502461

This is the (compilation + computation) time that was required to compute 50 kyr of viscous displacement with a domain of 128x128 points!

Implicit time stepping

If the Earth structure is laterally constant (i.e. the lithospheric thickness and the mantle viscosity do not vary in x and y), the performance can be improved by using an implicit time stepping, as derived by Bueler et al. (2007). This can be achieved by specifying the lithosphere as RigidLithosphere() or as LaterallyConstantLithosphere() and requires to set a fixed time step via the DiffEqOptions in the SolverOptions:

solidearth = SolidEarth(
    domain,
    lithosphere = RigidLithosphere(),
    layer_boundaries = [88f3],
    layer_viscosities = [1f21],
)
opts = SolverOptions(diffeq = DiffEqOptions(alg = Euler(), dt_min = 100f0))
sim_implicit = Simulation(domain, bcs, sealevel, solidearth; nout = nout, opts = opts)
run!(sim_implicit)
fig = plot_transect(sim_implicit, [:u], analytic_cylinder_solution = true)

The results appear to be comparable with the previously obtained ones! However, the computation time is lower:

println("Computation time using implicit time stepping: $(sim_implicit.nout.computation_time)")

Computation time using implicit time stepping: 2.745924

All of the computations shown above are performed with Float32 as floating point precision. It is easy to switch to Float64, by passing arguments as such (e.g. W = 3e6). This will increase the accuracy of the results, but also the computation time. The default Float32 is sufficient for most applications.

Of course, this example remains simple. If you want to learn how to increase the complexity of your simulations, go to the next examples!