diff --git a/.buildkite/pipeline.yml b/.buildkite/pipeline.yml index 13271e4e..30e8cc1b 100644 --- a/.buildkite/pipeline.yml +++ b/.buildkite/pipeline.yml @@ -12,7 +12,7 @@ steps: - label: ":scroll: Build docs and run tests" command: - - "srun --cpus-per-task=16 --mem=8G --time=1:00:00 --output=.buildkite/build_%j.log --unbuffered .buildkite/jobscript.sh" + - "srun --cpus-per-task=16 --mem=64G --time=3:00:00 --output=.buildkite/build_%j.log --unbuffered .buildkite/jobscript.sh" env: JULIA_PROJECT: "docs/" diff --git a/docs/src/Icethickness_Grigoriev_ice_cap_2021.tif b/docs/assets/Icethickness_Grigoriev_ice_cap_2021.tif similarity index 100% rename from docs/src/Icethickness_Grigoriev_ice_cap_2021.tif rename to docs/assets/Icethickness_Grigoriev_ice_cap_2021.tif diff --git a/docs/make.jl b/docs/make.jl index 6403164b..7b7dee57 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -16,23 +16,49 @@ if !(haskey(ENV, "GITHUB_ACTIONS") || haskey(ENV, "GITLAB_CI")) config["repo_root_url"] = "https://github.com/AlgebraicJulia/Decapodes.jl/blob/main/docs" end -# const literate_dir = joinpath(@__DIR__, "..", "examples") -# const generated_dir = joinpath(@__DIR__, "src", "examples") +const literate_dir = joinpath(@__DIR__, "..", "examples") +const generated_dir = joinpath(@__DIR__, "src", "examples") -# @info "Building literate files" -# for (root, dirs, files) in walkdir(literate_dir) -# out_dir = joinpath(generated_dir, relpath(root, literate_dir)) -# # @showprogress pmap(files) do file -# for file in files -# f,l = splitext(file) -# if l == ".jl" && !startswith(f, "_") -# Literate.markdown(joinpath(root, file), out_dir; -# config=config, documenter=true, credit=false) -# Literate.notebook(joinpath(root, file), out_dir; -# execute=true, documenter=true, credit=false) -# end -# end -# end +@info "Building literate files" +for (root, dirs, files) in walkdir(literate_dir) + out_dir = joinpath(generated_dir, relpath(root, literate_dir)) + pmap(files) do file + f,l = splitext(file) + if l == ".jl" && !startswith(f, "_") + Literate.markdown(joinpath(root, file), out_dir; + config=config, documenter=true, credit=false) + Literate.notebook(joinpath(root, file), out_dir; + execute=true, documenter=true, credit=false) + end + end +end +@info "Completed literate" + +pages = Any[] +push!(pages, "Decapodes.jl" => "index.md") +push!(pages, "Overview" => "overview.md") +push!(pages, "Equations" => "equations.md") +push!(pages, "BC Debug" => "bc_debug.md") +push!(pages, "ASCII Operators" => "ascii.md") +dirs = Dict("physics" => "Physics" + ,"biology" => "Biology" + ,"climate" => "Climate") +for d in keys(dirs) + dir = joinpath(@__DIR__, "src", d) + files = readdir(dir) + push!(pages, dirs[d] => joinpath.(d, files)) +end + +push!(pages, "Examples" => [ + "Brusselator" => "examples/chemistry/brusselator.md", + "Brusselator Teapot" => "examples/chemistry/brusselator_teapot.md", + "Gray-Scott" => "examples/chemistry/gray_scott.md", + "Budyko-Sellers" => "examples/climate/budyko_sellers.md", + "Budyko-Sellers-Halfar" => "examples/climate/budyko_sellers_halfar.md", + "Burger" => "examples/diff_adv/burger.md", +]) +push!(pages, "Canonical Models" => "canon.md") +push!(pages, "Library Reference" => "api.md") @info "Building Documenter.jl docs" makedocs( @@ -40,30 +66,12 @@ makedocs( format = Documenter.HTML( assets = ["assets/analytics.js"], ), + remotes = nothing, sitename = "Decapodes.jl", doctest = false, checkdocs = :none, - pages = Any[ - "Decapodes.jl" => "index.md", - "Halfar-NS" => "halmo.md", - "Overview" => "overview.md", - "Klausmeier" => "klausmeier.md", - "Glacial Flow" => "ice_dynamics.md", - "Grigoriev Ice Cap" => "grigoriev.md", - "Budyko-Sellers-Halfar" => "budyko_sellers_halfar.md", - "CISM v2.1" => "cism.md", - "NHS" => "nhs.md", - "Equations" => "equations.md", - "ASCII Operators" => "ascii.md", - "Misc Features" => "bc_debug.md", - "Pipe Flow" => "poiseuille.md", - # "Examples" => Any[ - # "examples/cfd_example.md" - # ], - "Canonical Models" => "canon.md", - "Library Reference" => "api.md" - ] -) + pages = pages) + @info "Deploying docs" deploydocs( diff --git a/docs/src/klausmeier.md b/docs/src/biology/klausmeier.md similarity index 100% rename from docs/src/klausmeier.md rename to docs/src/biology/klausmeier.md diff --git a/docs/src/budyko_sellers_halfar.md b/docs/src/climate/budyko_sellers_halfar.md similarity index 100% rename from docs/src/budyko_sellers_halfar.md rename to docs/src/climate/budyko_sellers_halfar.md diff --git a/docs/src/cism.md b/docs/src/climate/cism.md similarity index 100% rename from docs/src/cism.md rename to docs/src/climate/cism.md diff --git a/docs/src/grigoriev.md b/docs/src/climate/grigoriev.md similarity index 91% rename from docs/src/grigoriev.md rename to docs/src/climate/grigoriev.md index 736f0745..527aaffa 100644 --- a/docs/src/grigoriev.md +++ b/docs/src/climate/grigoriev.md @@ -34,7 +34,7 @@ Point3D = Point3{Float64}; # hide The ice thickness data is [stored in a TIF](https://zenodo.org/api/records/7735970/files-archive). We have downloaded it locally, and load it using basic `FileIO`. ``` @example DEC -file_name = "Icethickness_Grigoriev_ice_cap_2021.tif" +file_name = "../../assets/Icethickness_Grigoriev_ice_cap_2021.tif" ice_thickness_tif = load(file_name) ``` @@ -112,7 +112,7 @@ halfar_eq2 = @decapode begin n::Constant ḣ == ∂ₜ(h) - ḣ == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n-1)) * avg₀₁(h^(n+2))) + ḣ == ∘(⋆, d, ⋆)(Γ * d(h) ∧ (mag(♯(d(h)))^(n-1)) ∧ (h^(n+2))) end glens_law = @decapode begin @@ -138,36 +138,15 @@ to_graphviz(ice_dynamics) # Define our functions ``` @example DEC -include("sharp_op.jl") function generate(sd, my_symbol; hodge=GeometricHodge()) - ♯_m = ♯_mat(sd) - I = Vector{Int64}() - J = Vector{Int64}() - V = Vector{Float64}() - for e in 1:ne(s) - append!(J, [s[e,:∂v0],s[e,:∂v1]]) - append!(I, [e,e]) - append!(V, [0.5, 0.5]) - end - avg_mat = sparse(I,J,V) op = @match my_symbol begin - :♯ => x -> begin - ♯(sd, EForm(x)) - end :mag => x -> begin norm.(x) end - :avg₀₁ => x -> begin - avg_mat * x - end - :^ => (x,y) -> x .^ y - :* => (x,y) -> x .* y - :abs => x -> abs.(x) - :show => x -> begin - println(x) - x + :^ => (x,y) -> begin + x .^ y end - x => error("Unmatched operator $my_symbol") + _ => default_dec_matrix_generate(sd, my_symbol, hodge) end return (args...) -> op(args...) end diff --git a/docs/src/halmo.md b/docs/src/climate/halmo.md similarity index 100% rename from docs/src/halmo.md rename to docs/src/climate/halmo.md diff --git a/docs/src/ice_dynamics.md b/docs/src/climate/ice_dynamics.md similarity index 100% rename from docs/src/ice_dynamics.md rename to docs/src/climate/ice_dynamics.md diff --git a/docs/src/nhs.md b/docs/src/climate/nhs.md similarity index 100% rename from docs/src/nhs.md rename to docs/src/climate/nhs.md diff --git a/docs/src/poiseuille.md b/docs/src/physics/poiseuille.md similarity index 99% rename from docs/src/poiseuille.md rename to docs/src/physics/poiseuille.md index 53997083..d7e80bfc 100644 --- a/docs/src/poiseuille.md +++ b/docs/src/physics/poiseuille.md @@ -55,7 +55,7 @@ In order to solve our equations, we will need numerical linear operators that gi ```@example Poiseuille using MLStyle -include("../../examples/boundary_helpers.jl") +include("../../../examples/boundary_helpers.jl") function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin diff --git a/docs/src/sharp_op.jl b/docs/src/sharp_op.jl deleted file mode 100644 index 7da5a2d9..00000000 --- a/docs/src/sharp_op.jl +++ /dev/null @@ -1,60 +0,0 @@ -# TODO: Upstream these operators to CombinatorialSpaces.jl: -# https://github.com/AlgebraicJulia/CombinatorialSpaces.jl/pull/59/files - -using SparseArrays -using StaticArrays - -#""" Divided weighted normals by | σⁿ | . -#This weighting is that used in equation 5.8.1 from Hirani. -#See Hirani §5.8. -#""" -#♯_denominator(s::AbstractDeltaDualComplex2D, _::Int, t::Int) = -# volume(2,s,t) - -""" Divided weighted normals by | ⋆v | . -This weighting is NOT that of equation 5.8.1, but a different weighting scheme. -We essentially replace the denominator in equation 5.8.1 with | ⋆v | . This -may be what Hirani intended, and perhaps the denominator | σⁿ | in that equation -is either a mistake or clerical error. -See Hirani §5.8. -""" -♯_denominator(s::AbstractDeltaDualComplex2D, v::Int, _::Int) = - sum(dual_volume(2,s, elementary_duals(0,s,v))) - -""" Find a vector orthogonal to e pointing into the triangle shared with v. -""" -function get_orthogonal_vector(s::AbstractDeltaDualComplex2D, v::Int, e::Int) - e_vec = point(s, tgt(s, e)) - point(s, src(s, e)) - e_vec /= norm(e_vec) - e2_vec = point(s, v) - point(s, src(s, e)) - e2_vec - dot(e2_vec, e_vec)*e_vec -end - -function ♯_assign!(♯_mat::AbstractSparseMatrix, s::AbstractDeltaDualComplex2D, - v₀::Int, _::Int, t::Int, i::Int, tri_edges::SVector{3, Int}, tri_center::Int, - out_vec) - for e in deleteat(tri_edges, i) - v, sgn = src(s,e) == v₀ ? (tgt(s,e), -1) : (src(s,e), +1) - # | ⋆vₓ ∩ σⁿ | - dual_area = sum(dual_volume(2,s,d) for d in elementary_duals(0,s,v) - if s[s[d, :D_∂e0], :D_∂v0] == tri_center) - area = ♯_denominator(s, v, t) - ♯_mat[v,e] += sgn * sign(1,s,e) * (dual_area / area) * out_vec - end -end - -function ♯_mat(s::AbstractDeltaDualComplex2D) - #♯_mat = spzeros(attrtype_type(s, :Point), (nv(s), ne(s))) - ♯_mat = spzeros(Point3D, (nv(s), ne(s))) - for t in triangles(s) - tri_center, tri_edges = triangle_center(s,t), triangle_edges(s,t) - for (i, (v₀, e₀)) in enumerate(zip(triangle_vertices(s,t), tri_edges)) - out_vec = get_orthogonal_vector(s, v₀, e₀) - h = norm(out_vec) - #out_vec /= DS == DesbrunSharp() ? h : h^2 - out_vec /= h^2 - ♯_assign!(♯_mat, s, v₀, e₀, t, i, tri_edges, tri_center, out_vec) - end - end - ♯_mat -end diff --git a/examples/chemistry/_brusselator.jl b/examples/chemistry/_brusselator.jl new file mode 100644 index 00000000..bab7f0c6 --- /dev/null +++ b/examples/chemistry/_brusselator.jl @@ -0,0 +1,177 @@ +using CombinatorialSpaces +using Decapodes +using DiagrammaticEquations +using MLStyle +using OrdinaryDiffEq +using LinearAlgebra +using ComponentArrays +using CairoMakie +using GeometryBasics: Point2, Point3 +Point2D = Point2{Float64} +Point3D = Point3{Float64} + +Brusselator = @decapode begin + ## Values living on vertices. + (U, V)::Form0{X} # State variables. + (U2V)::Form0{X} # Named intermediate variables. + (U̇, V̇)::Form0{X} # Tangent variables. + ## Scalars. + (α)::Constant{X} + F::Parameter{X} + ## A named intermediate variable. + U2V == (U .* U) .* V + ## Specify how to compute the tangent variables. + U̇ == 1 + U2V - (4.4 * U) + (α * Δ(U)) + F + V̇ == (3.4 * U) - U2V + (α * Δ(U)) + ## Associate tangent variables with a state variable. + ∂ₜ(U) == U̇ + ∂ₜ(V) == V̇ +end + +# TODO: Create square domain of approximately 32x32 vertices. +# s = loadmesh(Rectangle_30x10()) +# Visualize the mesh. +s = triangulated_grid(200, 200, 1, 1, Point3D); +scaling_mat = Diagonal([1/maximum(x->x[1], s[:point]), + 1/maximum(x->x[2], s[:point]), + 1.0]); +s[:point] = map(x -> scaling_mat*x, s[:point]); +sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s); +subdivide_duals!(sd, Circumcenter()) + +# Define how operations map to Julia functions. +function generate(sd, my_symbol; hodge=DiagonalHodge()) + op = @match my_symbol begin + x => error("Unmatched operator $my_symbol") + end + return (args...) -> op(args...) +end + +U = map(sd[:point]) do (_,y) + 22 * (y *(1-y))^(3/2) +end + +V = map(sd[:point]) do (x,_) + 27 * (x *(1-x))^(3/2) +end + +F₁ = map(sd[:point]) do (x,y) + (x-0.3)^2 + (y-0.6)^2 ≤ (0.1)^2 ? 5.0 : 0.0 +end + +F₂ = zeros(nv(sd)) + +constants_and_parameters = ( + α = 0.001, + F = t -> t ≥ 1.1 ? F₂ : F₁) + +gensim(Brusselator) +sim = evalsim(Brusselator) +fₘ = sim(sd, generate) + +# Create problem and run sim for t ∈ [0,tₑ). +# Map symbols to data. +u₀ = ComponentArray(U=U, V=V) + +tₑ = 25 + +begin + @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()) + @info("Done") +end + +begin + frames = 300 + fig = Figure() + ax = CairoMakie.Axis(fig[1,1]) + msh = CairoMakie.mesh!(ax, s, color=soln(0).U, colormap=:jet, colorrange=extrema(soln(0).U)) + Colorbar(fig[1,2], msh) + CairoMakie.record(fig, "Brusselator_CPU.gif", range(0.0, tₑ; length=frames); framerate = 30) do t + msh.color = soln(t).U + end +end + +# Run on the sphere. +# You can use lower resolution meshes, such as Icosphere(3). +s = loadmesh(Icosphere(5)); + +# Visualize the mesh. +wireframe(s) + +sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s); +subdivide_duals!(sd, Circumcenter()) + +# Create initial data. +U = map(sd[:point]) do (_,y,_) + abs(y) +end + +V = map(sd[:point]) do (x,_,_) + abs(x) +end + +# TODO: Try making this sparse. +F₁ = map(sd[:point]) do (_,_,z) + z ≥ 0.8 ? 5.0 : 0.0 +end +mesh(s, color=F₁, colormap=:jet) + +# TODO: Try making this sparse. +F₂ = zeros(nv(sd)) + +constants_and_parameters = ( + α = 0.001, + F = t -> t ≥ 1.1 ? F₁ : F₂ + ) + +# Generate the simulation. +fₘ = sim(sd, generate) + +# Create problem and run sim for t ∈ [0,tₑ). +# Map symbols to data. +u₀ = ComponentArray(U=U, V=V) + +# Visualize the initial conditions. +# If GLMakie throws errors, then update your graphics drivers, +# or use an alternative Makie backend like CairoMakie. +fig_ic = Figure() +p1 = mesh(fig_ic[1,2], s, color=u₀.U, colormap=:jet) +p2 = mesh(fig_ic[1,3], s, color=u₀.V, colormap=:jet) +display(fig_ic) + +tₑ = 11.5 + +@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()) +@info("Done") + +# Visualize the final conditions. +mesh(s, color=soln(tₑ).U, colormap=:jet) + +# BEGIN Gif creation +begin +frames = 800 +fig = Figure(resolution = (1200, 1200)) +p1 = mesh(fig[1,1], s, color=soln(0).U, colormap=:jet, colorrange=extrema(soln(0).U)) +p2 = mesh(fig[2,1], s, color=soln(0).V, colormap=:jet, colorrange=extrema(soln(0).V)) +Colorbar(fig[1,2]) +Colorbar(fig[2,2]) + +record(fig, "brusselator_sphere.gif", range(0.0, tₑ; length=frames); framerate = 30) do t + p1.plot.color = soln(t).U + p2.plot.color = soln(t).V +end + +end +# END Gif creation diff --git a/examples/chemistry/brusselator_bounded.jl b/examples/chemistry/_brusselator_bounded.jl similarity index 99% rename from examples/chemistry/brusselator_bounded.jl rename to examples/chemistry/_brusselator_bounded.jl index 44f9405f..32400166 100644 --- a/examples/chemistry/brusselator_bounded.jl +++ b/examples/chemistry/_brusselator_bounded.jl @@ -9,7 +9,6 @@ using ComponentArrays using DiagrammaticEquations using DiagrammaticEquations.Deca using Decapodes -using MultiScaleArrays using MLStyle using OrdinaryDiffEq using LinearAlgebra diff --git a/examples/chemistry/brusselator.jl b/examples/chemistry/brusselator.jl deleted file mode 100644 index ba8a70a9..00000000 --- a/examples/chemistry/brusselator.jl +++ /dev/null @@ -1,267 +0,0 @@ -using Catlab -using Catlab.Graphics -using CombinatorialSpaces -using CombinatorialSpaces.ExteriorCalculus -using DiagrammaticEquations -using DiagrammaticEquations.Deca -using Decapodes -using MLStyle -using OrdinaryDiffEq -using LinearAlgebra -using Logging -using JLD2 -using Printf -using CairoMakie -import CairoMakie: wireframe, mesh, Figure, Axis -using ComponentArrays - - -using GeometryBasics: Point2, Point3 -Point2D = Point2{Float64} -Point3D = Point3{Float64} - -# Values living on vertices. -# State variables. -# Named intermediate variables. -Brusselator = @decapode begin - (U, V)::Form0{X} - (U2V, One)::Form0{X} - (U̇, V̇)::Form0{X} - - (α)::Constant{X} - F::Parameter{X} - - U2V == (U .* U) .* V - - U̇ == 1 + U2V - (4.4 * U) + (α * Δ(U)) + F - V̇ == (3.4 * U) - U2V + (α * Δ(U)) - ∂ₜ(U) == U̇ - ∂ₜ(V) == V̇ -end -# Visualize. You must have graphviz installed. -to_graphviz(Brusselator) - -# We resolve types of intermediate variables using sets of rules. - -infer_types!(Brusselator) -# Visualize. Note that variables now all have types. -to_graphviz(Brusselator) - -# Resolve overloads. i.e. ~dispatch -resolve_overloads!(Brusselator) -# Visualize. Note that functions are renamed. -to_graphviz(Brusselator) - -# TODO: Create square domain of approximately 32x32 vertices. -s = loadmesh(Rectangle_30x10()) -scaling_mat = Diagonal([1/maximum(x->x[1], s[:point]), - 1/maximum(x->x[2], s[:point]), - 1.0]) -s[:point] = map(x -> scaling_mat*x, s[:point]) -s[:edge_orientation] = false -orient!(s) -# Visualize the mesh. -wireframe(s) -sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s) -subdivide_duals!(sd, Circumcenter()) - -# Define how operations map to Julia functions. -hodge = GeometricHodge() -Δ₀ = δ(1, sd, hodge=hodge) * d(0, sd) -function generate(sd, my_symbol; hodge=GeometricHodge()) - op = @match my_symbol begin - - :Δ₀ => x -> Δ₀ * x - :.* => (x,y) -> x .* y - x => error("Unmatched operator $my_symbol") - end - return (args...) -> op(args...) -end - -# Create initial data. -@assert all(map(sd[:point]) do (x,y) - 0.0 ≤ x ≤ 1.0 && 0.0 ≤ y ≤ 1.0 -end) - -U = map(sd[:point]) do (_,y) - 22 * (y *(1-y))^(3/2) -end - -V = map(sd[:point]) do (x,_) - 27 * (x *(1-x))^(3/2) -end - -# TODO: Try making this sparse. -F₁ = map(sd[:point]) do (x,y) - (x-0.3)^2 + (y-0.6)^2 ≤ (0.1)^2 ? 5.0 : 0.0 -end -mesh(s, color=F₁, colormap=:jet) - -# TODO: Try making this sparse. -F₂ = zeros(nv(sd)) - -One = ones(nv(sd)) - -constants_and_parameters = ( - fourfour = 4.4, - threefour = 3.4, - α = 0.001, - F = t -> t ≥ 1.1 ? F₂ : F₁) - -# Generate the simulation. -gensim(expand_operators(Brusselator)) -sim = eval(gensim(expand_operators(Brusselator))) -fₘ = sim(sd, generate) - -# Create problem and run sim for t ∈ [0,tₑ). -# Map symbols to data. -u₀ = ComponentArray(U=U, V=V, One=One) - -# Visualize the initial conditions. -# If GLMakie throws errors, then update your graphics drivers, -# or use an alternative Makie backend like -fig_ic = Figure() -p1 = mesh(fig_ic[1,2], s, color=u₀.U, colormap=:jet) -p2 = mesh(fig_ic[1,3], s, color=u₀.V, colormap=:jet) - -tₑ = 11.5 - -@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()) -@info("Done") - -@save "brusselator.jld2" soln - -# Visualize the final conditions. -mesh(s, color=soln(tₑ).U, colormap=:jet) - -# BEGIN Gif creation -begin -frames = 100 -fig = Figure(resolution = (1200, 800)) -p1 = mesh(fig[1,2], s, color=soln(0).U, colormap=:jet, colorrange=extrema(soln(0).U)) -p2 = mesh(fig[1,4], s, color=soln(0).V, colormap=:jet, colorrange=extrema(soln(0).V)) -ax1 = Axis(fig[1,2], width = 400, height = 400) -ax2 = Axis(fig[1,4], width = 400, height = 400) -hidedecorations!(ax1) -hidedecorations!(ax2) -hidespines!(ax1) -hidespines!(ax2) -Colorbar(fig[1,1]) -Colorbar(fig[1,5]) -Label(fig[1,2,Top()], "U") -Label(fig[1,4,Top()], "V") -lab1 = Label(fig[1,3], "") - -## Animation -record(fig, "brusselator.gif", range(0.0, tₑ; length=frames); framerate = 15) do t - p1.plot.color = soln(t).U - p2.plot.color = soln(t).V - lab1.text = @sprintf("%.2f",t) -end - -end -# END Gif creation - -# Run on the sphere. -# You can use lower resolution meshes, such as Icosphere(3). -s = loadmesh(Icosphere(5)) -s[:edge_orientation] = false -orient!(s) -# Visualize the mesh. -wireframe(s) -sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s) -subdivide_duals!(sd, Circumcenter()) - -# Define how operations map to Julia functions. -hodge = GeometricHodge() -Δ₀ = δ(1, sd, hodge=hodge) * d(0, sd) -function generate(sd, my_symbol; hodge=GeometricHodge()) - op = @match my_symbol begin - :Δ₀ => x -> Δ₀ * x - :.* => (x,y) -> x .* y - x => error("Unmatched operator $my_symbol") - end - return (args...) -> op(args...) -end - -# Create initial data. -U = map(sd[:point]) do (_,y,_) - abs(y) -end - -V = map(sd[:point]) do (x,_,_) - abs(x) -end - -# TODO: Try making this sparse. -F₁ = map(sd[:point]) do (_,_,z) - z ≥ 0.8 ? 5.0 : 0.0 -end -mesh(s, color=F₁, colormap=:jet) - -# TODO: Try making this sparse. -F₂ = zeros(nv(sd)) - -One = ones(nv(sd)) - -constants_and_parameters = ( - fourfour = 4.4, - threefour = 3.4, - α = 0.001, - F = t -> t ≥ 1.1 ? F₁ : F₂ - ) - -# Generate the simulation. -fₘ = sim(sd, generate) - -# Create problem and run sim for t ∈ [0,tₑ). -# Map symbols to data. -u₀ = ComponentArray(U=U, V=V, One=One) - -# Visualize the initial conditions. -# If GLMakie throws errors, then update your graphics drivers, -# or use an alternative Makie backend like CairoMakie. -fig_ic = Figure() -p1 = mesh(fig_ic[1,2], s, color=u₀.U, colormap=:jet) -p2 = mesh(fig_ic[1,3], s, color=u₀.V, colormap=:jet) -display(fig_ic) - -tₑ = 11.5 - -@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()) -@info("Done") - -@save "brusselator_sphere.jld2" soln - -# Visualize the final conditions. -mesh(s, color=soln(tₑ).U, colormap=:jet) - -# BEGIN Gif creation -begin -frames = 800 -fig = Figure(resolution = (1200, 1200)) -p1 = mesh(fig[1,1], s, color=soln(0).U, colormap=:jet, colorrange=extrema(soln(0).U)) -p2 = mesh(fig[2,1], s, color=soln(0).V, colormap=:jet, colorrange=extrema(soln(0).V)) -Colorbar(fig[1,2]) -Colorbar(fig[2,2]) - -## Animation -record(fig, "brusselator_sphere.gif", range(0.0, tₑ; length=frames); framerate = 30) do t - p1.plot.color = soln(t).U - p2.plot.color = soln(t).V -end - -end -# END Gif creation diff --git a/examples/climate/budyko_sellers.jl b/examples/climate/budyko_sellers.jl index 4cc77cc2..79f53d73 100644 --- a/examples/climate/budyko_sellers.jl +++ b/examples/climate/budyko_sellers.jl @@ -16,7 +16,7 @@ using OrdinaryDiffEq using JLD2 # Uncomment to load GLMakie if your system supports it. # Otherwise, do using CairoMakie -#using GLMakie +using CairoMakie using GeometryBasics: Point2 Point2D = Point2{Float64} @@ -141,9 +141,6 @@ C = map(point(s′)) do ϕ end D = 0.6 -#Tₛ₀ = map(point(s′)) do ϕ -# 12 .- 40*((1/2)*(3*(sin(ϕ[1]))^2 - 1)) -#end # Isothermal initial conditions: Tₛ₀ = map(point(s′)) do ϕ 15 diff --git a/examples/climate/budyko_sellers_halfar.jl b/examples/climate/budyko_sellers_halfar.jl index b40649f2..6f9a8e06 100644 --- a/examples/climate/budyko_sellers_halfar.jl +++ b/examples/climate/budyko_sellers_halfar.jl @@ -1,5 +1,3 @@ -# Import Dependencies - # AlgebraicJulia Dependencies using Catlab using Catlab.Graphics @@ -7,39 +5,107 @@ using CombinatorialSpaces using DiagrammaticEquations using DiagrammaticEquations.Deca using Decapodes -using Decapodes: SchSummationDecapode # External Dependencies using MLStyle +using ComponentArrays using LinearAlgebra using OrdinaryDiffEq using JLD2 using CairoMakie +using SparseArrays using GeometryBasics: Point2 -using ComponentArray -Point2D = Point2{Float64} +Point2D = Point2{Float64}; -# Import prior models -## If the Budyko-Sellers and Halfar models are not already created, they can be -## with these scripts: -## include("budyko_sellers.jl") -## include("shallow_ice.jl") +halfar_eq2 = @decapode begin + h::Form0 + Γ::Form1 + n::Constant -budyko_sellers = apex(budyko_sellers_cospan) + ḣ == ∂ₜ(h) + ḣ == ∘(⋆, d, ⋆)(Γ * d(h) ∧ (mag(♯(d(h)))^(n-1)) ∧ (h^(n+2))) +end +glens_law = @decapode begin + Γ::Form1 + 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_cospan = oapply(ice_dynamics_composition_diagram, + [Open(halfar_eq2, [:Γ,:n]), + Open(glens_law, [:Γ,:n])]) halfar = apex(ice_dynamics_cospan) +to_graphviz(halfar, verbose=false) + +energy_balance = @decapode begin + (Tₛ, ASR, OLR, HT)::Form0 + (C)::Constant + + Tₛ̇ == ∂ₜ(Tₛ) + + 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 -to_graphviz(budyko_sellers) -to_graphviz(halfar) + Q == 450 * cosϕᵖ +end + +to_graphviz(oplus([energy_balance, absorbed_shortwave_radiation, outgoing_longwave_radiation, heat_transfer, insolation]), directed=false) + +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 +to_graphviz(budyko_sellers_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") + +budyko_sellers_cospan = 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ϕᵖ])]) + +budyko_sellers = apex(budyko_sellers_cospan) +write_json_acset(budyko_sellers, "budyko_sellers.json") # Save this Decapode as a JSON file +to_graphviz(budyko_sellers, verbose=false) -## Define the model -## Tₛ(ϕ,t) := Surface temperature -## A(ϕ) := Longwave emissions at 0°C warming = @decapode begin - (Tₛ)::Form0 - (A)::Form1 + Tₛ::Form0 + A::Form1 - ## A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e-17) A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e7) end @@ -56,34 +122,24 @@ to_graphviz(budyko_sellers_halfar_composition_diagram, box_labels=:name, junctio budyko_sellers_halfar_cospan = oapply(budyko_sellers_halfar_composition_diagram, [Open(budyko_sellers, [:Tₛ]), - Open(warming, [:A, :Tₛ]), - Open(halfar, [:stress_A])]) - + Open(warming, [:A, :Tₛ]), + Open(halfar, [:stress_A])]) budyko_sellers_halfar = apex(budyko_sellers_halfar_cospan) to_graphviz(budyko_sellers_halfar) budyko_sellers_halfar = expand_operators(budyko_sellers_halfar) infer_types!(budyko_sellers_halfar, op1_inf_rules_1D, op2_inf_rules_1D) -to_graphviz(budyko_sellers_halfar) - resolve_overloads!(budyko_sellers_halfar, op1_res_rules_1D, op2_res_rules_1D) to_graphviz(budyko_sellers_halfar) -## Define the mesh - s′ = EmbeddedDeltaSet1D{Bool, Point2D}() -## add_vertices!(s′, 30, point=Point2D.(range(-π/2 + π/32, π/2 - π/32, length=30), 0)) add_vertices!(s′, 100, point=Point2D.(range(-π/2 + π/32, π/2 - π/32, length=100), 0)) add_edges!(s′, 1:nv(s′)-1, 2:nv(s′)) orient!(s′) s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s′) subdivide_duals!(s, Circumcenter()) -# 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 @@ -104,24 +160,24 @@ C = map(point(s′)) do ϕ end D = 0.6 -# Isothermal initial conditions: Tₛ₀ = map(point(s′)) do ϕ 15 end +lines(map(x -> x[1], point(s′)), Tₛ₀) + n = 3 ρ = 910 g = 9.8 -# Ice height is a primal 0-form, with values at vertices. h₀ = map(point(s′)) do (x,_) (((x)^2)+2.5) / 1e3 end -# Visualize initial condition for ice sheet height. + lines(map(x -> x[1], point(s′)), h₀) # Store these values to be passed to the solver. -u₀ = ComponentArray(Tₛ=Tₛ₀, halfar_h=h₀) +u₀ = ComponentArray(Tₛ=Tₛ₀, halfar_dynamics_h=h₀) constants_and_parameters = ( budyko_sellers_absorbed_radiation_α = α, @@ -135,15 +191,11 @@ constants_and_parameters = ( halfar_stress_ρ = ρ, halfar_stress_g = g) -# Define how symbols map to Julia functions +# Symbols to functions function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin :♯ => x -> begin - ## This is an implementation of the "sharp" operator from the exterior - ## calculus, which takes co-vector fields to vector fields. - ## This could be up-streamed to the CombinatorialSpaces.jl library. (i.e. - ## this operation is not bespoke to this simulation.) e_vecs = map(edges(sd)) do e point(sd, sd[e, :∂v0]) - point(sd, sd[e, :∂v1]) end @@ -174,80 +226,64 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) end :^ => (x,y) -> x .^ y :* => (x,y) -> x .* y - :show => x -> begin - @show x - x - end x => error("Unmatched operator $my_symbol") end return (args...) -> op(args...) end -## Generate simulation +## Simulation generation sim = eval(gensim(budyko_sellers_halfar, dimension=1)) fₘ = sim(s, generate) -## Run simulation - tₑ = 1e6 -#tₑ = 5e13 # Julia will pre-compile 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 = solve(prob, Tsit5()) soln.retcode != :Unstable || error("Solver was not stable") -# This next run should be fast. @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.jld2" soln -# Visualize +## Visualize -lines(map(x -> x[1], point(s′)), soln(0.0).Tₛ) lines(map(x -> x[1], point(s′)), soln(tₑ).Tₛ) -lines(map(x -> x[1], point(s′)), soln(0.0).halfar_h) -lines(map(x -> x[1], point(s′)), soln(tₑ).halfar_h) +lines(map(x -> x[1], point(s′)), soln(tₑ).halfar_dynamics_h) begin -## Initial frame frames = 100 -fig = Figure(resolution = (800, 800)) +fig = Figure(; size = (800, 800)) ax1 = CairoMakie.Axis(fig[1,1]) xlims!(ax1, extrema(map(x -> x[1], point(s′)))) ylims!(ax1, extrema(soln(tₑ).Tₛ)) Label(fig[1,1,Top()], "Surface temperature, Tₛ, [C°]") Label(fig[2,1,Top()], "Line plot of temperature from North to South pole, every $(tₑ/frames) time units") -## Animation record(fig, "budyko_sellers_halfar_T.gif", range(0.0, tₑ; length=frames); framerate = 15) do t lines!(fig[1,1], map(x -> x[1], point(s′)), soln(t).Tₛ) end end begin -## Initial frame frames = 100 -fig = Figure(resolution = (800, 800)) +fig = Figure(; size = (800, 800)) ax1 = CairoMakie.Axis(fig[1,1]) xlims!(ax1, extrema(map(x -> x[1], point(s′)))) -ylims!(ax1, extrema(soln(tₑ).halfar_h)) +ylims!(ax1, extrema(soln(tₑ).halfar_dynamics_h)) Label(fig[1,1,Top()], "Ice height, h") Label(fig[2,1,Top()], "Line plot of ice height from North to South pole, every $(tₑ/frames) time units") -## Animation record(fig, "budyko_sellers_halfar_h.gif", range(0.0, tₑ; length=frames); framerate = 15) do t - lines!(fig[1,1], map(x -> x[1], point(s′)), soln(t).halfar_h) + lines!(fig[1,1], map(x -> x[1], point(s′)), soln(t).halfar_dynamics_h) end end + + diff --git a/examples/diff_adv/Project.toml b/examples/diff_adv/Project.toml deleted file mode 100644 index b080ae88..00000000 --- a/examples/diff_adv/Project.toml +++ /dev/null @@ -1,9 +0,0 @@ -[deps] -CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" -CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2" -Decapodes = "679ab3ea-c928-4fe6-8d59-fd451142d391" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" -JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" -OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" diff --git a/examples/diff_adv/diff_adv.jl b/examples/diff_adv/_diff_adv.jl similarity index 100% rename from examples/diff_adv/diff_adv.jl rename to examples/diff_adv/_diff_adv.jl diff --git a/examples/diff_adv/burger.jl b/examples/diff_adv/burger.jl index 86464ec4..d3eb43bd 100644 --- a/examples/diff_adv/burger.jl +++ b/examples/diff_adv/burger.jl @@ -12,11 +12,12 @@ using Logging: global_logger using GeometryBasics: Point2 Point2D = Point2{Float64} using Distributions -using CairoMakie using LinearAlgebra using MLStyle -using ComponentArrays using OrdinaryDiffEq +using CairoMakie +import CairoMakie: Axis +using ComponentArrays # Represent component Decapodes. Diffusion = @decapode begin diff --git a/examples/sw/pressure.jl b/examples/sw/_pressure.jl similarity index 100% rename from examples/sw/pressure.jl rename to examples/sw/_pressure.jl diff --git a/examples/sw/sw.jl b/examples/sw/_sw.jl similarity index 100% rename from examples/sw/sw.jl rename to examples/sw/_sw.jl diff --git a/examples/sw/sw_with_advection.jl b/examples/sw/_sw_with_advection.jl similarity index 100% rename from examples/sw/sw_with_advection.jl rename to examples/sw/_sw_with_advection.jl