From 58fb988f923633536f42815e14a58697c7646e85 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 3 Dec 2024 16:49:57 -0700 Subject: [PATCH 1/3] Initial halfar-melt-ebm coupling and simulation --- examples/climate/ebm_melt.jl | 284 +++++++++++++++++++++++++++++++++++ 1 file changed, 284 insertions(+) create mode 100644 examples/climate/ebm_melt.jl diff --git a/examples/climate/ebm_melt.jl b/examples/climate/ebm_melt.jl new file mode 100644 index 00000000..9f43829c --- /dev/null +++ b/examples/climate/ebm_melt.jl @@ -0,0 +1,284 @@ +# Import Dependencies + +# AlgebraicJulia Dependencies +using Catlab +using CombinatorialSpaces +using DiagrammaticEquations +using DiagrammaticEquations.Deca +using Decapodes +using Decapodes: SchSummationDecapode + +# External Dependencies +using ComponentArrays +using CoordRefSystems +using GLMakie +using JLD2 +using LinearAlgebra +using MLStyle +using NetCDF +using OrdinaryDiffEq +using ProgressBars +Point3D = Point3{Float64} + +## Load a mesh +s_plots = loadmesh(Icosphere(7)); +s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s_plots); +subdivide_duals!(s, Barycenter()); +wireframe(s_plots) + +## Load Data +# Effective sea ice thickness data can be downloaded here: +# https://pscfiles.apl.uw.edu/axel/piomas20c/v1.0/monthly/piomas20c.heff.1901.2010.v1.0.nc +ice_thickness_file = "examples/climate/piomas20c.heff.1901.2010.v1.0.nc" + +# Use ncinfo(ice_thickness_file) to get information on variables. +# Sea ice thickness ("sit") has dimensions of [y, x, time]. +# y,x index into "Latitude" and "Longitude" variables. +# Time is in units of days since 1901-01-01. +lat = ncread(ice_thickness_file, "Latitude") +lon = ncread(ice_thickness_file, "Longitude") +sit = ncread(ice_thickness_file, "sit") + +# Convert latitude from [90, -90] to [0, 180] for convenience. +lat .= -lat .+ 90 + +p_sph = map(point(s)) do p + p = convert(Spherical, Cartesian(p...)) + [rad2deg(p.θ).val, rad2deg(p.ϕ).val] +end + +# TODO: Do algebraic parameterization, rather than nearest-neighbor interpolation. +# TODO: You can set a value to 0.0 if the distance to the nearest-neighbor is greater than some threshold. +sit_sph_idxs = map(ProgressBar(p_sph)) do p + argmin(map(i -> sqrt((lat[i] - p[1])^2 + (lon[i] - p[2])^2), eachindex(lat))) +end + +sit_sph = map(sit_sph_idxs, p_sph) do i, p + ((p[1] > maximum(lat)) || isnan(sit[i])) ? 0.0f0 : sit[i] +end + +f = Figure() +ax = LScene(f[1,1], scenekw=(lights=[],)) +msh = mesh!(ax, s_plots, color=sit_sph) +Colorbar(f[1,2], msh) + +## Define the model + +halfar_eq2 = @decapode begin + (h, melt)::Form0 + Γ::Form1 + n::Constant + + ∂ₜ(h) == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n-1)) * avg₀₁(h^(n+2))) - melt +end + +glens_law = @decapode begin + (A,Γ)::Form1 + (ρ,g,n)::Constant + + Γ == (2/(n+2))*A*(ρ*g)^n +end + +ice_dynamics_composition_diagram = @relation () begin + dynamics(Γ,n) + stress(Γ,n) +end + +ice_dynamics = apex(oapply(ice_dynamics_composition_diagram, + [Open(halfar_eq2, [:Γ,:n]), + Open(glens_law, [:Γ,:n])])) + +energy_balance = @decapode begin + (Tₛ, ASR, OLR, HT)::Form0 + C::Constant + + ∂ₜ(Tₛ) == (ASR - OLR + HT) ./ C +end + +absorbed_shortwave_radiation = @decapode begin + (Q, ASR)::Form0 + α::Constant + + ASR == (1 .- α) .* Q +end + +outgoing_longwave_radiation = @decapode begin + (Tₛ, OLR)::Form0 + (A,B)::Constant + + OLR == A .+ (B .* Tₛ) +end + +heat_transfer = @decapode begin + (HT, Tₛ)::Form0 + (D,cosϕᵖ,cosϕᵈ)::Constant + + HT == (D ./ cosϕᵖ) .* ⋆(d(cosϕᵈ .* ⋆(d(Tₛ)))) +end + +insolation = @decapode begin + Q::Form0 + cosϕᵖ::Constant + + Q == 450 * cosϕᵖ +end + +budyko_sellers_composition_diagram = @relation () begin + energy(Tₛ, ASR, OLR, HT) + absorbed_radiation(Q, ASR) + outgoing_radiation(Tₛ, OLR) + diffusion(Tₛ, HT, cosϕᵖ) + insolation(Q, cosϕᵖ) +end + +budyko_sellers = apex(oapply(budyko_sellers_composition_diagram, + [Open(energy_balance, [:Tₛ, :ASR, :OLR, :HT]), + Open(absorbed_shortwave_radiation, [:Q, :ASR]), + Open(outgoing_longwave_radiation, [:Tₛ, :OLR]), + Open(heat_transfer, [:Tₛ, :HT, :cosϕᵖ]), + Open(insolation, [:Q, :cosϕᵖ])])) + +warming = @decapode begin + Tₛ::Form0 + A::Form1 + + #A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e7) + #A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.01e-13) + A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.01e-19) +end + +melting = @decapode begin + (Tₛ, h, melt, water)::Form0 + Dₕ₂ₒ::Constant + + melt == (Tₛ - 15)*1e-16*h + ∂ₜ(water) == melt + Dₕ₂ₒ*Δ(water) +end + +budyko_sellers_halfar_water_composition_diagram = @relation () begin + budyko_sellers(Tₛ) + + warming(A, Tₛ) + + melting(Tₛ, h, melt) + + halfar(A, h, melt) +end + +budyko_sellers_halfar_water = apex(oapply(budyko_sellers_halfar_water_composition_diagram, + [Open(budyko_sellers, [:Tₛ]), + Open(warming, [:A, :Tₛ]), + Open(melting, [:Tₛ, :h, :melt]), + Open(ice_dynamics, [:stress_A, :dynamics_h, :dynamics_melt])])) + +## Define constants, parameters, and initial conditions + +# This is a primal 0-form, with values at vertices. +cosϕᵖ = map(x -> cos(x[1]), point(s)) +# This is a dual 0-form, with values at edge centers. +cosϕᵈ = map(edges(s)) do e + (cos(point(s, src(s, e))[1]) + cos(point(s, tgt(s, e))[1])) / 2 +end + +α₀ = 0.354 +α₂ = 0.25 +α = map(point(s)) do ϕ + α₀ + α₂*((1/2)*(3*ϕ[1]^2 - 1)) +end +A = 210 +B = 2 +f = 0.70 +ρ = 1025 +cw = 4186 +H = 70 +C = map(point(s)) do ϕ + f * ρ * cw * H +end +D = 0.6 + +# Isothermal initial conditions: +Tₛ₀ = map(point(s)) do ϕ + 15.0 +end + +water = map(point(s)) do _ + 0.0 +end + +Dₕ₂ₒ = 1e-16 + +n = 3 +halfar_ρ = 910 +g = 9.8 + +h₀ = sit_sph +# Store these values to be passed to the solver. +u₀ = ComponentArray( + Tₛ = Tₛ₀, + h = h₀, + melting_water = water) + +constants_and_parameters = ( + budyko_sellers_absorbed_radiation_α = α, + budyko_sellers_outgoing_radiation_A = A, + budyko_sellers_outgoing_radiation_B = B, + budyko_sellers_energy_C = C, + budyko_sellers_diffusion_D = D, + budyko_sellers_cosϕᵖ = cosϕᵖ, + budyko_sellers_diffusion_cosϕᵈ = cosϕᵈ, + halfar_n = n, + halfar_stress_ρ = halfar_ρ, + halfar_stress_g = g, + melting_Dₕ₂ₒ = Dₕ₂ₒ) + +# Define how symbols map to Julia functions + +function generate(sd, my_symbol; hodge=GeometricHodge()) + op = @match my_symbol begin + :♯ => begin + sharp_mat = ♯_mat(sd, AltPPSharp()) + x -> sharp_mat * x + end + :mag => x -> begin + norm.(x) + end + :^ => (x,y) -> x .^ y + :* => (x,y) -> x .* y + x => error("Unmatched operator $my_symbol") + end + return (args...) -> op(args...) +end + +## Generate simulation + +sim = eval(gensim(budyko_sellers_halfar_water)) +fₘ = sim(s, generate) + +## Run simulation + +tₑ = 100.0 + +# Julia will precompile the generated simulation the first time it is run. +@info("Precompiling Solver") +prob = ODEProblem(fₘ, u₀, (0, 1e-4), constants_and_parameters) +soln = solve(prob, Tsit5()) +soln.retcode != :Unstable || error("Solver was not stable") + +@info("Solving") +prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) +soln = solve(prob, Tsit5()) +@show soln.retcode +@info("Done") + +extrema(soln(0.0).halfar_h) +extrema(soln(tₑ).halfar_h) + +@save "budyko_sellers_halfar_water.jld2" soln + +# Visualize + +g = Figure() +ax = LScene(g[1,1], scenekw=(lights=[],)) +msh = mesh!(ax, s_plots, color=soln.u[end].h) +Colorbar(g[1,2], msh) + From c152008402123ac99fee683f6117cce2468a3342 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 3 Dec 2024 21:42:47 -0700 Subject: [PATCH 2/3] Turn script into Documenter page --- docs/Project.toml | 1 + docs/make.jl | 3 +- .../src/ebm_melt/ebm_melt.md | 106 +++++++++++++----- 3 files changed, 81 insertions(+), 29 deletions(-) rename examples/climate/ebm_melt.jl => docs/src/ebm_melt/ebm_melt.md (65%) diff --git a/docs/Project.toml b/docs/Project.toml index b29fb873..1dda9289 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -18,6 +18,7 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" +NetCDF = "30363a11-5582-574a-97bb-aa9a979735b9" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" diff --git a/docs/make.jl b/docs/make.jl index e8d7e32b..6762c18c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -44,7 +44,7 @@ makedocs( pagesonly = true, linkcheck = true, linkcheck_ignore = [r"agupubs\.onlinelibrary\.wiley\.com", # This gives a 403 Forbidden - r"Decapodes\.jl/dev"], # 404, probably due to bad self-rerference + r"Decapodes\.jl/dev"], # 404, probably due to bad self-reference pages = Any[ "Decapodes.jl" => "index.md", "Overview" => "overview/overview.md", @@ -57,6 +57,7 @@ makedocs( "Glacial Flow" => "ice_dynamics/ice_dynamics.md", "Grigoriev Ice Cap" => "grigoriev/grigoriev.md", # Requires ice_dynamics "Budyko-Sellers-Halfar" => "bsh/budyko_sellers_halfar.md", # Requires ice_dynamics + "Halfar-EBM-Water" => "ebm_melt/ebm_melt.md", "Halfar-NS" => "halmo/halmo.md", # Requires grigoriev "NHS" => "nhs/nhs_lite.md", "Pipe Flow" => "poiseuille/poiseuille.md", diff --git a/examples/climate/ebm_melt.jl b/docs/src/ebm_melt/ebm_melt.md similarity index 65% rename from examples/climate/ebm_melt.jl rename to docs/src/ebm_melt/ebm_melt.md index 9f43829c..695458ff 100644 --- a/examples/climate/ebm_melt.jl +++ b/docs/src/ebm_melt/ebm_melt.md @@ -1,37 +1,49 @@ +# Haflar-EBM-Water + +This docs page demonsrates a composition of the Halfar model of ice dynamics with the Budyko-Sellers energy-balance model (EBM) defining temperature dynamics. Surface temperature affects the rate at which ice diffuses, and melts ice into a water density term. This water density then diffuses across the domain. + +We execute these dynamics on real sea ice thickness data provided by the Polar Science Center. + +``` @example DEC # Import Dependencies # AlgebraicJulia Dependencies using Catlab +using Catlab.Graphics using CombinatorialSpaces using DiagrammaticEquations -using DiagrammaticEquations.Deca using Decapodes -using Decapodes: SchSummationDecapode # External Dependencies using ComponentArrays using CoordRefSystems -using GLMakie -using JLD2 +using CairoMakie using LinearAlgebra using MLStyle using NetCDF using OrdinaryDiffEq -using ProgressBars Point3D = Point3{Float64} +``` -## Load a mesh +First, we load a mesh. We will execute these dynamics on the sphere: + +``` @example DEC +# Load a mesh s_plots = loadmesh(Icosphere(7)); s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s_plots); subdivide_duals!(s, Barycenter()); wireframe(s_plots) +``` + +## Load data + +The data provided by the Polar Science Center is given as a NetCDF file. Ice thickness is a matrix with the same dimensions as a matrix provided Latitude and Longitude of the associated point on the Earth's surface. We need to convert between polar and Cartesian coordinates to use this data on our mesh. -## Load Data -# Effective sea ice thickness data can be downloaded here: -# https://pscfiles.apl.uw.edu/axel/piomas20c/v1.0/monthly/piomas20c.heff.1901.2010.v1.0.nc -ice_thickness_file = "examples/climate/piomas20c.heff.1901.2010.v1.0.nc" +``` @example DEC +ice_thickness_file = "piomas20c.heff.1901.2010.v1.0.nc" +run(`curl -o $ice_thickness_file https://pscfiles.apl.uw.edu/axel/piomas20c/v1.0/monthly/piomas20c.heff.1901.2010.v1.0.nc`) -# Use ncinfo(ice_thickness_file) to get information on variables. +# Use ncinfo(ice_thickness_file) to interactively get information on variables. # Sea ice thickness ("sit") has dimensions of [y, x, time]. # y,x index into "Latitude" and "Longitude" variables. # Time is in units of days since 1901-01-01. @@ -42,14 +54,15 @@ sit = ncread(ice_thickness_file, "sit") # Convert latitude from [90, -90] to [0, 180] for convenience. lat .= -lat .+ 90 +# Convert mesh points from Cartesian to spherical coordinates. p_sph = map(point(s)) do p p = convert(Spherical, Cartesian(p...)) [rad2deg(p.θ).val, rad2deg(p.ϕ).val] end -# TODO: Do algebraic parameterization, rather than nearest-neighbor interpolation. -# TODO: You can set a value to 0.0 if the distance to the nearest-neighbor is greater than some threshold. -sit_sph_idxs = map(ProgressBar(p_sph)) do p +# Note: You can instead use an algebraic parameterization, rather than nearest-neighbor interpolation. +# Note: You can alternatively set a value to 0.0 if the distance to the nearest-neighbor is greater than some threshold. +sit_sph_idxs = map(p_sph) do p argmin(map(i -> sqrt((lat[i] - p[1])^2 + (lon[i] - p[2])^2), eachindex(lat))) end @@ -61,9 +74,14 @@ f = Figure() ax = LScene(f[1,1], scenekw=(lights=[],)) msh = mesh!(ax, s_plots, color=sit_sph) Colorbar(f[1,2], msh) +f +``` ## Define the model +Here, the Halfar model of ice dynamics is recalled, as well as the Budyko-Sellers EBM. These these models are composed individually. They are then coupled together via `warming` and `melting` components. + +``` @example DEC halfar_eq2 = @decapode begin (h, melt)::Form0 Γ::Form1 @@ -137,13 +155,14 @@ budyko_sellers = apex(oapply(budyko_sellers_composition_diagram, Open(outgoing_longwave_radiation, [:Tₛ, :OLR]), Open(heat_transfer, [:Tₛ, :HT, :cosϕᵖ]), Open(insolation, [:Q, :cosϕᵖ])])) +nothing # hide +``` +``` @example DEC warming = @decapode begin Tₛ::Form0 A::Form1 - #A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e7) - #A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.01e-13) A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.01e-19) end @@ -164,15 +183,21 @@ budyko_sellers_halfar_water_composition_diagram = @relation () begin halfar(A, h, melt) end +``` +``` @example DEC budyko_sellers_halfar_water = apex(oapply(budyko_sellers_halfar_water_composition_diagram, [Open(budyko_sellers, [:Tₛ]), Open(warming, [:A, :Tₛ]), Open(melting, [:Tₛ, :h, :melt]), Open(ice_dynamics, [:stress_A, :dynamics_h, :dynamics_melt])])) +``` + +## Define initial conditions -## Define constants, parameters, and initial conditions +The initial data must be specified for state variables, as well as constants and parameters. +``` @example DEC # This is a primal 0-form, with values at vertices. cosϕᵖ = map(x -> cos(x[1]), point(s)) # This is a dual 0-form, with values at edge centers. @@ -248,17 +273,20 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) end return (args...) -> op(args...) end +``` + +## Generate and run simulation -## Generate simulation +The composed model is generated and executed for 100 years. The model is run twice to demonstrate the speed of the model after the simulation code is precompiled by the first run. +We will save the final ice thickness data in a .jld2 file, an HDF5-compatible file format. We will also save the latitude and longitude of the points on the sphere. + +``` @example DEC sim = eval(gensim(budyko_sellers_halfar_water)) fₘ = sim(s, generate) -## Run simulation - tₑ = 100.0 -# Julia will precompile the generated simulation the first time it is run. @info("Precompiling Solver") prob = ODEProblem(fₘ, u₀, (0, 1e-4), constants_and_parameters) soln = solve(prob, Tsit5()) @@ -270,15 +298,37 @@ soln = solve(prob, Tsit5()) @show soln.retcode @info("Done") -extrema(soln(0.0).halfar_h) -extrema(soln(tₑ).halfar_h) +save("ice.jld2", + Dict("lat" => map(x -> x[1], p_sph), "lon" => map(x -> x[2], p_sph), "ice" => soln(tₑ).h)) + +(extrema(soln(0.0).h), extrema(soln(tₑ).h)) +``` -@save "budyko_sellers_halfar_water.jld2" soln +## Visualize -# Visualize +Let's visualize the initial conditions for ice height and the ice height after 100 years. -g = Figure() -ax = LScene(g[1,1], scenekw=(lights=[],)) +### Initial ice height + +``` @example DEC +f = Figure() +ax = LScene(f[1,1], scenekw=(lights=[],)) msh = mesh!(ax, s_plots, color=soln.u[end].h) -Colorbar(g[1,2], msh) +Colorbar(f[1,2], msh) +f +``` + +### Ice height after 100 years + +``` @example DEC +f = Figure() +ax = LScene(f[1,1], scenekw=(lights=[],)) +msh = mesh!(ax, s_plots, color=soln.u[end].h) +Colorbar(f[1,2], msh) +f +``` + +```@example DEC +run(`rm $ice_thickness_file) # hide +``` From abc59404f866728f76c71b54a638ba2af3e291d3 Mon Sep 17 00:00:00 2001 From: Luke Morris <70283489+lukem12345@users.noreply.github.com> Date: Tue, 3 Dec 2024 21:53:21 -0700 Subject: [PATCH 3/3] Add missing backtick --- docs/src/ebm_melt/ebm_melt.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/ebm_melt/ebm_melt.md b/docs/src/ebm_melt/ebm_melt.md index 695458ff..520a582f 100644 --- a/docs/src/ebm_melt/ebm_melt.md +++ b/docs/src/ebm_melt/ebm_melt.md @@ -329,6 +329,6 @@ f ``` ```@example DEC -run(`rm $ice_thickness_file) # hide +run(`rm $ice_thickness_file`) # hide ```