deglaciation compared to seakon (draft)

publication_v1.0/test3/test3_compare3DGIA.jl

@@ -57,7 +57,7 @@ function mainplot(n, heterogeneous)
         title1 = L"Gaussian decrease of viscosity $\,$"
         title2 = L"Gaussian increase of viscosity $\,$"
     elseif heterogeneous == "none"
-        seakon_files = ["E0L4V4", "E0L4V4"]
+        seakon_files = ["E0L4V4", "E0L0V1"]
         fastiso_files = ["ref_$suffix.jld2", "no_litho_$suffix.jld2"]
         elims = (-20, 45)
         title1 = L"Homogeneous PREM configuration $\,$"

publication_v1.0/test5/compute5C.jl

@@ -0,0 +1,87 @@
+push!(LOAD_PATH, "../")
+using FastIsostasy
+using JLD2, NCDatasets, CairoMakie, Interpolations, DelimitedFiles
+
+not(x::Bool) = !x
+
+function lon360tolon180(lon, X)
+    permidx = lon .> 180
+    lon180 = vcat(lon[permidx] .- 360, lon[not.(permidx)])
+    Hice180 = cat(X[permidx, :, :], X[not.(permidx), :, :], dims=1)
+    return lon180, Hice180
+end
+
+function compute5C()
+    file = "../data/ICE6Gzip/IceT.I6F_C.131QB_VM5a_1deg.nc"
+    ds = NCDataset(file)
+    Hice = copy(ds["stgit"][:, :])
+    t = copy(ds["Time"][:, :])
+    lat = copy(ds["Lat"][:, :])
+    lon = copy(ds["Lon"][:, :])
+    close(ds)
+
+    lon180, Hice180 = lon360tolon180(lon, Hice)
+    fig, ax, hm = heatmap(Hice180[:, :, 2], colormap = :ice, colorrange = (1e-8, 4e3),
+        lowclip = :transparent)
+    Colorbar(fig[1, 2], hm)
+    hidedecorations!(ax)
+    fig
+
+    lon, Hice = lon180, Hice180
+    t .*= -1
+    nlon, nlat, nt = size(Hice)
+    Hsum = [sum(Hice[:, :, k]) for k in eachindex(t)]
+
+    Hlim = (1e-8, 4e3)
+    fig = Figure(resolution = (1600, 1000), fontsize = 24)
+    ax1 = Axis(fig[1:3, 1], aspect = DataAspect())
+    ax2 = Axis(fig[4, 1])
+    hidedecorations!(ax1)
+    hideydecorations!(ax2)
+    xlims!(ax2, extrema(t))
+    ylims!(ax2, extrema(Hsum) .+ (-2e6, 2e6))
+
+    kobs = Observable(1)
+    Hobs = @lift( Hice[:, :, $kobs] )
+    nobs = @lift( length(t[1:$kobs]) )
+    tobs = @lift( t[1:$nobs] )
+    Vproxy = @lift( Hsum[1:$nobs] )
+    hm = heatmap!(ax1, Hobs, colormap = :ice, colorrange = Hlim, lowclip = :transparent)
+    Colorbar(fig[1:3, 2], hm, height = Relative(0.6))
+    points = Observable(Point2f[(t[1], Hsum[1])])
+    lines!(ax2, points)
+
+    record(fig, "plots/test5/ICE6G-cylce.mp4", eachindex(t), framerate = 10) do k
+        kobs[] = k
+        new_point = Point2f(t[k], Hsum[k])
+        points[] = push!(points[], new_point)
+    end
+
+    Hitp = linear_interpolation((lon, lat, t), Hice, extrapolation_bc = Flat())
+
+    Omega = ComputationDomain(3500e3, 7)
+    c = PhysicalConstants()
+    lbmantle = c.r_equator .- [5.721, 5.179] .* 1e6
+    lb = vcat(96e3, lbmantle)
+    lv = [0.5, 1.58, 3.16] .* 1e21
+    p = LayeredEarth(Omega, layer_boundaries = lb[1:2], layer_viscosities = lv[1:2])
+
+    Hice_vec = [copy(Omega.null) for _ in t]
+    for k in eachindex(t)
+        for IJ in CartesianIndices(Omega.Lat)
+            i, j = Tuple(IJ)
+            # println(Omega.Lon[i, j], "  ", Omega.Lat[i, j], "  ", t[k])
+            Hice_vec[k][i, j] = Hitp(Omega.Lon[i, j], Omega.Lat[i, j], t[k])
+        end
+    end
+    deltaH = [Hice_vec[k] - Hice_vec[1] for k in eachindex(t)]
+
+    tsec = years2seconds.(t .* 1e3)
+    interactive_sl = true
+    fip = FastIsoProblem(Omega, c, p, tsec, interactive_sl, tsec, deltaH)
+    solve!(fip)
+    println("Computation took $(fip.out.computation_time) s")
+
+    @save "../data/test5/ICE6G/homogeneous-interactivesl=$interactive_sl-N="*
+        "$(Omega.Nx).jld2" t fip Hitp Hice_vec deltaH
+end

