Solve typo in Project

Project.toml

@@ -18,7 +18,7 @@ Infiltrator = "5903a43b-9cc3-4c30-8d17-598619ec4e9b"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-Lux = "b210s8857-7c20-44ae-9111-449ecde12c47"
+Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"

examples/SIA2D.jl

@@ -0,0 +1,91 @@
+using OrdinaryDiffEq
+# using DifferentialEquations
+
+####  Non-allocating functions  ####
+
+nx, ny = 200, 200   # elements in the spatial grid
+Δx, Δy = 100, 100     # size in meters of each box in the grid
+n = 3.0
+ρ = 900.0
+g = 9.81
+A = 8e-18
+Γ = 2.0 * A * (ρ * g)^n / (n+2) # 1 / m^3 s 
+
+diff_x(I, Δx) = (I[begin+1:end,:] - I[1:end-1,:]) / Δx
+diff_y(I, Δy) = (I[:,begin+1:end] - I[:,1:end - 1]) / Δy
+avg(I) = (I[1:end-1,1:end-1] + I[2:end,1:end-1] + I[1:end-1,2:end] + I[2:end,2:end]) * 0.25
+avg_x(I) = (I[1:end-1,:] + I[2:end,:]) * 0.5
+avg_y(I) = (I[:,1:end-1] + I[:,2:end]) * 0.5
+inn(I) = I[2:end-1,2:end-1]
+
+year2seg = 365.25 * 24 * 60 * 60
+# tspan = [0.0, 10 * year2seg]
+tspan = [0.0, 30.0]
+
+B = zeros((nx, ny))  # elevation of the bed of the glacier
+H₀ = [ .15 * ( 70^2 - ( (i - nx/2)^2 + (j - ny/2)^2 )  ) for i in 1:nx, j in 1:ny ] # initial ice thickness
+H₀ .= max.(0.0, H₀)
+dH = zeros(size(H₀)) 
+p = nothing
+
+function SIA!(dH, H, p, t)
+    S = B .+ H
+
+    # All grid variables computed in a staggered grid
+    # Compute surface gradients on edges
+    dSdx = diff_x(S, Δx)  
+    dSdy = diff_y(S, Δy) 
+    ∇Sx = avg_y(dSdx)
+    ∇Sy = avg_x(dSdy)
+    ∇S  = (∇Sx.^2 .+ ∇Sy.^2).^((n - 1)/2) 
+
+    H̄ = avg(H)
+
+    D = Γ .* H̄.^(n + 2) .* ∇S
+
+    # Compute flux components
+    dSdx_edges = diff_x(S[:,2:end - 1], Δx)
+    dSdy_edges = diff_y(S[2:end - 1,:], Δy)
+
+    # Cap surface elevaton differences with the upstream ice thickness to 
+    # imporse boundary condition of the SIA equation
+    # η₀ = params.physical.η₀
+    # dSdx_edges .= @views @. min(dSdx_edges,  η₀ * H[2:end, 2:end-1]/Δx)
+    # dSdx_edges .= @views @. max(dSdx_edges, -η₀ * H[1:end-1, 2:end-1]/Δx)
+    # dSdy_edges .= @views @. min(dSdy_edges,  η₀ * H[2:end-1, 2:end]/Δy)
+    # dSdy_edges .= @views @. max(dSdy_edges, -η₀ * H[2:end-1, 1:end-1]/Δy)
+
+    Dx = avg_y(D)
+    Dy = avg_x(D)
+    Fx = -Dx .* dSdx_edges
+    Fy = -Dy .* dSdy_edges 
+
+    #  Flux divergence
+    Fxx = diff_x(Fx, Δx)
+    Fyy = diff_y(Fy, Δy)
+    dH .= zeros(size(H))
+    dH[2:end-1, 2:end-1] .= -(Fxx .+ Fyy) 
+    @show t
+end
+
+prob = ODEProblem(SIA!, H₀, tspan, p)
+sol = solve(prob, RDPK3Sp35(), reltol = 1e-8, abstol = 1e-8, saveat=[0.0, 15.0, 30.0])
+
+using Plots
+
+H₁ = sol.u[end]
+
+
+heatmap_diff = Plots.heatmap(H₀ .- H₁, title="Variation in Ice Thickness", aspect_ratio=1, colorrange=(-200,200), colormap=:balance, clim=(-200, 200), axis=([], false))
+Plots.savefig(heatmap_diff, "variation.pdf")
+
+surface = Plots.surface(Δx * collect(1:nx), Δy * collect(1:ny), H₀, alpha=0.5, color=:green)
+
+x = Δx * collect(1:nx) / 1000 # [km]
+y1 = H₀[:, 100]
+y2 = sol.u[2][:,100]
+y3 = H₁[:, 100]
+line_plot = plot(x, [y1 y2 y3], label=["Initial condition (t=0)" "Intermediate solution (t=15)" "Final solution (t=30)"], title="Example of solutions of SIA equation", linewidth=3)
+xlabel!("Distance (km)")
+ylabel!("Heigth (m)")
+Plots.pdf(line_plot, "line_plot.pdf")

examples/Torsvik_2012/plots/plot_results.jl