publication_v1.0/test5/laty.jl

@@ -0,0 +1,32 @@
+function load_laty_ICE6G()
+    latydir = "../data/Latychev/ICE6G/r1d_ICE6G"
+    latyfiles = readdir(latydir)
+    tlaty = [-120, -26, -21, -16, -11, -6, -2, -1, -0.125, 0, 0.125] .* 1e3
+
+    latyfiles = latyfiles[1:end-1]
+    tlaty = tlaty[1:end-1]
+    nlon, nlat = 512, 256
+    ulaty = zeros(nlon, nlat, length(tlaty))
+
+    for k in eachindex(latyfiles)
+        ulaty[:, :, k] = reshape(vec(readdlm(joinpath(latydir, latyfiles[k]))), nlon, nlat)
+    end
+
+    gllatlon_file = "../data/Latychev/ICE6G/gl256_LatLon"
+    LonLat, header = readdlm(gllatlon_file, header = true)
+    Lon = reshape(LonLat[:, 1], nlon, nlat)
+    Lat = reshape(LonLat[:, 2], nlon, nlat)
+    Lon, Lat = reverse(Lon, dims = 2), reverse(Lat, dims = 2)
+    ulaty = reverse(ulaty, dims = 2)
+
+    Lon, ulaty = lon360tolon180(Lon[:, 1], ulaty)
+
+    for k in eachindex(tlaty)
+        fig = heatmap(ulaty[:, :, k], colorrange = (-500, 500), colormap = :PuOr)
+        save("plots/test5/laty-ICE6G-global$(tlaty[k]).png", fig)
+    end
+    itp = linear_interpolation((Lon[:, 1], Lat[1, :], tlaty), ulaty,
+        extrapolation_bc = Flat())
+
+    return tlaty, ulaty, Lon, Lat, itp
+end

publication_v1.0/test5/plot5C.jl

@@ -0,0 +1,75 @@
+using FastIsostasy
+using JLD2, NCDatasets, CairoMakie, Interpolations, DelimitedFiles
+include("laty.jl")
+include("compute5C.jl")
+
+# @load "../data/test5/ICE6G/homogeneous-N=200.jld2" t fip Hitp Hice_vec deltaH
+@load "../data/test5/ICE6G/homogeneous-interactivesl=true-N=128.jld2" t fip Hitp Hice_vec deltaH
+
+Hlim = (1e-8, 4e3)
+kobs = Observable(1)
+fig = Figure(resolution = (1600, 1000), fontsize = 30)
+axs = [Axis(fig[1, j], title = @lift("t = $(t[$kobs]) ka"),
+    aspect = DataAspect()) for j in 1:2]
+[hidedecorations!(ax) for ax in axs]
+hm1 = heatmap!(axs[1], @lift(Hice_vec[$kobs]), colorrange = Hlim, lowclip = :transparent,
+    colormap = :ice)
+hm2 = heatmap!(axs[2], @lift(fip.out.u[$kobs]), colormap = :vik, colorrange = (-800, 800))
+Colorbar(fig[2, 1], hm1, vertical = false, flipaxis = false, width = Relative(0.8))
+Colorbar(fig[2, 2], hm2, vertical = false, flipaxis = false, width = Relative(0.8))
+record(fig, "plots/test5/ICE6G-cycle-displacement.mp4", eachindex(t), framerate = 10) do k
+    kobs[] = k
+end
+
+dHlims = (-2600, 2600)
+ulims = (-500, 500)
+kobs = Observable(1)
+fig = Figure(resolution = (1600, 1000), fontsize = 30)
+vars = [deltaH, fip.out.u .+ fip.out.ue]
+axs = [Axis3(fig[1, j], title = @lift("t = $(t[$kobs]) ka,"*
+ "  range = $(round.(extrema(vars[j][$kobs])))")) for j in 1:2]
+# [hidedecorations!(ax) for ax in axs]
+sf1 = surface!(axs[1], @lift(deltaH[$kobs]), colorrange = dHlims, colormap = :vik)
+sf2 = surface!(axs[2], @lift(fip.out.u[$kobs] + fip.out.ue[$kobs]), colormap = :PuOr,
+    colorrange = ulims)
+zlims!(axs[1], dHlims)
+zlims!(axs[2], ulims)
+Colorbar(fig[2, 1], sf1, vertical = false, flipaxis = false, width = Relative(0.8))
+Colorbar(fig[2, 2], sf2, vertical = false, flipaxis = false, width = Relative(0.8))
+record(fig, "plots/test5/ICE6G-cycle-surface.mp4", eachindex(t), framerate = 10) do k
+    kobs[] = k
+end
+
+tlaty, ulaty, Lon, Lat, itp = load_laty_ICE6G()
+
+mean_error = fill(Inf, length(tlaty))
+max_error = fill(Inf, length(tlaty))
+for k in eachindex(tlaty)
+    k_fastiso = argmin( (t .- tlaty[k]/1e3) .^ 2 )
+    ulatyitp = itp.(fip.Omega.Lon, fip.Omega.Lat, tlaty[k])
+    ufastiso = fip.out.u[k_fastiso] + fip.out.ue[k_fastiso]
+    mean_error[k] = mean( abs.(ulatyitp - ufastiso) )
+    max_error[k] = maximum( abs.(ulatyitp - ufastiso) )
+    tmpfig = Figure(resolution = (1800, 700), fontsize = 30)
+    axs = [Axis(tmpfig[1, j]) for j in 1:3]
+    [hidedecorations!(ax) for ax in axs]
+    hm1 = heatmap!(axs[1], ulatyitp, colorrange = (-500, 500), colormap = :PuOr)
+    hm2 = heatmap!(axs[2], ufastiso, colorrange = (-500, 500), colormap = :PuOr)
+    hm3 = heatmap!(axs[3], ulatyitp - ufastiso, colorrange = (-100, 100),
+        colormap = :lighttemperaturemap)
+    Colorbar(tmpfig[2, 1:2], hm1, label = "vertical displacement (m)", vertical = false,
+        width = Relative(0.4))
+    Colorbar(tmpfig[2, 3], hm3, label = L"$u_\mathrm{sk} - u_\mathrm{fi} $ (m)",
+        vertical = false, width = Relative(0.8))
+    save("plots/test5/laty-ICE6G-local$(tlaty[k]).png", tmpfig)
+end
+
+fig = Figure()
+ax = Axis(fig[1, 1], xlabel = "time (kyr)", ylabel = "absolute difference (m)")
+lines!(ax, tlaty, mean_error, label = "mean")
+lines!(ax, tlaty, max_error, label = "max")
+axislegend(ax)
+save("plots/test5/ICE6G-difftimeseries.png", fig)
+
+barplot(eachindex(tlaty), max_error)
+barplot(eachindex(tlaty), mean_error)