@@ -0,0 +1,140 @@
+using Pkg; Pkg.activate(".")
+
+using SphereUDE
+using Serialization
+using Plots 
+using Plots.PlotMeasures
+using JLD2
+using LinearAlgebra
+
+JLD2.@load "examples/Torsvik_2012/results/results_dict.jld2" 
+
+data_directions_sph = cart2sph(results_dict["directions"], radians=false)
+fit_directions_sph = cart2sph(results_dict["fit_directions"], radians=false)
+
+data_lat = data_directions_sph[1,:]
+fit_lat  = fit_directions_sph[1,:]
+
+data_lon = data_directions_sph[2,:]
+fit_lon  = fit_directions_sph[2,:]
+
+α_error = 140 ./ sqrt.(results_dict["kappas"])
+
+# Latitude Plot 
+
+Plots.scatter(results_dict["times"], data_lat, label="Paleopole Latitudes", yerr=α_error, c=:lightsteelblue2, markersize=5)
+plot_latitude = Plots.plot!(results_dict["fit_times"], fit_lat, label="Estimated APWP using SphereUDE", 
+                    xlabel="Age (Myr)", yticks=[-90, -60, -30, 0, 30, 60],
+                    ylabel="Latitude (degrees)", ylims=(-90,60), lw = 4, c=:brown,
+                    legend=:topleft)
+plot!(fontfamily="Computer Modern",
+     #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    legend=true,
+    margin= 7mm,
+    size=(1200,500),
+    dpi=600)
+
+Plots.savefig(plot_latitude, "examples/Torsvik_2012/plots/latitude.pdf")
+
+### Longitude Plot 
+
+α_error_lon = α_error ./ cos.(π ./ 180. .* data_lat)
+
+Plots.scatter(results_dict["times"], data_lon, label="Paleopole Longitudes", yerr=α_error_lon, c=:lightsteelblue2, markersize=5)
+plot_longitude = Plots.plot!(results_dict["fit_times"], fit_lon, label="Estimated APWP using SphereUDE", 
+                    xlabel="Age (Myr)", yticks=[-180, -120, -60 , 0, 60, 120, 180],
+                    ylabel="Longitude (degrees)", ylims=(-180,180), lw = 4, c=:brown,
+                    legend=:bottomright)
+plot!(fontfamily="Computer Modern",
+     #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 7mm,
+    size=(1200,500),
+    dpi=600)
+
+Plots.savefig(plot_longitude, "examples/Torsvik_2012/plots/longitude.pdf")
+
+### Lat and long combined
+
+combo_plot = plot(plot_latitude, plot_longitude, layout = (2, 1))
+plot!(fontfamily="Computer Modern",
+     #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 7mm,
+    size=(1200,800),
+    dpi=600)
+Plots.savefig(combo_plot, "examples/Torsvik_2012/plots/latitude_longitude.pdf")
+
+
+### Angular velocity Plot
+
+angular_velocity = mapslices(x -> norm(x), results_dict["fit_rotations"], dims=1)[:]
+angular_velocity_path = [norm(cross(results_dict["fit_directions"][:,i], results_dict["fit_rotations"][:,i] )) for i in axes(results_dict["fit_rotations"],2)]
+
+plot_angular = Plots.plot(results_dict["fit_times"], angular_velocity, label="Maximum total angular velocity", 
+                    xlabel="Age (Myr)", 
+                    ylabel="Angular velocity (degrees/My)", lw = 5, c=:darkgreen,
+                    legend=:topleft)
+plot!(results_dict["fit_times"], angular_velocity_path, label="Pole angular velocity", lw = 4, c=:lightseagreen, ls=:dot)
+plot!(fontfamily="Computer Modern",
+    #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 10mm,
+    size=(1200,500),
+    dpi=300)
+
+
+Plots.savefig(plot_angular, "examples/Torsvik_2012/plots/angular.pdf")
+
+
+### Loss function
+
+losses = results_dict["losses"]
+
+plot_loss = Plots.plot(1:length(losses), losses, label="Loss Function", 
+                    xlabel="Epoch", 
+                    ylabel="Loss", lw = 5, c=:indigo,
+                    yaxis=:log,
+                    # yticks=[1,10,100],
+                    xlimits=(0,10000),
+                    ylimits=(10, 1000),
+                    xticks=[0,2500,5000,7500,10000],
+                    legend=:topright)
+
+vspan!(plot_loss, [0,5000], color = :navajowhite3, alpha = 0.2, labels = "ADAM");
+vspan!(plot_loss, [5000,10000], color = :navajowhite4, alpha = 0.2, labels = "BFGS");
+
+plot!(fontfamily="Computer Modern",
+    #title="PIL51",
+    titlefontsize=18,
+    tickfontsize=15,
+    legendfontsize=15,
+    guidefontsize=18,
+    #ylimits=(0.1,10),
+    #xlimits=(10^(-4),10^(-1)),
+    margin= 10mm,
+    size=(700,500),
+    dpi=300)
+
+Plots.savefig(plot_loss, "examples/Torsvik_2012/plots/loss.pdf")

examples/double_rotation/2rotations.ipynb

@@ -29,7 +29,7 @@
     "using SciMLSensitivity\n",
     "using Optimization, OptimizationOptimisers, OptimizationOptimJL\n",
     "\n",
-    "using SphereFit"
+    "using SphereUDE"
    ]
   },
   {