diff --git a/.gitignore b/.gitignore index 0d9d8cf9..615d61e1 100644 --- a/.gitignore +++ b/.gitignore @@ -6,3 +6,4 @@ .vscode .buildkite/toolbox/config.json .buildkite/toolbox/username +docs/Manifest.toml diff --git a/docs/Project.toml b/docs/Project.toml index a1cacd15..8f660d9f 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,4 @@ [deps] -AlgebraicPetri = "4f99eebe-17bf-4e98-b6a1-2c4f205a959b" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" @@ -10,10 +9,9 @@ DiagrammaticEquations = "6f00c28b-6bed-4403-80fa-30e0dc12f317" Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" +Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6" FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" -HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" @@ -21,5 +19,5 @@ Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d" -ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/docs/docinfo.jl b/docs/docinfo.jl new file mode 100644 index 00000000..9f61e9f7 --- /dev/null +++ b/docs/docinfo.jl @@ -0,0 +1,26 @@ +module DocInfo + +using Dates + +mutable struct Info + start_time::DateTime + finish_time::DateTime +end + +function Info() + Info(now(), now()) +end + +function get_elapsed(info::Info) + info.finish_time - info.start_time +end + +function get_report(info::Info) + info.finish_time = now() + elapsed = get_elapsed(info) + elapsed_sec = round(elapsed, Dates.Second(1)) + @info "Page built in $(elapsed_sec)." + @info "This page was last built at $(info.finish_time)." +end + +end diff --git a/docs/make.jl b/docs/make.jl index 59ded021..5968bc13 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,13 +1,10 @@ +@info "Loading Documenter" using Documenter using Literate using Distributed @info "Loading Decapodes" using Decapodes -using Catlab -using Catlab.WiringDiagrams -using AlgebraicPetri -using CairoMakie # Set Literate.jl config if not being compiled on recognized service. config = Dict{String,String}() @@ -43,24 +40,25 @@ makedocs( sitename = "Decapodes.jl", doctest = false, checkdocs = :none, + pagesonly = true, + linkcheck = true, + linkcheck_ignore = [r"agupubs.onlinelibrary.wiley.com"], # This gives a 403 Forbidden pages = Any[ "Decapodes.jl" => "index.md", + "Overview" => "overview/overview.md", + "Equations" => "equations/equations.md", "Vortices" => "navier_stokes/ns.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", + "Cahn-Hilliard" => "ch/cahn-hilliard.md", + "Klausmeier" => "klausmeier/klausmeier.md", + "CISM v2.1" => "cism/cism.md", + "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-NS" => "halmo/halmo.md", # Requires grigoriev + "NHS" => "nhs/nhs_lite.md", + "Pipe Flow" => "poiseuille/poiseuille.md", + "Misc Features" => "bc/bc_debug.md", # Requires overview "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" ] diff --git a/docs/src/api.md b/docs/src/api.md index d3b4d3c9..21b0ba3f 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -1,7 +1,17 @@ # Library Reference +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "docinfo.jl")) +info = DocInfo.Info() +``` + ## Decapodes + ```@autodocs Modules = [ Decapodes ] Private = false ``` + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/ascii.md b/docs/src/ascii.md index 7f9c06b8..8972f4c0 100644 --- a/docs/src/ascii.md +++ b/docs/src/ascii.md @@ -1,5 +1,10 @@ # ASCII and Vector Calculus Operators +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "docinfo.jl")) +info = DocInfo.Info() +``` + Some users may have trouble entering unicode characters like ⋆ or ∂ in their development environment. So, we offer the following ASCII equivalents. Further, some users may like to use vector calculus symbols instead of exterior calculus symbols where possible. We offer support for such symbols as well. ## ASCII Equivalents @@ -25,3 +30,6 @@ Some users may have trouble entering unicode characters like ⋆ or ∂ in their | ∇x | ∘(d,⋆) | \nabla \ x | | adv(X,Y) | ∘(⋆,d,⋆)(X∧Y) | adv | +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/bc_debug.md b/docs/src/bc/bc_debug.md similarity index 81% rename from docs/src/bc_debug.md rename to docs/src/bc/bc_debug.md index 40925f16..022c380a 100644 --- a/docs/src/bc_debug.md +++ b/docs/src/bc/bc_debug.md @@ -1,22 +1,23 @@ # Simulation Setup +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + This tutorial showcases some of the other features included in the Decapodes.jl package. Currently, these features are the treatment of boundary conditions and the simulation debugger interface. To begin, we set up the same -advection-diffusion problem presented in the Overview section. +advection-diffusion problem presented in the [Overview](../overview/overview.md) section. As before, we define the Diffusion, Advection, and Superposition components, -and now include a BC (Bounday Condition) component. Decapodes.jl interprets any -`Hom` which begins with a `∂` as a boundary condition. These boundary -conditions recieve special treatment at the scheduling step. Below we show the +and now include a Boundary Condition (BC) component. By convention, BCs are encoded in Decapodes by using a `∂` symbol. Below we show the graphical rendering of this boundary condition diagram, which we will use to -impose a Dirichlet condition on the time derivative of concentration at the +impose a [Dirichlet condition](https://en.wikipedia.org/wiki/Dirichlet_boundary_condition) on the time derivative of concentration at the mesh boundary. ```@example Debug using Catlab -using Catlab.Graphics using DiagrammaticEquations -using DiagrammaticEquations.Deca using Decapodes Diffusion = @decapode begin @@ -58,7 +59,6 @@ to_graphviz(BoundaryConditions) As before, we compose these physics components over our wiring diagram. - ```@example Debug compose_diff_adv = @relation (C, V) begin diffusion(C, ϕ₁) @@ -67,7 +67,7 @@ compose_diff_adv = @relation (C, V) begin superposition(ϕ₁, ϕ₂, ϕ, C_up, C) end -to_graphviz(compose_diff_adv, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(compose_diff_adv) ``` ```@example Debug @@ -87,11 +87,10 @@ ensures that this boundary condition holds true for any variables dependent on this variable, though also means that the boundary conditions on a variable have no immediate impact on the variables this variable is dependent on. -In the visualization below, wee see that the final operation +In the visualization below, we see that the final operation executed on the data is the boundary condition we are enforcing on the change in concentration. - ```@example Debug to_graphviz(DiffusionAdvection) ``` @@ -101,16 +100,15 @@ boundary conditions and so we will use the `plot_mesh` from the previous example instead of the mesh with periodic boundary conditions. Because the mesh is only a primal mesh, we also generate and subdivide the dual mesh. -`Rectangle_30x10` is a default mesh that is downloaded via `Artifacts.jl` when a user installs Decapodes. Via CombinatorialSpaces.jl, we can instantiate any `.obj` file of triangulated faces as a simplicial set. - ```@example Debug -using CombinatorialSpaces, CombinatorialSpaces.DiscreteExteriorCalculus +using CombinatorialSpaces using CairoMakie plot_mesh = loadmesh(Rectangle_30x10()) # Generate the dual mesh plot_mesh_dual = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3{Float64}}(plot_mesh) + # Calculate distances and subdivisions for the dual mesh subdivide_duals!(plot_mesh_dual, Circumcenter()) @@ -124,15 +122,14 @@ Finally, we define our operators, generate the simulation function, and compute the simulation. Note that when we define the boundary condition operator, we hardcode the boundary indices and values into the operator itself. We also move the initial concentration to the left, so that we are able to see a constant -concentration on the left boundary which will act as a source in the flow. The +concentration on the left boundary which will act as a source in the flow. You can find the file for boundary conditions [here](../boundary_helpers.jl). The modified initial condition is shown below: ```@example Debug using LinearAlgebra using ComponentArrays using MLStyle -using CombinatorialSpaces.DiscreteExteriorCalculus: ∧ -include("../../examples/boundary_helpers.jl") +include("../boundary_helpers.jl") function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin @@ -142,12 +139,9 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) x[boundary] .= 0 x end - :∧₀₁ => (x,y) -> begin - ∧(Tuple{0,1}, sd, x,y) - end x => error("Unmatched operator $my_symbol") end - return (args...) -> op(args...) + return op end using Distributions @@ -194,4 +188,8 @@ record(fig, "diff_adv_right.gif", range(0.0, 100.0; length=150); framerate = 30) end ``` -![](diff_adv_right.gif) +![Diffusion-Advection result and your first BC Decapode!](diff_adv_right.gif) + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/examples/boundary_helpers.jl b/docs/src/boundary_helpers.jl similarity index 100% rename from examples/boundary_helpers.jl rename to docs/src/boundary_helpers.jl diff --git a/docs/src/budyko_sellers_halfar.md b/docs/src/bsh/budyko_sellers_halfar.md similarity index 76% rename from docs/src/budyko_sellers_halfar.md rename to docs/src/bsh/budyko_sellers_halfar.md index 7768df8d..9aa58f67 100644 --- a/docs/src/budyko_sellers_halfar.md +++ b/docs/src/bsh/budyko_sellers_halfar.md @@ -1,29 +1,33 @@ # Budko-Sellers-Halfar +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..","docinfo.jl")) +info = DocInfo.Info() +``` + In this example, we will compose the Budyko-Sellers 1D energy balance model of the Earth's surface temperature with the Halfar model of glacial dynamics. Note that each of these components models is itself a composition of smaller physical models. In this walkthrough, we will compose them together using the same techniques. ``` @example DEC # AlgebraicJulia Dependencies using Catlab -using Catlab.Graphics using CombinatorialSpaces using DiagrammaticEquations -using DiagrammaticEquations.Deca using Decapodes # External Dependencies -using MLStyle +using CairoMakie using ComponentArrays +using GeometryBasics: Point2 +using JLD2 using LinearAlgebra +using MLStyle using OrdinaryDiffEq -using JLD2 -using CairoMakie using SparseArrays -using GeometryBasics: Point2 Point2D = Point2{Float64}; +nothing # hide ``` -We have defined the Halfar ice model in other docs pages, and so will quickly define it here. +We have defined the [Halfar ice model](../ice_dynamics/ice_dynamics.md) in other docs pages, and so will quickly define it here. ``` @example DEC halfar_eq2 = @decapode begin @@ -34,6 +38,7 @@ halfar_eq2 = @decapode begin ḣ == ∂ₜ(h) ḣ == ∘(⋆, d, ⋆)(Γ * d(h) ∧ (mag(♯(d(h)))^(n-1)) ∧ (h^(n+2))) end + glens_law = @decapode begin Γ::Form1 A::Form1 @@ -41,14 +46,17 @@ glens_law = @decapode begin Γ == (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) ``` @@ -105,7 +113,8 @@ budyko_sellers_composition_diagram = @relation () begin diffusion(Tₛ, HT, cosϕᵖ) insolation(Q, cosϕᵖ) end -to_graphviz(budyko_sellers_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") + +draw_composition(budyko_sellers_composition_diagram) ``` ``` @example DEC @@ -117,7 +126,10 @@ budyko_sellers_cospan = oapply(budyko_sellers_composition_diagram, 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 + +# Save this Decapode as a JSON file +write_json_acset(budyko_sellers, "budyko_sellers.json") + to_graphviz(budyko_sellers, verbose=false) ``` @@ -133,6 +145,7 @@ warming = @decapode begin A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e7) end + to_graphviz(warming) ``` @@ -143,12 +156,11 @@ Observe that Decapodes composition is hierarchical. This composition technique i ``` @example DEC budyko_sellers_halfar_composition_diagram = @relation () begin budyko_sellers(Tₛ) - warming(A, Tₛ) - halfar(A) end -to_graphviz(budyko_sellers_halfar_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") + +draw_composition(budyko_sellers_halfar_composition_diagram) ``` We apply a composition by plugging in a Decapode for each component. We also specify the internal name of the variables to be used in combining. @@ -159,6 +171,7 @@ budyko_sellers_halfar_cospan = oapply(budyko_sellers_halfar_composition_diagram, Open(warming, [:A, :Tₛ]), Open(halfar, [:stress_A])]) budyko_sellers_halfar = apex(budyko_sellers_halfar_cospan) + to_graphviz(budyko_sellers_halfar) ``` @@ -176,12 +189,12 @@ to_graphviz(budyko_sellers_halfar) These dynamics will occur on a 1-D manifold (a line). Points near +-π/2 will represent points near the North/ South poles. Points near 0 represent those at the equator. ``` @example DEC -s′ = EmbeddedDeltaSet1D{Bool, Point2D}() -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()) +s = EmbeddedDeltaSet1D{Bool, Point2D}() +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) +sd = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s) +subdivide_duals!(sd, Circumcenter()) ``` ## Define input data @@ -190,15 +203,16 @@ We need to supply initial conditions to our model. We will use synthetic data he ``` @example DEC # This is a primal 0-form, with values at vertices. -cosϕᵖ = map(x -> cos(x[1]), point(s′)) +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 +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 ϕ +α = map(point(s)) do ϕ α₀ + α₂*((1/2)*(3*ϕ[1]^2 - 1)) end A = 210 @@ -207,18 +221,18 @@ f = 0.70 ρ = 1025 cw = 4186 H = 70 -C = map(point(s′)) do ϕ +C = map(point(s)) do ϕ f * ρ * cw * H end D = 0.6 # Isothermal initial conditions: -Tₛ₀ = map(point(s′)) do ϕ +Tₛ₀ = map(point(s)) do ϕ 15 end # Visualize initial condition for temperature. -lines(map(x -> x[1], point(s′)), Tₛ₀) +lines(map(x -> x[1], point(s)), Tₛ₀) ``` ``` @example DEC @@ -227,12 +241,12 @@ n = 3 g = 9.8 # Ice height is a primal 0-form, with values at vertices. -h₀ = map(point(s′)) do (x,_) +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₀) +lines(map(x -> x[1], point(s)), h₀) ``` ``` @example DEC @@ -252,11 +266,23 @@ constants_and_parameters = ( halfar_stress_g = g) ``` -# Symbols to functions +## Symbols to functions -The symbols along edges in our Decapode must be mapped to executable functions. In the Discrete Exterior Calculus, all our operators are defined as relations bewteen points, lines, and triangles on meshes known as simplicial sets. Thus, DEC operators are re-usable across any simplicial set. +The symbols along edges in our Decapode must be mapped to executable functions. In the Discrete Exterior Calculus, all our operators are defined as relations between points, lines, and triangles on meshes known as simplicial sets. Thus, DEC operators are re-usable across any simplicial set. ``` @example DEC +function create_average_matrix(sd) + I = Vector{Int64}() + J = Vector{Int64}() + V = Vector{Float64}() + for e in 1:ne(sd) + append!(J, [sd[e,:∂v0],sd[e,:∂v1]]) + append!(I, [e,e]) + append!(V, [0.5, 0.5]) + end + avg_mat = sparse(I,J,V) +end + function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin :♯ => x -> begin @@ -277,23 +303,11 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs)) end end - :mag => x -> begin - norm.(x) + :mag => x -> norm.(x) + :avg₀₁ => begin + avg_mat = create_average_matrix(sd) + x -> avg_mat * x end - :avg₀₁ => x -> begin - 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) - avg_mat * x - end - :^ => (x,y) -> x .^ y - :* => (x,y) -> x .* y x => error("Unmatched operator $my_symbol") end return (args...) -> op(args...) @@ -306,7 +320,7 @@ From our Decapode, we automatically generate a finite difference method solver t ``` @example DEC sim = eval(gensim(budyko_sellers_halfar, dimension=1)) -fₘ = sim(s, generate) +fₘ = sim(sd, generate) ``` ## Run simulation @@ -316,12 +330,6 @@ We wrap our simulator and initial conditions and solve them with the stability-d ``` @example DEC tₑ = 1e6 -# 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.retcode != :Unstable || error("Solver was not stable") - @info("Solving") prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) soln = solve(prob, Tsit5()) @@ -329,7 +337,7 @@ soln = solve(prob, Tsit5()) @info("Done") ``` -We can save the solution file to examine later. +We can save the [solution file](budyko_sellers_halfar.jld2) to examine later. ``` @example DEC @save "budyko_sellers_halfar.jld2" soln @@ -338,46 +346,47 @@ We can save the solution file to examine later. ## Visualize Quickly examine the final conditions for temperature: + ``` @example DEC -lines(map(x -> x[1], point(s′)), soln(tₑ).Tₛ) +lines(map(x -> x[1], point(s)), soln(tₑ).Tₛ) ``` Quickly examine the final conditions for ice height: + ``` @example DEC -lines(map(x -> x[1], point(s′)), soln(tₑ).halfar_dynamics_h) +lines(map(x -> x[1], point(s)), soln(tₑ).halfar_dynamics_h) ``` -Create animated GIFs of the temperature and ice height dynamics: -``` @example DEC +``` @setup DEC begin # Initial frame frames = 100 -fig = Figure(; size = (800, 800)) +fig = Figure() ax1 = CairoMakie.Axis(fig[1,1]) -xlims!(ax1, extrema(map(x -> x[1], point(s′)))) +xlims!(ax1, extrema(map(x -> x[1], point(s)))) ylims!(ax1, extrema(soln(tₑ).Tₛ)) +ax1.xlabel = "Line plot of temperature from North to South pole, every $(tₑ/frames) time units" 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ₛ) + lines!(fig[1,1], map(x -> x[1], point(s)), soln(t).Tₛ) end end begin # Initial frame frames = 100 -fig = Figure(; size = (800, 800)) +fig = Figure() ax1 = CairoMakie.Axis(fig[1,1]) -xlims!(ax1, extrema(map(x -> x[1], point(s′)))) +xlims!(ax1, extrema(map(x -> x[1], point(s)))) ylims!(ax1, extrema(soln(tₑ).halfar_dynamics_h)) +ax1.xlabel = "Line plot of temperature from North to South pole, every $(tₑ/frames) time units" 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_dynamics_h) + lines!(fig[1,1], map(x -> x[1], point(s)), soln(t).halfar_dynamics_h) end end ``` @@ -385,3 +394,7 @@ end ![BSH_Temperature](budyko_sellers_halfar_T.gif) ![BSH_IceHeight](budyko_sellers_halfar_h.gif) + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/canon.md b/docs/src/canon.md index 3b153432..562d0ddc 100644 --- a/docs/src/canon.md +++ b/docs/src/canon.md @@ -1,24 +1,38 @@ +# Canon + +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "docinfo.jl")) +info = DocInfo.Info() +``` + ## Physics + ```@autodocs Modules = [ Decapodes.Canon.Physics ] Private = false ``` ## Chemistry + ```@autodocs Modules = [ Decapodes.Canon.Chemistry ] Private = false ``` ## Biology + ```@autodocs Modules = [ Decapodes.Canon.Biology ] Private = false ``` ## Environment + ```@autodocs Modules = [ Decapodes.Canon.Environment ] Private = false ``` +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/ch/CahnHilliard_Final.jpg b/docs/src/ch/CahnHilliard_Final.jpg new file mode 100644 index 00000000..08a79d65 Binary files /dev/null and b/docs/src/ch/CahnHilliard_Final.jpg differ diff --git a/docs/src/ch/cahn-hilliard.md b/docs/src/ch/cahn-hilliard.md new file mode 100644 index 00000000..421876a1 --- /dev/null +++ b/docs/src/ch/cahn-hilliard.md @@ -0,0 +1,123 @@ +# The Cahn-Hilliard Equation + +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + +For this example, Decapodes will model the Cahn-Hilliard equation. This equation describes the evolution of a binary fluid as its two phases separate out into distinct domains. Below is a high resolution preview of this model. Notice how the fluid has separated into distinct regions (blue and red) as well as the presence of a transition region. + +!["Cahn Hilliard sample"](CahnHilliard_Final.jpg) + +## Formulating the Equation + +We first load in our dependencies + +```@example DEC +# AlgebraicJulia Dependencies +using Catlab +using CombinatorialSpaces +using Decapodes +using DiagrammaticEquations + +# External Dependencies +using CairoMakie +using ComponentArrays +using GeometryBasics +using LinearAlgebra +using MLStyle +using OrdinaryDiffEq +using Random +Point3D = Point3{Float64}; +nothing #hide +``` + +and then proceed to describe our physics using Decapodes. + +```@example DEC +CahnHilliard = @decapode begin + C::Form0 + (D, γ)::Constant + ∂ₜ(C) == D * Δ(C.^3 - C - γ * Δ(C)) +end + +to_graphviz(CahnHilliard) +``` + +In this equation `C` will represent the concentration of the binary fluid, ranging from `-1` to `1` to differentiate between different phases. We also have a diffusion constant `D` and a constant `γ` whose square root is the length of the transition regions. This formulation of the Cahn-Hilliard equation was drawn from the Wikipedia page on the topic found [here](https://en.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation). + +## Loading the Data + +We now generate the mesh information. We'll run the equation on a triangulated grid. We hide the mesh visualization code for clarity. + +```@example DEC +s = triangulated_grid(100, 100, 0.5, 0.5, Point3D); +sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s); +subdivide_duals!(sd, Circumcenter()); +fig = Figure() # hide +ax = CairoMakie.Axis(fig[1,1], aspect=1) # hide +wf = wireframe!(ax, s; linewidth=1) # hide +save("CahnHilliard_Rect.png", fig) # hide +nothing # hide +``` + +!["CahnHilliardRect"](CahnHilliard_Rect.png) + +The Cahn-Hilliard equation starts with a random concentration holding values between `-1` and `1`. For both `D` and `γ` constants we choose 0.5. + +```@example DEC +Random.seed!(0) + +C = rand(Float64, nv(sd)) * 2 .- 1 +u₀ = ComponentArray(C=C) +constants = (D = 0.5, γ = 0.5); + +fig = Figure() # hide +ax = CairoMakie.Axis(fig[1,1], aspect=1) # hide +msh = CairoMakie.mesh!(ax, s, color=C, colormap=:jet, colorrange=extrema(C)) # hide +Colorbar(fig[1,2], msh) +save("CahnHilliard_initial.png", fig) # hide +nothing # hide +``` + +!["Initial conditions"](CahnHilliard_initial.png) + +We'll now create the simulation code representing the Cahn-Hilliard equation. We pass `nothing` in the second argument to sim since we have no custom functions to pass in. + +```@example DEC +sim = eval(gensim(CahnHilliard)) +fₘ = sim(sd, nothing, DiagonalHodge()); +``` + +## Getting the Solution + +Now that everything is set up and ready, we can solve the equation. We run the simulation for 200 time units to see the long-term evolution of the fluid. Note we only save the solution at intervals of 0.1 time units in order to reduce the memory-footprint of the solve. + +```@example DEC +tₑ = 200 +prob = ODEProblem(fₘ, u₀, (0, tₑ), constants) +soln = solve(prob, Tsit5(), saveat=0.1); +soln.retcode +``` + +And we can see the result as a gif. + +```@setup DEC +function create_gif(solution, file_name) + frames = 200 + fig = Figure() + ax = CairoMakie.Axis(fig[1,1]) + msh = CairoMakie.mesh!(ax, s, color=solution(0).C, colormap=:jet, colorrange=extrema(solution(0).C)) + Colorbar(fig[1,2], msh) + CairoMakie.record(fig, file_name, range(0.0, tₑ; length=frames); framerate = 15) do t + msh.color = solution(t).C + end +end +create_gif(soln, "CahnHilliard_Rect.gif") +``` + +!["CahnHilliardRes"](CahnHilliard_Rect.gif) + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/cism.md b/docs/src/cism/cism.md similarity index 83% rename from docs/src/cism.md rename to docs/src/cism/cism.md index 14f7b2db..f490307d 100644 --- a/docs/src/cism.md +++ b/docs/src/cism/cism.md @@ -1,41 +1,44 @@ # Replicating the Community Ice Sheet Model v2.1 Halfar Dome Benchmark with Decapodes +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + The Decapodes framework takes high-level representations of physics equations and automatically generates solvers. We do so by translating equations from vector calculus notation to the "discrete exterior calculus" (DEC). This process is roughly about recognizing whether physical quantities represent scalar or vector quantities, and recognizing whether differential operators represent gradient, divergence, and so on. -In this benchmark, we will implement the "small slope approximation" of glacial dynamics used by P. Halfar in his 1981 work "On the dynamics of the ice sheets" by taking his original formulation, translating it into the DEC, then providing a mesh and initial conditions. - -The initial conditions used here are exactly those considered by W. H. Lipscomb et al. in "Description And Evaluation of the Community Ice Sheet Model (CISM) v2.1" (2019). +In this benchmark, we will implement the "small slope approximation" of glacial dynamics used by P. Halfar in his 1981 work ["On the dynamics of the ice sheets"](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JC086iC11p11065) by taking his original formulation, translating it into the DEC, then providing a mesh and initial conditions. +The initial conditions used here are exactly those considered by W. H. Lipscomb et al. in ["Description And Evaluation of the Community Ice Sheet Model (CISM) v2.1" (2019)](https://gmd.copernicus.org/articles/12/387/2019/). ```@example DEC # AlgebraicJulia Dependencies using Catlab -using Catlab.Graphics using CombinatorialSpaces using Decapodes using DiagrammaticEquations -using DiagrammaticEquations.Deca # External Dependencies +using BenchmarkTools +using CairoMakie using ComponentArrays -using MLStyle +using GeometryBasics: Point2, Point3 +using JLD2 using LinearAlgebra +using MLStyle using OrdinaryDiffEq -using JLD2 using SparseArrays using Statistics -using CairoMakie -using BenchmarkTools -using GeometryBasics: Point2, Point3 Point2D = Point2{Float64} -Point3D = Point3{Float64}; # hide +Point3D = Point3{Float64} +nothing # hide ``` ## Specifying and Composing Physics -!["Halfar Equation 2"](https://cise.ufl.edu/~luke.morris/halfar_eq2.jpg) +!["Halfar Equation 2"](halfar_eq2.png) We will translate Halfar's equation into the DEC below. Although the equation given by Halfar is dense, this notation does not allow you to see which operators represent divergence, which components represent diffusivity constants, and so on. In the DEC, there is a small pool of operators, ⋆, d, ∧, ♯, and ♭, which combine according to set rules to encode all of these notions. @@ -60,7 +63,7 @@ end to_graphviz(halfar_eq2) ``` -!["Glen's Law"](https://cise.ufl.edu/~luke.morris/glen_law1.jpg) +!["Glen's Law"](glens_law.png) Here, we recognize that Gamma is in fact what glaciologists call "Glen's Flow Law." It states that the strain rate of a sheet of ice can be related to applied stress via a power law. Below, we encode the formulation as it is usually given in the literature, depending explicitly on the gravitational constant, g. @@ -89,7 +92,7 @@ ice_dynamics_composition_diagram = @relation () begin dynamics(Γ,n) stress(Γ,n) end -to_graphviz(ice_dynamics_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(ice_dynamics_composition_diagram) ``` ```@example DEC @@ -120,16 +123,17 @@ We are done encoding our dynamics. Now, we need to provide a mesh, initial data Our mesh library, CombinatorialSpaces, can interpret arbitrary .OBJ files to run our dynamics on. Here, we use a script that generates a triangulated grid of the resolution used in the CISM benchmark. Note though that the "resolution" of a triangulated and non-triangulated grid is difficult to directly compare. ```@example DEC -s′ = triangulated_grid(60_000,100_000,2_000,2_000,Point3D) -s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s′) -subdivide_duals!(s, Barycenter()) -x̄ = mean(p -> p[1], point(s)) -ȳ = mean(p -> p[2], point(s)) +s = triangulated_grid(60_000,100_000,2_000,2_000,Point3D) +sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s) +subdivide_duals!(sd, Barycenter()) +x̄ = mean(p -> p[1], point(sd)) +ȳ = mean(p -> p[2], point(sd)) fig = Figure() ax = CairoMakie.Axis(fig[1,1], aspect=0.6, xticks = [0, 3e4, 6e4]) -wf = wireframe!(ax, s) +wf = wireframe!(ax, sd; linewidth=0.5) save("ice_mesh.png", fig) +nothing # hide ``` !["Wireframe of the Domain"](ice_mesh.png) @@ -150,7 +154,7 @@ g = 9.8101 alpha = 1/9 beta = 1/18 flwa = 1e-16 -A = fill(1e-16, ne(s)) +A = fill(1e-16, ne(sd)) Gamma = 2.0/(n+2) * flwa * (ρ * g)^n t0 = (beta/Gamma) * (7.0/4.0)^3 * (R₀^4 / H^7) @@ -168,11 +172,12 @@ end # Set the initial conditions for ice sheet height: # Ice height is a primal 0-form. i.e. valued at vertices. -h₀ = map(x -> height_at_p(x[1], x[2], 0), point(s′)) +h₀ = map(x -> height_at_p(x[1], x[2], 0), point(s)) fig = Figure() ax = CairoMakie.Axis(fig[1,1], aspect=0.6, xticks = [0, 3e4, 6e4]) -msh = mesh!(ax, s′, color=h₀, colormap=:jet) +msh = mesh!(ax, s, color=h₀, colormap=:jet) save("ice_initial_conditions.png", fig) +nothing # hide ``` !["Initial Conditions"](ice_initial_conditions.png) @@ -197,13 +202,12 @@ We provide here the mapping from symbols to differential operators. As more of t function generate(sd, my_symbol; hodge=GeometricHodge()) # We pre-allocate matrices that encode differential operators. op = @match my_symbol begin - :mag => x -> begin - norm.(x) + :mag => x -> norm.(x) + :♯ => begin + sharp_mat = ♯_mat(sd, AltPPSharp()) + x -> sharp_mat * x end - :^ => (x,y) -> begin - x .^ y - end - _ => default_dec_matrix_generate(sd, my_symbol, hodge) + x => error("Unmatched operator $my_symbol") end return (args...) -> op(args...) end @@ -217,7 +221,7 @@ The `gensim` function takes our high-level representation of the physics equatio ####################### sim = eval(gensim(ice_dynamics2D)) -fₘ = sim(s, generate) +fₘ = sim(sd, generate) ``` Julia is a "Just-In-Time" compiled language. That means that functions are compiled the first time they are called, and later calls to those functions skip this step. To get a feel for just how fast this simulation is, we will run the dynamics twice, once for a very short timespan to trigger pre-compilation, and then again for the actual dynamics. @@ -240,7 +244,7 @@ tₑ = 200 # This next run should be fast. @info("Solving") prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) -soln = solve(prob, Tsit5()) +soln = solve(prob, Tsit5(), saveat=0.1) @show soln.retcode @info("Done") ``` @@ -250,18 +254,19 @@ We can benchmark the compiled simulation with `@benchmarkable`. This macro runs ```@example DEC # Time the simulation -b = @benchmarkable solve(prob, Tsit5()) +b = @benchmarkable solve(prob, Tsit5(), saveat=0.1) c = run(b) ``` +Here we save the solution information to a [file](ice_dynamics2D.jld2). + ```@example DEC -# Save the solution information to a file. @save "ice_dynamics2D.jld2" soln ``` We recall that these dynamics are of the "shallow slope" and "shallow ice" approximations. So, at the edge of our parabolic dome of ice, we expect increased error as the slope increases. On the interior of the dome, we expect the dynamics to match more closely that given by the analytic model. We will see that the CISM results likewise accumulate error in the same neighborhood. -!["Halfar Small Ice Approximation Quote"](https://cise.ufl.edu/~luke.morris/halfar_small_grad_approximation.jpg) +!["Halfar Small Ice Approximation Quote"](halfar_quote.png) ```@example DEC # Plot the final conditions @@ -270,12 +275,13 @@ function plot_final_conditions() ax = CairoMakie.Axis(fig[1,1], title="Modeled thickness (m) at time 200.0", aspect=0.6, xticks = [0, 3e4, 6e4]) - msh = mesh!(ax, s′, color=soln(200.0).dynamics_h, colormap=:jet) + msh = mesh!(ax, s, color=soln(200.0).dynamics_h, colormap=:jet) Colorbar(fig[1,2], msh) fig end fig = plot_final_conditions() save("ice_numeric_solution.png", fig) +nothing # hide ``` !["Numerical Solution"](ice_numeric_solution.png) @@ -283,17 +289,18 @@ save("ice_numeric_solution.png", fig) ```@example DEC # Plot the final conditions according to the analytic solution. function plot_analytic() - hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s′)) + hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s)) fig = Figure() ax = CairoMakie.Axis(fig[1,1], title="Analytic thickness (m) at time 200.0", aspect=0.6, xticks = [0, 3e4, 6e4]) - msh = mesh!(ax, s′, color=hₐ, colormap=:jet) + msh = mesh!(ax, s, color=hₐ, colormap=:jet) Colorbar(fig[1,2], msh) fig end fig = plot_analytic() save("ice_analytic_solution.png", fig) +nothing # hide ``` !["Analytic Solution](ice_analytic_solution.png) @@ -301,19 +308,20 @@ save("ice_analytic_solution.png", fig) ```@example DEC # Plot the error. function plot_error() - hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s′)) + hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s)) h_diff = soln(tₑ).dynamics_h - hₐ extrema(h_diff) fig = Figure() ax = CairoMakie.Axis(fig[1,1], title="Modeled thickness - Analytic thickness at time 200.0", aspect=0.6, xticks = [0, 3e4, 6e4]) - msh = mesh!(ax, s′, color=h_diff, colormap=:jet) + msh = mesh!(ax, s, color=h_diff, colormap=:jet) Colorbar(fig[1,2], msh) fig end fig = plot_error() save("ice_error.png", fig) +nothing # hide ``` !["Numeric Solution - Analytic Solution"](ice_error.png) @@ -322,24 +330,24 @@ We compute below that the maximum absolute error is approximately 89 meters. We ```@example DEC # Compute max absolute error: -hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s′)) +hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s)) h_diff = soln(tₑ).dynamics_h - hₐ maximum(abs.(h_diff)) ``` -We compute root meand square error (RMSE) as well, both over the entire domain, and *excluding where the ice distribution is 0 in the analytic solution.* (Considering the entire domain decreases the RMSE, while not telling you much about the dynamics in the area of interest.) Note that the official CISM benchmark reports 6.43 and 9.06 RMSE for their two solver implementations. +We compute root-mean-square (RMS) error as well, both over the entire domain, and *excluding where the ice distribution is 0 in the analytic solution.* This is done since considering the entire domain decreases the RMS while not telling you much about the area of interest. Note that the official CISM benchmark reports `6.43` and `9.06` RMS for their two solver implementations. ```@example DEC -# Compute RMSE not considering the "outside". -hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s′)) +# Compute RMS not considering the "outside". +hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s)) nonzeros = findall(!=(0), hₐ) h_diff = soln(tₑ).dynamics_h - hₐ rmse = sqrt(sum(map(x -> x*x, h_diff[nonzeros])) / length(nonzeros)) ``` ```@example DEC -# Compute RMSE of the entire domain. -hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s′)) +# Compute RMS of the entire domain. +hₐ = map(x -> height_at_p(x[1], x[2], 200.0), point(s)) h_diff = soln(tₑ).dynamics_h - hₐ rmse = sqrt(sum(map(x -> x*x, h_diff)) / length(h_diff)) ``` @@ -350,7 +358,7 @@ begin frames = 100 fig = Figure() ax = CairoMakie.Axis(fig[1,1], aspect=0.6, xticks = [0, 3e4, 6e4]) - msh = mesh!(ax, s′, color=soln(0).dynamics_h, colormap=:jet, colorrange=extrema(soln(tₑ).dynamics_h)) + msh = mesh!(ax, s, color=soln(0).dynamics_h, colormap=:jet, colorrange=extrema(soln(tₑ).dynamics_h)) Colorbar(fig[1,2], msh) record(fig, "ice_dynamics_cism.gif", range(0.0, tₑ; length=frames); framerate = 30) do t msh.color = soln(t).dynamics_h @@ -364,11 +372,14 @@ For comparison's sake, we paste the results produced by CISM below. We observe t Not that since the DEC is based on triangulated meshes, the "resolution" of the CISM benchmark and the Decapodes implementation cannot be directly compared. An advantage of the DEC is that we do not need to operate on uniform grids. For example, you could construct a mesh that is finer along the dome edge, where you need more resolution, and coarser as you are farther away from the reach of the ice. -![CISM Results](https://cise.ufl.edu/~luke.morris/cism_results.png) +![CISM Results](official_res.png) -We saw in this document how to create performant and accurate simulations in the Decapodes framework, and compared against the CISM library . Although we do not expect to be both more performant and accurate compared to every hand-crafted simulation, Decapodes makes up for this difference in terms of development time, flexibility, and composition. For example, the original implementation of the Decapodes shallow ice model took place over a couple of afternoons during a hackathon. +We saw in this document how to create performant and accurate simulations in the Decapodes framework, and compared against the CISM library . Although we do not expect to be both more performant and accurate compared to every hand-crafted simulation, Decapodes makes up for this difference in terms of development time, flexibility, and composition. For example, the original implementation of the Decapodes shallow ice model took place over a couple of afternoons. Since Decapodes targets high-level representations of physics, it is uniquely suited to incorporating knowledge from subject matter experts to increase simulation accuracy. This process does not require an ice dynamics expert to edit physics equations that have already been weaved into solver code. Further improvements to the Decapodes library are made continuously. We are creating implementations of DEC operators that are constructed and execute faster. And we are in the beginning stages of 3D simulations using the DEC. +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/cism/glens_law.png b/docs/src/cism/glens_law.png new file mode 100644 index 00000000..3c0da180 Binary files /dev/null and b/docs/src/cism/glens_law.png differ diff --git a/docs/src/cism/halfar_eq2.png b/docs/src/cism/halfar_eq2.png new file mode 100644 index 00000000..a0d05a77 Binary files /dev/null and b/docs/src/cism/halfar_eq2.png differ diff --git a/docs/src/cism/halfar_quote.png b/docs/src/cism/halfar_quote.png new file mode 100644 index 00000000..ed189305 Binary files /dev/null and b/docs/src/cism/halfar_quote.png differ diff --git a/docs/src/cism/official_res.png b/docs/src/cism/official_res.png new file mode 100644 index 00000000..c2bd3836 Binary files /dev/null and b/docs/src/cism/official_res.png differ diff --git a/docs/src/equations.md b/docs/src/equations/equations.md similarity index 95% rename from docs/src/equations.md rename to docs/src/equations/equations.md index ca075f26..804f3946 100644 --- a/docs/src/equations.md +++ b/docs/src/equations/equations.md @@ -1,20 +1,21 @@ # Simple Equations +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..","docinfo.jl")) +info = DocInfo.Info() +``` + This tutorial shows how to use Decapodes to represent simple equations. These aren't using any of the Discrete Exterior Calculus or CombinatorialSpaces features of Decapodes. They just are a reference for how to build equations with the `@decapodes` macro and see how they are stored as ACSets. ```@example Harmonic using Catlab -using Catlab.Graphics using CombinatorialSpaces -using CombinatorialSpaces.ExteriorCalculus using DiagrammaticEquations -using DiagrammaticEquations.Deca using Decapodes ``` The harmonic oscillator can be written in Decapodes in at least three different ways. - ```@example Harmonic oscillator = @decapode begin X::Form0 @@ -96,3 +97,7 @@ to_graphviz(oscillator) ``` Often you will have a linear material where you are scaling by a constant, and a nonlinear version of that material where that scaling is replaced by a generic nonlinear function. This is why we allow Decapodes to represent both of these types of equations. + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/Icethickness_Grigoriev_ice_cap_2021.tif b/docs/src/grigoriev/Icethickness_Grigoriev_ice_cap_2021.tif similarity index 100% rename from docs/src/Icethickness_Grigoriev_ice_cap_2021.tif rename to docs/src/grigoriev/Icethickness_Grigoriev_ice_cap_2021.tif diff --git a/docs/src/grigoriev.md b/docs/src/grigoriev/grigoriev.md similarity index 64% rename from docs/src/grigoriev.md rename to docs/src/grigoriev/grigoriev.md index 736f0745..a79db269 100644 --- a/docs/src/grigoriev.md +++ b/docs/src/grigoriev/grigoriev.md @@ -1,44 +1,49 @@ # Halfar's model of glacial flow -Let's model glacial flow using a model of how ice height of a glacial sheet changes over time, from P. Halfar's 1981 paper: "On the dynamics of the ice sheets". +```@setup INFO +include(joinpath(Base.@__DIR__, ".." , "..", "docinfo.jl")) +info = DocInfo.Info() +``` + +Let's model glacial flow using a model of how ice height of a glacial sheet changes over time, from P. Halfar's 1981 paper: ["On the dynamics of the ice sheets"](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JC086iC11p11065). -Let's run the Halfar shallow ice/ shallow slope model on some "real world" data for ice thickness. Van Tricht et al. in their 2023 communication [Measuring and modelling the ice thickness of the Grigoriev ice cap (Kyrgyzstan) and comparison with global dataset](https://tc.copernicus.org/articles/17/4315/2023/tc-17-4315-2023.html) published ice thickness data on an ice cap and stored their data in a TIF. In this document, we will demonstrate how to parse such data and execute a Decapodes model on these initial conditions. +Let's run the Halfar shallow ice / shallow slope model on some "real world" data for ice thickness. Van Tricht et al. in their 2023 communication [Measuring and modelling the ice thickness of the Grigoriev ice cap (Kyrgyzstan) and comparison with global dataset](https://tc.copernicus.org/articles/17/4315/2023/tc-17-4315-2023.html) published ice thickness data on an ice cap and stored their data in a TIF. In this document, we will demonstrate how to parse such data and execute a Decapodes model on these initial conditions. For the parameters to Glen's law, we will use those used in the [Community Ice Sheet Model benchmark](https://cise.ufl.edu/~luke.morris/cism.html). Of course, the parameters of this Kyrgyzstani ice cap likely differ from these by quite some amount, but they are a good place to start. Further, this ice cap does not satisfy the "shallow slope" assumption across the entire domain. ``` @example DEC # AlgebraicJulia Dependencies using Catlab -using Catlab.Graphics using CombinatorialSpaces using DiagrammaticEquations -using DiagrammaticEquations.Deca using Decapodes # External Dependencies +using CairoMakie +using ComponentArrays using FileIO +using GeometryBasics: Point2 using Interpolations -using MLStyle -using ComponentArrays +using JLD2 using LinearAlgebra +using MLStyle using OrdinaryDiffEq -using JLD2 using SparseArrays -using CairoMakie -using GeometryBasics: Point2 Point2D = Point2{Float64} -Point3D = Point3{Float64}; # hide +Point3D = Point3{Float64}; +nothing # hide ``` -# Loading a Scientific Dataset -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`. +## Loading a Scientific Dataset + +The ice thickness data is stored in a TIF that can be downloaded [here](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" ice_thickness_tif = load(file_name) ``` -This data may appear to be a simple binary mask, but that is only because values with no ice are set to `-Inf`. We will account for this we interpolate our data. +This data may visually appear to be a binary mask but that is only because values with no ice are set to `-Inf`. We will account for this when interpolate our data. We use the `Interpolations.jl` library to interpolate this dataset: @@ -50,7 +55,8 @@ const MIN_Y = 243504.5 const MAX_Y = 245599.8 ice_coords = (range(MIN_X, MAX_X, length=size(ice_thickness_tif,1)), range(MIN_Y, MAX_Y, length=size(ice_thickness_tif,2))) -# Note that the tif is set to -floatmax(Float32) where there is no ice. + +# Note that the TIF is set to -floatmax(Float32) where there is no ice. # For our purposes, this is equivalent to 0.0. ice_interp = LinearInterpolation(ice_coords, Float32.(ice_thickness_tif)) ``` @@ -63,31 +69,34 @@ Let's generate a triangulated grid located at the appropriate coordinates: # Specify a resolution: RES_Y = (MAX_Y-MIN_Y)/30.0 RES_X = RES_Y + # Generate the mesh with appropriate dimensions and resolution: -s′ = triangulated_grid( - MAX_X-MIN_X, MAX_Y-MIN_Y, - RES_X, RES_Y, Point3D) +s = triangulated_grid(MAX_X-MIN_X, MAX_Y-MIN_Y, RES_X, RES_Y, Point3D) + # Shift it into place: -s′[:point] = map(x -> x + Point3D(MIN_X, MIN_Y, 0), s′[:point]) -s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s′) -subdivide_duals!(s, Barycenter()) +s[:point] = map(x -> x + Point3D(MIN_X, MIN_Y, 0), s[:point]) +sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s) +subdivide_duals!(sd, Barycenter()) fig = Figure() ax = CairoMakie.Axis(fig[1,1]) wf = wireframe!(ax, s) -display(fig) +save("Grigoriev_IceMesh.png", fig) +nothing # hide ``` -The coordinates of a vertex are stored in `s[:point]`. Let's use our interpolator to assign ice thickness values to each vertex in the mesh: +!["Grigoriev_IceMesh"](Grigoriev_IceMesh.png) + +The coordinates of a vertex are stored in `sd[:point]`. Let's use our interpolator to assign ice thickness values to each vertex in the mesh: ``` @example DEC # These are the values used by the CISM benchmark: n = 3 ρ = 910 g = 9.8101 -A = fill(1e-16, ne(s)) +A = fill(1e-16, ne(sd)) -h₀ = map(s[:point]) do (x,y,_) +h₀ = map(sd[:point]) do (x,y,_) tif_val = ice_interp(x,y) # Accommodate for the -∞'s that encode "no ice". tif_val < 0.0 ? 0.0 : tif_val @@ -95,15 +104,16 @@ end # Store these values to be passed to the solver. u₀ = ComponentArray(h=h₀, stress_A=A) -constants_and_parameters = ( - n = n, - stress_ρ = ρ, - stress_g = g, - stress_A = A) +constants_and_parameters = (n = n, + stress_ρ = ρ, + stress_g = g, + stress_A = A) +nothing # hide ``` -# Defining and Composing Models -For exposition on this Halfar Decapode, see our [Glacial Flow](https://algebraicjulia.github.io/Decapodes.jl/dev/ice_dynamics) docs page. You can skip ahead to the next section. +## Defining and Composing Models + +For exposition on this Halfar Decapode, see our [Glacial Flow](../ice_dynamics/ice_dynamics.md) docs page. Otherwise, you may skip ahead to the next section. ``` @example DEC halfar_eq2 = @decapode begin @@ -112,7 +122,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 @@ -135,67 +145,46 @@ ice_dynamics = apex(ice_dynamics_cospan) to_graphviz(ice_dynamics) ``` -# Define our functions +## 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 + :mag => x -> norm.(x) + :♯ => begin + sharp_mat = ♯_mat(sd, AltPPSharp()) + x -> sharp_mat * x end x => error("Unmatched operator $my_symbol") end - return (args...) -> op(args...) + return op end ``` -# Generate simulation +## Generate simulation ``` @example DEC sim = eval(gensim(ice_dynamics, dimension=2)) -fₘ = sim(s, generate) +fₘ = sim(sd, generate) ``` -# Run +## Run ``` @example DEC -tₑ = 1e1 +tₑ = 10 @info("Solving Grigoriev Ice Cap") prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) soln = solve(prob, Tsit5()) @show soln.retcode @info("Done") + @save "grigoriev.jld2" soln ``` -# Results and Discussion +## Results and Discussion -``` @example DEC +``` @setup DEC # Visualize the initial conditions. function plot_ic() f = Figure() @@ -203,7 +192,7 @@ function plot_ic() title="Grigoriev Ice Cap Initial Thickness [m]", xticks = range(MIN_X, MAX_X; length=5), yticks = range(MIN_Y, MAX_Y; length=5)) - msh = mesh!(ax, s′, color=soln(0.0).h, colormap=:jet) + msh = mesh!(ax, s, color=soln(0.0).h, colormap=:jet) Colorbar(f[1,2], msh) f end @@ -217,7 +206,7 @@ function plot_fc() title="Grigoriev Ice Cap Final Thickness [m]", xticks = range(MIN_X, MAX_X; length=5), yticks = range(MIN_Y, MAX_Y; length=5)) - msh = mesh!(ax, s′, color=soln(tₑ).h, colormap=:jet) + msh = mesh!(ax, s, color=soln(tₑ).h, colormap=:jet) Colorbar(f[1,2], msh) f end @@ -230,7 +219,7 @@ function save_dynamics(save_file_name) h = @lift(soln($time).h) f = Figure() ax = CairoMakie.Axis(f[1,1], title = @lift("Grigoriev Ice Cap Dynamic Thickness [m] at time $($time)")) - gmsh = mesh!(ax, s′, color=h, colormap=:jet, + gmsh = mesh!(ax, s, color=h, colormap=:jet, colorrange=extrema(soln(tₑ).h)) #Colorbar(f[1,2], gmsh, limits=extrema(soln(tₑ).h)) Colorbar(f[1,2], gmsh) @@ -245,5 +234,11 @@ save_dynamics("grigoriev.gif") We observe the usual Halfar model phenomena of ice "melting". Note that since the "shallow slope" approximation does not hold on the boundaries (due to the so-called "ice cliffs" described in the Van Tricht et al. paper), we do not expect the "creep" effect to be physical in this region of the domain. Rather, the Halfar model's predictive power is tuned for the interiors of ice caps and glaciers. Note that we also assume here that the bedrock that the ice rests on is flat. We may in further documents demonstrate how to use topographic data from Digital Elevation Models to inform the elevation of points in the mesh itself. ![Grigoriev_ICs](grigoriev_ic.png) + ![Grigoriev_FCs](grigoriev_fc.png) + ![Grigoriev_Dynamics](grigoriev.gif) + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/halmo.md b/docs/src/halmo/halmo.md similarity index 63% rename from docs/src/halmo.md rename to docs/src/halmo/halmo.md index 509a30e3..889c1b1c 100644 --- a/docs/src/halmo.md +++ b/docs/src/halmo/halmo.md @@ -1,10 +1,15 @@ # Couple Ice and Water Dynamics -Let's use Decapodes to implement the incompressible Navier-Stokes as given by [Mohamed et al.](https://arxiv.org/abs/1508.01166). We will run these dynamics [on the sphere](https://youtu.be/k0hFhAvhHvs?si=Wi9-OgBbAODtxMtb). We will couple this model with Halfar glacier dynamics [on the sphere](https://algebraicjulia.github.io/Decapodes.jl/dev/ice_dynamics/#2-Manifold-in-3D). For the initial conditions of the Halfar ice thickness, we will use scientific dataset for Greenland, much like the scientific dataset used for the [Grigoriev ice cap](https://algebraicjulia.github.io/Decapodes.jl/dev/grigoriev/) simulation. +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + +Let's use Decapodes to implement the incompressible Navier-Stokes as given by [Mohamed et al.](https://arxiv.org/abs/1508.01166). We will run these dynamics [on the sphere](https://youtu.be/k0hFhAvhHvs?si=Wi9-OgBbAODtxMtb). We will couple this model with Halfar glacier dynamics [on the sphere](https://algebraicjulia.github.io/Decapodes.jl/dev/ice_dynamics/#2-Manifold-in-3D). For the initial conditions of the Halfar ice thickness, we will use an idealized polar ice cap. Note that the time scale at which ice creeps is much larger than the time scale at which the water in the ocean would flow. So we can either choose to model a very slow moving fluid around the ice (like a storm on a gas giant), or we can choose to model on a shorter timescale, on which the ice does not move very much. -```@example DEC +```@example DEC_halmo # AlgebraicJulia Dependencies using Catlab using CombinatorialSpaces @@ -20,14 +25,16 @@ using LinearAlgebra using MLStyle using OrdinaryDiffEq Point3D = Point3{Float64}; +nothing # hide ``` -## Specify our models. +## Specify our models Our first component is the Mohamed et al. formulation of the incompressible Navier-Stokes equations. We will call the flow here "w". This will be the flow after collisions with glaciers are considered. This is [Equation 10](https://arxiv.org/abs/1508.01166) for N=2. -```@example DEC + +```@example DEC_halmo eq10forN2 = @decapode begin (𝐮,w)::DualForm1 (P, 𝑝ᵈ)::DualForm0 @@ -42,7 +49,7 @@ to_graphviz(eq10forN2) Halfar's equation and Glen's law are composed like so: -```@example DEC +```@example DEC_halmo halfar_eq2 = @decapode begin h::Form0 Γ::Form1 @@ -73,24 +80,23 @@ to_graphviz(ice_dynamics, verbose=false) We now have our dynamics that govern glaciers, and our dynamics that govern water. We need to specify the physics of what happens when glaciers and water interact. There are many options, and the choice you make depends on the time-scale and resolution of the dynamics that you are interested in. -An interaction between glacier and water dynamics can look like this: +An interaction between glacier and water dynamics can look like the following, where `flow_after` is the flow of water after interaction with ice is considered. -```@example DEC +```@example DEC_halmo ice_water_composition_diagram = @relation () begin glacier_dynamics(ice_thickness) water_dynamics(flow, flow_after) interaction(ice_thickness, flow, flow_after) end -to_graphviz(ice_water_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(ice_water_composition_diagram) ``` -, where `flow_after` is the flow of water after interaction with ice is considered. We will use the language of Decapodes to encode the dynamics that ice blocks water from flowing. -We can detect the ice with a sigmoid function. Where there is ice, we want the flow to be 0, and where there is no ice, we will not impede the flow. We won't consider any further special boundary conditions between ice and water here. Since h is a scalar-like quantity, and flow is a vector-like quantity, we can relate them using the wedge product operator from the exterior calculus. We can state these dynamics using the language of Decapodes like so: +We can detect the ice with a sigmoid function. Where there is ice, we want the flow to be 0, and where there is no ice, we will not impede the flow. We won't consider any further special boundary conditions between ice and water here. Since `h` is a scalar-like quantity, and flow is a vector-like quantity, we can relate them using the wedge product operator from the exterior calculus. We can state these dynamics using the language of Decapodes like so: -```@example DEC +```@example DEC_halmo blocking = @decapode begin h::Form0 (𝐮,w)::DualForm1 @@ -99,10 +105,12 @@ blocking = @decapode begin end to_graphviz(blocking) ``` -, where σ is a sigmoid function that is 0 when d(h) is 0, and goes to 1 otherwise. We see that w is indeed defined as 𝐮, after interacting with the ice boundary is considered. + +Here, `σ` is a sigmoid function that is 0 when d(h) is 0, and goes to 1 otherwise. We see that `w` is indeed defined as `𝐮`, after interacting with the ice boundary is considered. We can apply our composition diagram to generate our physics: -```@example DEC + +```@example DEC_halmo ice_water = apex(oapply(ice_water_composition_diagram, [Open(ice_dynamics, [:dynamics_h]), Open(eq10forN2, [:𝐮, :w]), @@ -113,72 +121,67 @@ to_graphviz(ice_dynamics, verbose=false) We can now generate our simulation: -```@example DEC -open("ice_water.jl", "w") do f - write(f, string(gensim(expand_operators(ice_water)))) -end -sim = include("ice_water.jl") +```@example DEC_halmo +sim = eval(gensim(ice_water)) ``` ## Meshes and Initial Conditions Since we want to demonstrate these physics on the Earth, we will use one of our icosphere discretizations with the appropriate radius. -``` @example DEC +``` @example DEC_halmo rₑ = 6378e3 # [km] -s′ = loadmesh(Icosphere(5, rₑ)) -s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s′) -subdivide_duals!(s, Barycenter()) -wireframe(s) +s = loadmesh(Icosphere(5, rₑ)) +sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s) +subdivide_duals!(sd, Barycenter()) +wireframe(sd) ``` Let's demonstrate how to add operators by providing the definition of a sigmoid function: -```@example DEC +```@example DEC_halmo sigmoid(x) = (2 ./ (1 .+ exp.(-x*1e2)) .- 1) -function generate(s, my_symbol; hodge=GeometricHodge()) +function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin # This is a new function. :σ => sigmoid :mag => x -> norm.(x) - :^ => (x,y) -> x .^ y # Remaining operations (such as our differential operators) are built-in. - _ => default_dec_matrix_generate(s, my_symbol, hodge) + _ => default_dec_matrix_generate(sd, my_symbol, hodge) end - return (args...) -> op(args...) + return op end; ``` Let's combine our mesh with our physics to instantiate our simulation: -```@example DEC -fₘ = sim(s, generate); +```@example DEC_halmo +fₘ = sim(sd, generate); ``` We can now supply initial conditions: -```@example DEC -#TODO: Grab ice data from interpolator. - -ice_thickness = map(s[:point]) do (_,_,z) +```@example DEC_halmo +ice_thickness = map(sd[:point]) do (_,_,z) z < 0.8*rₑ ? 0 : 1 end -flow = dec_dual_derivative(0,s) * - map(s[s[:tri_center], :dual_point]) do (_,_,z) +flow = dec_dual_derivative(0,sd) * + map(sd[sd[:tri_center], :dual_point]) do (_,_,z) (rₑ-abs(z))/rₑ end + # There is no water "under" the ice: -flow = dec_wedge_product_pd(Tuple{0,1},s)(1 .- sigmoid(ice_thickness), flow) +flow = dec_wedge_product_pd(Tuple{0,1},sd)(1 .- sigmoid(ice_thickness), flow) u₀ = ComponentArray( ice_thickness = ice_thickness, flow = flow, - water_dynamics_P = zeros(ntriangles(s))) + water_dynamics_P = zeros(ntriangles(sd))) constants_and_parameters = ( glacier_dynamics_n = 3, - glacier_dynamics_stress_A = fill(1e-16, ne(s)), + glacier_dynamics_stress_A = fill(1e-16, ne(sd)), glacier_dynamics_stress_ρ = 910, glacier_dynamics_stress_g = 9.8101, water_dynamics_μ = 0.01); @@ -188,8 +191,8 @@ constants_and_parameters = ( We specified our physics, our mesh, and our initial conditions. We have everything we need to execute the simulation. -```@example DEC -tₑ = 1e2 +```@example DEC_halmo +tₑ = 100 # Julia will pre-compile the generated simulation the first time it is run. @info("Precompiling Solver") @@ -206,50 +209,14 @@ soln = solve(prob, Vern7()) ## Results -In the DEC, vorticity is encoded with `d⋆`, and speed can be encoded with `norm ♯`. We can use our operators from CombinatorialSpaces.jl to create our GIFs. - -```@example DEC -ihs0 = dec_inv_hodge_star(Val{0}, s, GeometricHodge()) -dd1 = dec_dual_derivative(1, s) -♯_m = ♯_mat(s, LLSDDSharp()) -using LinearAlgebra: norm -function vorticity(α) - ihs0*dd1*α -end -function speed(α) - norm.(♯_m * α) -end +Let's look at the dynamics of the ice: -function save_vorticity(is_2d=false) - frames = 200 - time = Observable(0.0) - fig = Figure(title = @lift("Vorticity at $($time)")) - ax = is_2d ? - CairoMakie.Axis(fig[1,1]) : - LScene(fig[1,1], scenekw=(lights=[],)) - msh = CairoMakie.mesh!(ax, s′, - color=@lift(vorticity(soln($time).flow)), - colorrange=extrema(vorticity(soln(tₑ).flow)).*.9, - colormap=:jet) - - Colorbar(fig[1,2], msh) - record(fig, "vorticity_ice_water.gif", range(0.0, tₑ; length=frames); framerate = 20) do t - time[] = t - end -end -save_vorticity(false) -``` - -![Vorticity](vorticity_ice_water.gif) - -Let's look at the dynamics of the ice as well: - -``` @example DEC +``` @example DEC_halmo begin frames = 200 fig = Figure() ax = LScene(fig[1,1], scenekw=(lights=[],)) - msh = CairoMakie.mesh!(ax, s′, color=soln(0).ice_thickness, colormap=:jet, colorrange=extrema(soln(0).ice_thickness)) + msh = CairoMakie.mesh!(ax, s, color=soln(0).ice_thickness, colormap=:jet, colorrange=extrema(soln(0).ice_thickness)) Colorbar(fig[1,2], msh) record(fig, "halmo_ice.gif", range(0.0, tₑ; length=frames); framerate = 20) do t @@ -260,3 +227,6 @@ end ![HalfarMohamedIce](halmo_ice.gif) +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/ice_dynamics.md b/docs/src/ice_dynamics/ice_dynamics.md similarity index 74% rename from docs/src/ice_dynamics.md rename to docs/src/ice_dynamics/ice_dynamics.md index 14b002c9..b880cd5d 100644 --- a/docs/src/ice_dynamics.md +++ b/docs/src/ice_dynamics/ice_dynamics.md @@ -1,48 +1,53 @@ # Halfar's model of glacial flow -Let's model glacial flow using a model of how ice height of a glacial sheet changes over time, from P. Halfar's 1981 paper: "On the dynamics of the ice sheets". +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + +Let's model glacial flow using a model of how ice height of a glacial sheet changes over time, from P. Halfar's 1981 paper: ["On the dynamics of the ice sheets"](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JC086iC11p11065). ``` @example DEC # AlgebraicJulia Dependencies using Catlab -using Catlab.Graphics using CombinatorialSpaces using DiagrammaticEquations -using DiagrammaticEquations.Deca using Decapodes # External Dependencies -using MLStyle +using CairoMakie using ComponentArrays +using GeometryBasics: Point2, Point3 +using JLD2 using LinearAlgebra +using MLStyle using OrdinaryDiffEq -using JLD2 using SparseArrays using Statistics -using CairoMakie -using GeometryBasics: Point2, Point3 Point2D = Point2{Float64}; Point3D = Point3{Float64}; ``` -# Defining the models +## Defining the models The first step is to find a suitable equation for our model, and translate it into the Discrete Exterior Calculus. The Exterior Calculus is a generalization of vector calculus, so for low-dimensional spaces, this translation is straightforward. For example, divergence is typically written as (⋆, d, ⋆). Scalar fields are typically interpreted as "0Forms", i.e. values assigned to vertices of a mesh. We use the `@decapode` macro to interpret the equations. Here, we have equation 2 from Halfar: + ```math \frac{\partial h}{\partial t} = \frac{2}{n + 2} (\frac{\rho g}{A})^n \frac{\partial}{\partial x}(\frac{\partial h}{\partial x} |\frac{\partial h}{\partial x}| ^{n-1} h^{n+2}). ``` + We'll change the term out front to Γ so we can demonstrate composition in a moment. In the exterior calculus, we could write the above equations like so: + ```math \partial_t(h) = \circ(\star, d, \star)(\Gamma\quad d(h)\quad \text{avg}_{01}|d(h)^\sharp|^{n-1} \quad \text{avg}_{01}(h^{n+2})). ``` `avg` here is an operator that performs the midpoint rule, setting the value at an edge to be the average of the values at its two vertices. - ``` @example DEC halfar_eq2 = @decapode begin h::Form0 @@ -56,7 +61,8 @@ end to_graphviz(halfar_eq2) ``` -And here, a formulation of Glen's law from J.W. Glen's 1958 "The flow law of ice". +And here, a formulation of Glen's law from J.W. Glen's 1958 ["The flow law of ice"](http://go.owu.edu/~chjackso/Climate/papers/Glen_1958_The%20flow%20law%20of%20ice.pdf). + ``` @example DEC glens_law = @decapode begin #Γ::Form0 @@ -69,9 +75,9 @@ end to_graphviz(glens_law) ``` -# Composing models +## Composing models -We can use operadic composition to specify how our models come together. In this example, we have two Decapodes, and two quantities that are shared between them. +We can use [operadic composition](https://arxiv.org/abs/2105.12282) to specify how our models come together. In this example, we have two Decapodes, and two quantities that are shared between them. ``` @example DEC ice_dynamics_composition_diagram = @relation () begin @@ -79,7 +85,7 @@ ice_dynamics_composition_diagram = @relation () begin stress(Γ,n) end -to_graphviz(ice_dynamics_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(ice_dynamics_composition_diagram) ``` To a apply a composition, we specify which Decapodes to plug into those boxes, and what each calls the corresponding shared variables internally. @@ -93,7 +99,7 @@ ice_dynamics = apex(ice_dynamics_cospan) to_graphviz(ice_dynamics) ``` -# Provide a semantics +## Provide a semantics To interpret our composed Decapode, we need to specify what Discrete Exterior Calculus to interpret our quantities in. Let's choose the 1D Discrete Exterior Calculus: @@ -105,22 +111,25 @@ resolve_overloads!(ice_dynamics1D, op1_res_rules_1D, op2_res_rules_1D) to_graphviz(ice_dynamics1D) ``` -# Define a mesh +## Define a mesh We'll need a mesh to simulate on. Since this is a 1D mesh, we can go ahead and make one right now: ``` @example DEC -# This is a 1D mesh, consisting of edges and vertices. -s′ = EmbeddedDeltaSet1D{Bool, Point2D}() +# This is an empty 1D mesh. +s = EmbeddedDeltaSet1D{Bool, Point2D}() + # 20 vertices along a line, connected by edges. -add_vertices!(s′, 20, point=Point2D.(range(0, 10_000, length=20), 0)) -add_edges!(s′, 1:nv(s′)-1, 2:nv(s′)) -orient!(s′) -s = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s′) -subdivide_duals!(s, Circumcenter()) +add_vertices!(s, 20, point=Point2D.(range(0, 10_000, length=20), 0)) +add_edges!(s, 1:nv(s)-1, 2:nv(s)) +orient!(s) + +# The dual 1D mesh +sd = EmbeddedDeltaDualComplex1D{Bool, Float64, Point2D}(s) +subdivide_duals!(sd, Circumcenter()) ``` -# Define input data +## Define input data We need initial conditions to use for our simulation. @@ -132,12 +141,12 @@ A = 1e-16 # Ice height is a primal 0-form, with values at vertices. # We choose a distribution that obeys the shallow height and shallow slope conditions. -h₀ = map(point(s′)) do (x,_) +h₀ = map(point(s)) do (x,_) ((7072-((x-5000)^2))/9e3+2777)/2777e-1 end # Visualize initial conditions for ice sheet height. -lines(map(x -> x[1], point(s′)), h₀, linewidth=5) +lines(map(x -> x[1], point(s)), h₀, linewidth=5) ``` We need to tell our Decapode which data maps to which symbols. We can wrap up our data like so: @@ -152,12 +161,23 @@ constants_and_parameters = ( stress_A = A) ``` -# Define functions +## Define functions In order to solve our equations, we will need numerical linear operators that give meaning to our symbolic operators. In the DEC, there are a handful of operators that one uses to construct all the usual vector calculus operations, namely: ♯, ♭, ∧, d, ⋆. The CombinatorialSpaces.jl library specifies many of these for us. - ``` @example DEC +function create_average_matrix(sd) + I = Vector{Int64}() + J = Vector{Int64}() + V = Vector{Float64}() + for e in 1:ne(sd) + append!(J, [sd[e,:∂v0],sd[e,:∂v1]]) + append!(I, [e,e]) + append!(V, [0.5, 0.5]) + end + avg_mat = sparse(I,J,V) +end + function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin :♯ => x -> begin @@ -178,26 +198,10 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) sum([nv*norm(nv)*x[e] for (e,nv) in zip(es,nvs)]) / sum(norm.(nvs)) end end - :mag => x -> begin - norm.(x) - end - :avg₀₁ => x -> begin - 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) - avg_mat * x - end - :^ => (x,y) -> x .^ y - :* => (x,y) -> x .* y - :show => x -> begin - @show x - x + :mag => x -> norm.(x) + :avg₀₁ => begin + avg_mat = create_average_matrix(sd) + x -> avg_mat * x end x => error("Unmatched operator $my_symbol") end @@ -205,28 +209,26 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) end ``` -# Generate the simulation +## Generate the simulation Now, we have everything we need to generate our simulation: - ``` @example DEC sim = eval(gensim(ice_dynamics1D, dimension=1)) -fₘ = sim(s, generate) +fₘ = sim(sd, generate) ``` -# Pre-compile and run +## Pre-compile and run The first time that you run a function, Julia will pre-compile it, so that later runs will be fast. We'll solve our simulation for a short time span, to trigger this pre-compilation, and then run it. - ``` @example DEC @info("Precompiling Solver") prob = ODEProblem(fₘ, u₀, (0, 1e-8), constants_and_parameters) soln = solve(prob, Tsit5()) soln.retcode != :Unstable || error("Solver was not stable") -tₑ = 8e3 +tₑ = 8_000 # This next run should be fast. @info("Solving") @@ -236,18 +238,18 @@ soln = solve(prob, Tsit5()) @info("Done") ``` -We can save our solution file in case we want to examine its contents when this Julia session ends. +We can save our [solution file](ice_dynamics1D.jld2) in case we want to examine its contents when this Julia session ends. ``` @example DEC @save "ice_dynamics1D.jld2" soln ``` -# Visualize +## Visualize Let's examine the final conditions: ``` @example DEC -fig,ax,ob = lines(map(x -> x[1], point(s′)), soln(tₑ).dynamics_h, linewidth=5) +fig,ax,ob = lines(map(x -> x[1], point(s)), soln(tₑ).dynamics_h, linewidth=5) ylims!(ax, extrema(h₀)) fig ``` @@ -260,17 +262,17 @@ Let's create a GIF to examine an animation of these dynamics: # Create a gif begin frames = 100 - fig, ax, ob = lines(map(x -> x[1], point(s′)), soln(0).dynamics_h) + fig, ax, ob = lines(map(x -> x[1], point(s)), soln(0).dynamics_h) ylims!(ax, extrema(h₀)) record(fig, "ice_dynamics1D.gif", range(0.0, tₑ; length=frames); framerate = 15) do t - lines!(map(x -> x[1], point(s′)), soln(t).dynamics_h) + lines!(map(x -> x[1], point(s)), soln(t).dynamics_h) end end ``` ![IceDynamics1D](ice_dynamics1D.gif) -# 2D Re-interpretation +## 2D Re-interpretation The first, one-dimensional, semantics that we provided to our Decapode restricted the kinds of glacial sheets that we could model. (i.e. We could only look at glacial sheets which were constant along y). We can give our Decapode an alternate semantics, as some physics on a 2-dimensional manifold. @@ -284,7 +286,7 @@ resolve_overloads!(ice_dynamics2D) to_graphviz(ice_dynamics2D) ``` -# Store as JSON +## Store as JSON We quickly demonstrate how to serialize a Decapode to JSON and read it back in: @@ -292,6 +294,8 @@ We quickly demonstrate how to serialize a Decapode to JSON and read it back in: write_json_acset(ice_dynamics2D, "ice_dynamics2D.json") ``` +You can view the JSON file [here](ice_dynamics2D.json). + ``` @example DEC # When reading back in, we specify that all attributes are "Symbol"s. ice_dynamics2 = read_json_acset(SummationDecapode{Symbol,Symbol,Symbol}, "ice_dynamics2D.json") @@ -301,20 +305,20 @@ ice_dynamics3 = read_json_acset(SummationDecapode{String,String,String}, "ice_dy to_graphviz(ice_dynamics3) ``` -# Define our mesh +## Define our mesh ``` @example DEC -s′ = triangulated_grid(10_000,10_000,800,800,Point3D) -s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s′) -subdivide_duals!(s, Barycenter()) +s = triangulated_grid(10_000,10_000,800,800,Point3D) +sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s) +subdivide_duals!(sd, Barycenter()) fig = Figure() ax = CairoMakie.Axis(fig[1,1]) -wf = wireframe!(ax, s′) +wf = wireframe!(ax, s) fig ``` -# Define our input data +## Define our input data ``` @example DEC n = 3 @@ -323,15 +327,15 @@ g = 9.8 A = 1e-16 # Ice height is a primal 0-form, with values at vertices. -h₀ = map(point(s′)) do (x,y) +h₀ = map(point(s)) do (x,y) (7072-((x-5000)^2 + (y-5000)^2)^(1/2))/9e3+10 end # Visualize initial condition for ice sheet height. -mesh(s′, color=h₀, colormap=:jet) +mesh(s, color=h₀, colormap=:jet) fig = Figure() ax = CairoMakie.Axis(fig[1,1]) -msh = mesh!(ax, s′, color=h₀, colormap=:jet) +msh = mesh!(ax, s, color=h₀, colormap=:jet) Colorbar(fig[1,2], msh) fig ``` @@ -346,35 +350,19 @@ constants_and_parameters = ( stress_A = A) ``` -# Define our functions +## Define our functions ``` @example DEC function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin - :♯ => x -> begin - ♯(sd, EForm(x)) - end - :mag => x -> begin - norm.(x) + :♯ => begin + sharp_mat = ♯_mat(sd, AltPPSharp()) + x -> sharp_mat * x end - :avg₀₁ => x -> begin - 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) - avg_mat * x - end - :^ => (x,y) -> x .^ y - :* => (x,y) -> x .* y - :show => x -> begin - @show x - @show length(x) - x + :mag => x -> norm.(x) + :avg₀₁ => begin + avg_mat = create_average_matrix(sd) + x -> avg_mat * x end x => error("Unmatched operator $my_symbol") end @@ -382,14 +370,14 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) end ``` -# Generate simulation +## Generate simulation ``` @example DEC sim = eval(gensim(ice_dynamics2D, dimension=2)) -fₘ = sim(s, generate) +fₘ = sim(sd, generate) ``` -# Pre-compile and run +## Pre-compile and run 2D ``` @example DEC @info("Precompiling Solver") @@ -413,13 +401,13 @@ soln = solve(prob, Tsit5()) @save "ice_dynamics2D.jld2" soln ``` -# Visualize +## Visualize 2D ``` @example DEC # Final conditions: fig = Figure() ax = CairoMakie.Axis(fig[1,1]) -msh = mesh!(ax, s′, color=soln(tₑ).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) +msh = mesh!(ax, s, color=soln(tₑ).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) Colorbar(fig[1,2], msh) fig ``` @@ -429,7 +417,7 @@ begin frames = 100 fig = Figure() ax = CairoMakie.Axis(fig[1,1]) - msh = CairoMakie.mesh!(ax, s′, color=soln(0).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) + msh = CairoMakie.mesh!(ax, s, color=soln(0).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) Colorbar(fig[1,2], msh) record(fig, "ice_dynamics2D.gif", range(0.0, tₑ; length=frames); framerate = 15) do t msh.color = soln(t).dynamics_h @@ -439,16 +427,15 @@ end ![IceDynamics2D](ice_dynamics2D.gif) -# 2-Manifold in 3D +## 2-Manifold in 3D We note that just because our physics is happening on a 2-manifold, (a surface), this doesn't restrict us to the 2D plane. In fact, we can "embed" our 2-manifold in 3D space to simulate a glacial sheets spread across the globe. ``` @example DEC -s′ = loadmesh(Icosphere(3, 10_000)) -s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s′) -subdivide_duals!(s, Barycenter()) -wireframe(s) - +s = loadmesh(Icosphere(3, 10_000)) +sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s) +subdivide_duals!(sd, Barycenter()) +wireframe(sd) ``` ``` @example DEC @@ -458,7 +445,7 @@ g = 9.8 A = 1e-16 # Ice height is a primal 0-form, with values at vertices. -h₀ = map(point(s′)) do (x,y,z) +h₀ = map(point(s)) do (x,y,z) (z*z)/(10_000*10_000) end @@ -466,7 +453,7 @@ end # There is lots of ice at the poles, and no ice at the equator. fig = Figure() ax = LScene(fig[1,1], scenekw=(lights=[],)) -msh = CairoMakie.mesh!(ax, s′, color=h₀, colormap=:jet) +msh = CairoMakie.mesh!(ax, s, color=h₀, colormap=:jet) Colorbar(fig[1,2], msh) fig ``` @@ -483,12 +470,12 @@ constants_and_parameters = ( ``` @example DEC sim = eval(gensim(ice_dynamics2D, dimension=2)) -fₘ = sim(s, generate) +fₘ = sim(sd, generate) ``` -``` @example DEC -# For brevity's sake, we'll skip the pre-compilation cell. +For brevity's sake, we'll skip the pre-compilation cell. +``` @example DEC tₑ = 5e25 @info("Solving") @@ -504,7 +491,7 @@ extrema(soln(0).dynamics_h), extrema(soln(tₑ).dynamics_h) ``` @example DEC fig = Figure() ax = LScene(fig[1,1], scenekw=(lights=[],)) -msh = CairoMakie.mesh!(ax, s′, color=soln(tₑ).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) +msh = CairoMakie.mesh!(ax, s, color=soln(tₑ).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) Colorbar(fig[1,2], msh) fig ``` @@ -514,7 +501,7 @@ begin frames = 200 fig = Figure() ax = LScene(fig[1,1], scenekw=(lights=[],)) - msh = CairoMakie.mesh!(ax, s′, color=soln(0).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) + msh = CairoMakie.mesh!(ax, s, color=soln(0).dynamics_h, colormap=:jet, colorrange=extrema(soln(0).dynamics_h)) Colorbar(fig[1,2], msh) # These particular initial conditions diffuse quite quickly, so let's just look at @@ -526,3 +513,7 @@ end ``` ![IceDynamics2DSphere](ice_dynamics2D_sphere.gif) + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/index.md b/docs/src/index.md index 46418ac3..77545be5 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -3,18 +3,18 @@ Decapodes.jl is a framework for developing, composing, and simulating physical systems. Decapodes.jl is the synthesis of Applied Category Theory (ACT) techniques for formalizing and composing physics equations, and Discrete Exterior Calculus (DEC) techniques for formalizing differential operators. -[CombinatorialSpaces.jl](https://algebraicjulia.github.io/CombinatorialSpaces.jl/dev/) hosts tools for discretizing space and defining DEC operators on simplicial complexes, and [DiagrammaticEquations.jl](https://github.com/AlgebraicJulia/DiagrammaticEquations.jl) hosts tooling for representing the equations as formal ACT diagrams. This repository combines these two packages, compiling diagrams down to simulatable code. +[CombinatorialSpaces.jl](https://algebraicjulia.github.io/CombinatorialSpaces.jl/dev/) hosts tools for discretizing space and defining DEC operators on simplicial complexes, and [DiagrammaticEquations.jl](https://github.com/AlgebraicJulia/DiagrammaticEquations.jl) hosts tooling for representing the equations as formal ACT diagrams. This repository combines these two packages, compiling diagrams down to runnable simulation code. By combining the power of ACT and the DEC, we seek to improve the scientific computing workflow. Decapodes simulations are [hierarchically composable](https://algebraicjulia.github.io/Decapodes.jl/dev/budyko_sellers_halfar/), generalize over [any type of manifold](https://algebraicjulia.github.io/Decapodes.jl/dev/ice_dynamics/), and are [performant and accurate](https://www.cise.ufl.edu/~luke.morris/cism.html) with a declarative domain specific language (DSL) that is [human-readable](https://algebraicjulia.github.io/Decapodes.jl/dev/klausmeier/#Model-Representation). ![Grigoriev Ice Cap Dynamics](https://algebraicjulia.github.io/Decapodes.jl/dev/grigoriev.gif) -# Getting started +## Getting started -Walkthroughs creating, composing, and solving Decapodes are available in the side-bar of this documentation page. Further example scripts are avaiable in the `examples` folder of the Decapodes.jl GitHub repo, and will be added to this documentation site soon. +Walkthroughs creating, composing, and solving Decapodes are available in the side-bar of this documentation page. Further example scripts are available in the `examples` folder of the Decapodes.jl GitHub repo, and will be added to this documentation site soon. -# NOTE -This library is currently under active development, and so is not yet at a -point where a constant API/behavior can be assumed. That being said, if this -project looks interesting/relevant please contact us and -[let us know](https://www.algebraicjulia.org/#contributing)! +!!! note "Under Active Development" + This library is currently under active development, and so is not yet at a + point where a constant API/behavior can be assumed. That being said, if this + project looks interesting/relevant please contact us and + [let us know](https://www.algebraicjulia.org/#contributing)! diff --git a/docs/src/klausmeier.md b/docs/src/klausmeier/klausmeier.md similarity index 91% rename from docs/src/klausmeier.md rename to docs/src/klausmeier/klausmeier.md index 1f3f9e81..d5723e0c 100644 --- a/docs/src/klausmeier.md +++ b/docs/src/klausmeier/klausmeier.md @@ -1,10 +1,15 @@ # Klausmeier +```@setup INFO +include(joinpath(Base.@__DIR__, ".." , "..", "docinfo.jl")) +info = DocInfo.Info() +``` + +![Somaliland Vegetation](somalia.png) + ```@raw html
- Somaliland Vegetation -
One of the first aerial photographs of British Somaliland (now Somaliland) investigated by W.A. Macfadyen in his 1950 "Vegetation Patterns in the Semi-Desert Plains of British Somaliland"\[1\]. From this point of view, Macfadyen's "vegetation arcs" are plainly visible.
+
One of the first aerial photographs of British Somaliland (now Somaliland) investigated by W.A. Macfadyen in his 1950 "Vegetation Patterns in the Semi-Desert Plains of British Somaliland" [1]. From this point of view, Macfadyen's "vegetation arcs" are plainly visible.
``` @@ -12,7 +17,7 @@ From aerial photographs in the late 1940s, British ecologist W.A. Macfadyen discovered that vegetation in semi-arid environments often grows in striping patterns, but was unaware of the exact mechanism that causes them. What is especially striking about these "vegetation arcs" is that these stripes appear to climb uphill, with denser plant growth at the leading edge of these traveling waves. Much like how the Mandelbrot set and other interesting fractal patterns can arise from simple sets of rules, these vegetation dynamics can be explained by simple sets of partial differential equations. -The Klausmeier model, given by Christopher Klausmeier in his 1999 paper "Regular and Irregular Patterns in Semiarid Vegetation"\[2\], models such dynamics. Although Macfadyen had discovered these vegetation patterns 50s years prior\[1,3\], defining these dynamics through accessible and physically-meaningful PDEs proved a catalyst for further research. At the time of writing, Klausmeier's paper has been cited 594 times. +The Klausmeier model, given by Christopher Klausmeier in his 1999 paper *Regular and Irregular Patterns in Semiarid Vegetation*\[2\], models such dynamics. Although Macfadyen had discovered these vegetation patterns 50s years prior\[1,3\], defining these dynamics through accessible and physically-meaningful PDEs proved a catalyst for further research. At the time of writing, Klausmeier's paper has been cited 594 times. In this document, we will use Decapodes to formally represent these equations. Moreover, we will demonstrate how one can automatically generate simulation that reproduces the dynamics given by a scientist, simply by reading in the equations given in their original publication. @@ -21,23 +26,24 @@ The lofty goal of this document, and of Decapodes itself, is that through both e ![Klausmeier GIF](klausmeier.gif) ## using Decapodes + ```@example DEC # Load Dependencies -using DiagrammaticEquations -using DiagrammaticEquations.Deca -using Decapodes +using CairoMakie using Catlab using CombinatorialSpaces - +using ComponentArrays +using Decapodes +using DiagrammaticEquations +using DiagrammaticEquations.Deca using Distributions -using CairoMakie +using GeometryBasics: Point2 using JLD2 using LinearAlgebra using MLStyle -using ComponentArrays using OrdinaryDiffEq -using GeometryBasics: Point2 Point2D = Point2{Float64} +nothing # hide ``` ## Model Representation @@ -53,7 +59,7 @@ Hydrodynamics = @decapode begin dX::Form1 (a,ν)::Constant - ∂ₜ(w) == a - w - w * n^2 + ν * ℒ(dX, w) + ∂ₜ(w) == a - w - w * n^2 + ν * L(dX, w) end # See Klausmeier Equation 2.b @@ -62,10 +68,12 @@ Phytodynamics = @decapode begin m::Constant ∂ₜ(n) == w * n^2 - m*n + Δ(n) -end # hide +end +nothing # hide ``` Now that we have our two component models, we can specify a means of composing them via a composition pattern. + ```@example DEC # Specify Composition compose_klausmeier = @relation () begin @@ -73,10 +81,11 @@ compose_klausmeier = @relation () begin hydro(N, W) end -to_graphviz(compose_klausmeier, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(compose_klausmeier) ``` We apply our composition pattern by plugging in component Decapodes, and specifying which internal quantities to share along edges. Decapodes are formalized via the field of Applied Category Theory. A practical consequence here is that we can view a Decapode as a sort of computation graph. + ```@example DEC # Apply Composition klausmeier_cospan = oapply(compose_klausmeier, @@ -87,34 +96,9 @@ to_graphviz(Klausmeier) ``` With our model now explicitly represented, we have everything we need to automatically generate simulation code. We could write this to an intermediate file and use it later, or we can go ahead and `eval`uate the code in this session. -```@example DEC -sim = eval(gensim(Klausmeier, dimension=1)) -``` - -We discretize our differential operators using the Discrete Exterior Calculus. The DEC is an elegant way of building up more complex differential operators from simpler ones. To demonstrate, we will define the Δ operator by building it up with matrix multiplication of simpler operators. Since most operators in the DEC are matrices, most simulations consist mainly of matrix-vector multiplications, and are thus very fast. We can also cache these - -If this code seems too low level, do not worry. Decapodes defines and caches for you many differential operators behind the scenes, so you do not have to worry about defining your own. ```@example DEC -function generate(sd, my_symbol; hodge=DiagonalHodge()) - lap_mat = hodge_star(1,sd) * d(0,sd) * inv_hodge_star(0,sd) * dual_derivative(0,sd) - dd_mat = dual_derivative(0,sd) - ih_mat = inv_hodge_star(0,sd) - h_mat = hodge_star(1,sd) - ih_dd_mat = ih_mat * dd_mat - op = @match my_symbol begin - :Δ => x -> begin - lap_mat * x - end - :ℒ => (x,y) -> begin - h_mat * ∧(0,1,sd, ih_dd_mat * y, x) - end - :^ => (x,y) -> begin - x .^ y - end - end - return (args...) -> op(args...) -end +sim = eval(gensim(Klausmeier, dimension=1)) ``` We now need a mesh to define our domain. In the 2D case, our CombinatorialSpaces library can read in arbitrary .OBJ files. In 1D, it is often simpler to just generate a mesh on the fly. Since we are running our physics on a circle - i.e periodic boundaries - we will use a simple function that generates it. @@ -134,12 +118,30 @@ function circle(n, c) subdivide_duals!(sd, Circumcenter()) s,sd end -s,sd = circle(7, 500) +s,sd = circle(9, 500) scatter(sd[:point]) ``` +We discretize our differential operators using the Discrete Exterior Calculus. The DEC is an elegant way of building up more complex differential operators from simpler ones. To demonstrate, we will define the Δ operator by building it up with matrix multiplication of simpler operators. Since most operators in the DEC are matrices, most simulations consist mainly of matrix-vector multiplications, and are thus very fast. + +If this code seems too low level, do not worry. Decapodes defines and caches for you many differential operators behind the scenes, so you do not have to worry about defining your own. + +```@example DEC +lap_mat = dec_hodge_star(1,sd) * dec_differential(0,sd) * dec_inv_hodge_star(0,sd) * dec_dual_derivative(0,sd) + +function generate(sd, my_symbol; hodge=DiagonalHodge()) + op = @match my_symbol begin + :Δ => x -> begin + lap_mat * x + end + end + return (args...) -> op(args...) +end +``` + Let's pass our mesh and methods of generating operators to our simulation code. + ```@example DEC # Instantiate Simulation fₘ = sim(sd, generate, DiagonalHodge()) @@ -148,6 +150,7 @@ fₘ = sim(sd, generate, DiagonalHodge()) With our simulation now ready, let's specify initial data to pass to it. We'll define them with plain Julia code. The most interesting parameter here is our "downhill gradient" `dX`. This parameter defines how steep our slope is. Since our mesh is a circle, and we are setting `dX` to a constant value, this means that "downhill" always points counter-clockwise. Essentially, this is an elegant way of encoding an infinite hill. + ```@example DEC # Define Initial Conditions n_dist = Normal(pi) @@ -166,11 +169,12 @@ cs_ps = (phyto_m = 0.45, ``` Let's execute our simulation. + ```@example DEC # Run Simulation -tₑ = 600.0 +tₑ = 300.0 prob = ODEProblem(fₘ, u₀, (0.0, tₑ), cs_ps) -sol = solve(prob, Tsit5()) +sol = solve(prob, Tsit5(), saveat=0.1, save_idxs=[:N, :W]) sol.retcode ``` @@ -182,7 +186,8 @@ Let's perform some basic visualization and analysis of our results to verify our n = sol(0).N nₑ = sol(tₑ).N w = sol(0).W -wₑ = sol(tₑ).W # hide +wₑ = sol(tₑ).W +nothing # hide ``` ```@example DEC @@ -205,6 +210,7 @@ save_dynamics(:N, 20, "klausmeier.gif") ![Klausmeier](klausmeier.gif) We can observe a few interesting phenomena that we wanted to capture: + - The vegetation density bands move uphill in traveling waves. - The leading edge of the waves is denser than the rest of the band. - Over time, the periodicity of the vegetation bands stabilizes. @@ -216,6 +222,7 @@ We can observe a few interesting phenomena that we wanted to capture: Due to the ease of composition of Decapodes, representing the Klausmeier model opens up many areas for future work. For example, we can now compose these dynamics with a model of temperature dynamics informed by the Budyko-Sellers model. We can take advantage of the fact that the Lie derivative generalizes partial derivatives, and model the flow of water according to any vector field. Or, we can extend this model by composing it with a model that can recreate the so-called "leopard pattern" of vegetation, such as an "Interaction-Dispersion" model of vegetation dynamics given by Lejeune et al\[4\]. ## References + \[1\] W. A. Macfadyen, “Vegetation Patterns in the Semi-Desert Plains of British Somaliland,” The Geographical Journal, vol. 116, no. 4/6, p. 199, Oct. 1950, doi: 10.2307/1789384. \[2\] C. A. Klausmeier, “Regular and Irregular Patterns in Semiarid Vegetation,” Science, vol. 284, no. 5421, pp. 1826–1828, Jun. 1999, doi: 10.1126/science.284.5421.1826. @@ -223,3 +230,7 @@ Due to the ease of composition of Decapodes, representing the Klausmeier model o \[3\] W. A. Macfadyen, “Soil and Vegetation in British Somaliland,” Nature, vol. 165, no. 4186, Art. no. 4186, Jan. 1950, doi: 10.1038/165121a0. \[4\] O. Lejeune and M. Tlidi, “A Model for the Explanation of Vegetation Stripes (Tiger Bush),” Journal of Vegetation Science, vol. 10, no. 2, pp. 201–208, 1999, doi: 10.2307/3237141. + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/klausmeier/somalia.png b/docs/src/klausmeier/somalia.png new file mode 100644 index 00000000..8ef3e87b Binary files /dev/null and b/docs/src/klausmeier/somalia.png differ diff --git a/docs/src/navier_stokes/ns.md b/docs/src/navier_stokes/ns.md index 32bf0361..7fb9e648 100644 --- a/docs/src/navier_stokes/ns.md +++ b/docs/src/navier_stokes/ns.md @@ -1,5 +1,10 @@ # Navier Stokes Vorticity Model +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + This is a discretization of the incompressible Navier Stokes equations using the Discrete Exterior Calculus. The formulations are based on those given by [Mohamed, Hirani, Samtaney](https://arxiv.org/abs/1508.01166) (in turn from [Marsden, Ratiu, Abraham](https://link.springer.com/book/10.1007/978-1-4612-1029-0)). @@ -87,4 +92,8 @@ Here is one set of results from using the inviscid Poisson formulation: These vortices should be stable so we should see the same periodic function for both lines here. The difference between the lines is the accumulated error. -![Azimuth Profile](azimuth.png) \ No newline at end of file +![Azimuth Profile](azimuth.png) + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/nhs/nhs.jl b/docs/src/nhs/nhs.jl new file mode 100644 index 00000000..b3fc2fb5 --- /dev/null +++ b/docs/src/nhs/nhs.jl @@ -0,0 +1,238 @@ +# AlgebraicJulia Dependencies +using Catlab +using CombinatorialSpaces +using Decapodes +using DiagrammaticEquations + +# External Dependencies +using CairoMakie +using ComponentArrays +using GeometryBasics: Point3 +using JLD2 +using LinearAlgebra +using MLStyle +using OrdinaryDiffEq +Point3D = Point3{Float64}; + +@info "Creating Decapode" +momentum = @decapode begin + (v,V)::DualForm1 + f::Form0 + uˢ::DualForm1 + ∂tuˢ::DualForm1 + p::DualForm0 + b::DualForm0 + ĝ::DualForm1 + Fᵥ::DualForm1 + StressDivergence::DualForm1 + + ∂ₜ(v) == + -ℒ₁(v,v) + 0.5*d(ι₁₁(v,v)) - + d(ι₁₁(v,V)) + ι₁₂(v,d(V)) + ι₁₂(V,d(v)) - + (f - ∘(d,⋆)(uˢ)) ∧ᵖᵈ₀₁ v - + d(p) + + b ∧ᵈᵈ₀₁ ĝ - + StressDivergence + + ∂tuˢ + + Fᵥ +end + +tracer_conservation = @decapode begin + (c,C,F,FluxDivergence)::DualForm0 + (v,V)::DualForm1 + + ∂ₜ(c) == + -1*ι₁₁(v,d(c)) - + ι₁₁(V,d(c)) - + ι₁₁(v,d(C)) - + FluxDivergence + + F +end + +equation_of_state = @decapode begin + (b,T,S)::DualForm0 + (g,α,β)::Constant + + b == g*(α*T - β*S) +end + +isotropic_diffusivity = @decapode begin + v::DualForm1 + c::DualForm0 + StressDivergence::DualForm1 + FluxDivergence::DualForm0 + (κ,nu)::Constant + + StressDivergence == nu*Δᵈ₁(v) + FluxDivergence == κ*Δᵈ₀(c) +end + +tracer_composition = @relation () begin + # "The turbulence closure selected by the user determines the form of ... diffusive flux divergence" + turbulence(FD,v,c) + + continuity(FD,v,c) +end + +isotropic_tracer = apex(oapply(tracer_composition, [ + Open(isotropic_diffusivity, [:FluxDivergence, :v, :c]), + Open(tracer_conservation, [:FluxDivergence, :v, :c])])) + + nonhydrostatic_composition = @relation () begin + # "The turbulence closure selected by the user determines the form of stress divergence" + # => Note that the StressDivergence term, SD, is shared by momentum and all the tracers. + momentum(V, v, b, SD) + + # "Both T and S obey the tracer conservation equation" + # => Temperature and Salinity both receive a copy of the tracer physics. + temperature(V, v, T, SD, nu) + salinity(V, v, S, SD, nu) + + # "Buoyancy is determined from a linear equation of state" + # => The b term in momentum is that described by the equation of state here. + eos(b, T, S) +end + +isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ + Open(momentum, [:V, :v, :b, :StressDivergence]), + Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_StressDivergence, :turbulence_nu]), + Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_StressDivergence, :turbulence_nu]), + Open(equation_of_state, [:b, :T, :S])])); + +@info "Creating Mesh" +# This is a torus with resolution of its dual mesh similar to that +# used by Oceananigans (explicitly represented as a torus, not as a +# square with periodic boundary conditions!) +download("https://cise.ufl.edu/~luke.morris/torus.obj", "torus.obj") +s = EmbeddedDeltaSet2D("torus.obj") +sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s) +subdivide_duals!(sd, Barycenter()) + +@info "Generating Simulation" +function generate(sd, my_symbol; hodge=GeometricHodge()) + op = @match my_symbol begin + _ => default_dec_matrix_generate(sd, my_symbol, hodge) + end + return op +end + +open("nhs_sim.jl", "w") do f + write(f, string(gensim(isotropic_nonhydrostatic_buoyancy))) +end +sim = include("nhs_sim.jl") + +@info "Creating Initial Conditions" +fₘ = sim(sd, generate) + +S = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) + 0.0 +end +T = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) + 0.0 +end +p = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) + 0.0 +end +f = zeros(nv(sd)) +Fₛ = zeros(ntriangles(sd)) +Fₜ = zeros(ntriangles(sd)) +Cₛ = zeros(ntriangles(sd)) +Cₜ = zeros(ntriangles(sd)) +V = zeros(ne(sd)) +v = rand(ne(sd)) * 1e-8 +ĝ = ♭(sd, DualVectorField(fill(Point3D(0,0,0), ntriangles(sd)))).data +Fᵥ = zeros(ne(sd)) +qₛ = zeros(ne(sd)) +qₜ = zeros(ne(sd)) +uˢ = zeros(ne(sd)) +∂tuˢ = zeros(ne(sd)) + +u₀ = ComponentArray( + v = v, + V = V, + momentum_f = f, + momentum_uˢ = uˢ, + momentum_∂tuˢ = ∂tuˢ, + momentum_p = p, + momentum_ĝ = ĝ, + momentum_Fᵥ = Fᵥ, + T = T, + temperature_continuity_C = Cₜ, + temperature_continuity_F = Fₜ, + S = S, + salinity_continuity_C = Cₛ, + salinity_continuity_F = Fₛ) + +gᶜ = 9.81 +α = 2e-3 +β = 5e-4 +constants_and_parameters = ( + temperature_turbulence_κ = 0.0, + salinity_turbulence_κ = 0.0, + nu = 1e-5, + eos_g = gᶜ, + eos_α = α, + eos_β = β); + +# 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, Vern7()) +soln.retcode != :Unstable || error("Solver was not stable") + +tₑ = 50 + +@info("Solving") +prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) +soln = solve(prob, Vern7(), force_dtmin=true, dtmax=0.2) +@show soln.retcode +@info("Done") + +ihs0 = dec_inv_hodge_star(Val{0}, sd, GeometricHodge()) +dd1 = dec_dual_derivative(1, sd) +♯_m = ♯_mat(sd, LLSDDSharp()) +using LinearAlgebra: norm +function vorticity(α) + ihs0*dd1*α +end +function speed(α) + norm.(♯_m * α) +end + +function save_vorticity(is_2d=false) + frames = 200 + time = Observable(0.0) + fig = Figure(title = @lift("Vorticity at $($time)")) + ax = is_2d ? + CairoMakie.Axis(fig[1,1]) : + LScene(fig[1,1], scenekw=(lights=[],)) + msh = CairoMakie.mesh!(ax, s, + color=@lift(vorticity(soln($time).v)), + colorrange=extrema(vorticity(soln(tₑ).v)).*.9, + colormap=:jet) + + Colorbar(fig[1,2], msh) + record(fig, "vorticity.gif", range(0.0, tₑ; length=frames); framerate = 20) do t + time[] = t + end +end +save_vorticity(false) + +function save_speed(is_2d=false) frames = 200 + time = Observable(0.0) + fig = Figure(title = @lift("Speed at $($time)")) + ax = is_2d ? + CairoMakie.Axis(fig[1,1]) : + LScene(fig[1,1], scenekw=(lights=[],)) + msh = CairoMakie.scatter!(ax, sd[sd[:tri_center], :dual_point], + color=@lift(speed(soln($time).v)), + colorrange=extrema(speed(soln(tₑ).v)).*.9, + colormap=:jet, + markersize=5) + + Colorbar(fig[1,2], msh) + record(fig, "speed.gif", range(0.0, tₑ; length=frames); framerate = 20) do t + time[] = t + end +end +save_speed(false) diff --git a/docs/src/nhs/nhs.sh b/docs/src/nhs/nhs.sh new file mode 100644 index 00000000..a4b30a12 --- /dev/null +++ b/docs/src/nhs/nhs.sh @@ -0,0 +1,27 @@ +#!/bin/bash +#SBATCH --job-name=NHS # Job name +#SBATCH --mail-type=ALL # Mail events (NONE, BEGIN, END, FAIL, ALL) +#SBATCH --mail-user=luke.morris@ufl.edu # Where to send mail +#SBATCH --ntasks=1 # Run on a single CPU +#SBATCH --cpus-per-task=16 # Number of cores you can use for threads. +#SBATCH --mem=32gb # Job memory request +#SBATCH --time=08:00:00 # Time limit hrs:min:sec +#SBATCH --output=nhs_%j.log # Standard output and error log +#SBATCH --export=ALL # Give the compute node your environment variables (e.g. JULIA_PKG_PRECOMPILE_AUTO) +pwd; hostname; date + +module load julia + +mkdir -p "/blue/fairbanksj/fairbanksj/jldepot" +export JULIA_DEPOT_PATH="/blue/fairbanksj/fairbanksj/jldepot" +echo "JULIA_DEPOT_PATH:" +echo "$JULIA_DEPOT_PATH" + +echo "Launching script" +date + +julia --threads=auto --proj=. -e 'using Pkg; Pkg.instantiate(); Pkg.precompile()' +julia --threads=auto --proj=. ./nhs.jl + +date +echo "Exiting script" diff --git a/docs/src/nhs.md b/docs/src/nhs/nhs_full.md similarity index 88% rename from docs/src/nhs.md rename to docs/src/nhs/nhs_full.md index d01d0172..1fba9d20 100644 --- a/docs/src/nhs.md +++ b/docs/src/nhs/nhs_full.md @@ -1,6 +1,6 @@ # Implement Oceananigans.jl's NonhydrostaticModel in the Discrete Exterior Calculus -Let's use Decapodes to implement the [NonhydrostaticModel](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/) from Oceananigans.jl. We will take the opportunity to demonstrate how we can use our "algebra of model compositions" to encode certain guarantees on the models we generate. We will use the [2D Turbulence](https://clima.github.io/OceananigansDocumentation/stable/generated/two_dimensional_turbulence/) as a guiding example, and use only equations found in the Oceananigans docs to construct our model. +Let's use Decapodes to implement the [NonhydrostaticModel](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/) from Oceananigans.jl. We will take the opportunity to demonstrate how we can use our "algebra of model compositions" to encode certain guarantees on the models we generate. We will use the [2D Turbulence](https://clima.github.io/OceananigansDocumentation/stable/literated/two_dimensional_turbulence/) as a guiding example, and use only equations found in the Oceananigans docs to construct our model. ```@example DEC # AlgebraicJulia Dependencies @@ -18,11 +18,12 @@ using LinearAlgebra using MLStyle using OrdinaryDiffEq Point3D = Point3{Float64}; +nothing # hide ``` -## Specify our models. +## Specify our models -This is [Equation 1: "The momentum conservation equation"](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-momentum-conservation-equation). This is the first formulation of mutual advection (of v along V, and V along v) that we could find in the exterior calculus. +This is [Equation 1: "The momentum conservation equation"](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-momentum-conservation-equation). This is the first formulation of mutual advection (of `v` along `V`, and `V` along `v`) that we could find in the exterior calculus. ```@example DEC momentum = @decapode begin @@ -49,9 +50,10 @@ end to_graphviz(momentum) ``` -Why did we write "StressDivergence" instead of ∇⋅τ, as in the linked equation? According to [this docs page](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/), the user makes a selection of what model to insert in place of the term ∇⋅τ. For example, in [the isotropic case](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity), Oceananigans.jl replaces this term with: ∇⋅τ = nuΔv. Thus, we write StressDivergence, and replace this term with a choice of "turbulence closure" model. Using the "constant isotropic diffusivity" case, we can operate purely in terms of scalar-valued forms. +Why did we write "StressDivergence" instead of ∇⋅τ, as in the linked equation? According to [this docs page](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/), the user makes a selection of what model to insert in place of the term ∇⋅τ. For example, in [the isotropic case](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity), Oceananigans.jl replaces this term with: ∇⋅τ = *ν*Δv. Thus, we write StressDivergence, and replace this term with a choice of "turbulence closure" model. Using the "constant isotropic diffusivity" case, we can operate purely in terms of scalar-valued forms. This is [Equation 2: "The tracer conservation equation"](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-tracer-conservation-equation). + ```@example DEC tracer_conservation = @decapode begin (c,C,F,FluxDivergence)::DualForm0 @@ -68,6 +70,7 @@ to_graphviz(tracer_conservation) ``` This is [Equation 2: "Linear equation of state"](https://clima.github.io/OceananigansDocumentation/stable/physics/buoyancy_and_equations_of_state/#Linear-equation-of-state) of seawater buoyancy. + ```@example DEC equation_of_state = @decapode begin (b,T,S)::DualForm0 @@ -79,6 +82,7 @@ to_graphviz(equation_of_state) ``` This is [Equation 2: "Constant isotropic diffusivity"](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity). + ```@example DEC isotropic_diffusivity = @decapode begin v::DualForm1 @@ -95,11 +99,12 @@ to_graphviz(isotropic_diffusivity) ## Compatibility Guarantees via Operadic Composition -Decapodes composition is formally known as an "operad algebra". That means that we don't have to encode our composition in a single UWD and then apply it. Rather, we can define several UWDs, compose those, and then apply those. Of course, since the output of oapply is another Decapode, we could perform an intermediate oapply, if that is convenient. +Decapodes composition is formally known as an "operad algebra". That means that we don't have to encode our composition in a single undirected wiring diagram (UWD) and then apply it. Rather, we can define several UWDs, compose those, and then apply those. Of course, since the output of oapply is another Decapode, we could perform an intermediate oapply, if that is convenient. Besides it being convenient to break apart large UWDs into component UWDs, this hierarchical composition can enforce rules on our physical quantities. For example: + 1. We want all the tracers (salinity, temperature, etc.) in our physics to obey the same conservation equation. 2. We want them to obey the same "turbulence closure", which affects their flux-divergence term. 3. At the same time, a choice of turbulence closure doesn't just affect (each of) the flux-divergence terms, it also constrains which stress-divergence is physically valid in the momentum equation. @@ -107,6 +112,7 @@ For example: We will use our operad algebra to guarantee model compatibility and physical consistency, guarantees that would be burdensome to fit into a one-off type system. Here, we specify the equations that any tracer obeys: + ```@example DEC tracer_composition = @relation () begin # "The turbulence closure selected by the user determines the form of ... diffusive flux divergence" @@ -114,10 +120,11 @@ tracer_composition = @relation () begin continuity(FD,v,c) end -to_graphviz(tracer_composition, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(tracer_composition) ``` - Let's "lock in" isotropic diffusivity by doing an intermediate oapply. +Let's "lock in" isotropic diffusivity by doing an intermediate oapply. + ```@example DEC isotropic_tracer = apex(oapply(tracer_composition, [ Open(isotropic_diffusivity, [:FluxDivergence, :v, :c]), @@ -142,10 +149,9 @@ nonhydrostatic_composition = @relation () begin # => The b term in momentum is that described by the equation of state here. eos(b, T, S) end -to_graphviz(nonhydrostatic_composition, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(nonhydrostatic_composition) ``` - ```@example DEC isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ Open(momentum, [:V, :v, :b, :StressDivergence]), @@ -170,16 +176,16 @@ subdivide_duals!(sd, Barycenter()) function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin - _ => default_dec_matrix_generate(sd, my_symbol, hodge) + _ => default_dec_matrix_generate(sd, my_symbol, hodge) # TODO: Don't do this! end - return (args...) -> op(args...) + return op end open("nhs.jl", "w") do f - write(f, string(gensim(expand_operators(isotropic_nonhydrostatic_buoyancy)))) + write(f, string(gensim(isotropic_nonhydrostatic_buoyancy))) end sim = include("nhs.jl") -fₘ = sim(sd, generate) +fₘ = sim(sd, generate) # TODO: Slow because dual operators take too long S = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) 0.0 @@ -237,17 +243,17 @@ constants_and_parameters = ( We specified our physics, our mesh, and our initial conditions. We have everything we need to execute the simulation. ```@example DEC -tₑ = 50 - # 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, Vern7()) soln.retcode != :Unstable || error("Solver was not stable") +tₑ = 50 + @info("Solving") prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) -soln = solve(prob, Vern7(), force_dtmin=true, dtmax=0.2, progress=true,progress_steps=1) +soln = solve(prob, Vern7(), force_dtmin=true, dtmax=0.2) @show soln.retcode @info("Done") ``` @@ -310,4 +316,3 @@ save_speed(false) ![Vorticity](vorticity.gif) ![Speed](speed.gif) - diff --git a/docs/src/nhs/nhs_lite.md b/docs/src/nhs/nhs_lite.md new file mode 100644 index 00000000..9644cc6b --- /dev/null +++ b/docs/src/nhs/nhs_lite.md @@ -0,0 +1,203 @@ +# Implement Oceananigans.jl's NonhydrostaticModel in the Discrete Exterior Calculus + +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + +Let's use Decapodes to implement the [NonhydrostaticModel](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/) from Oceananigans.jl. We will take the opportunity to demonstrate how we can use our "algebra of model compositions" to encode certain guarantees on the models we generate. We will use the [2D Turbulence](https://clima.github.io/OceananigansDocumentation/stable/literated/two_dimensional_turbulence/) as a guiding example, and use only equations found in the Oceananigans docs to construct our model. + +The full code that generated these results is available in a [julia script](nhs.jl). + +```@example DEC +# AlgebraicJulia Dependencies +using Catlab +using CombinatorialSpaces +using Decapodes +using DiagrammaticEquations + +# External Dependencies +using CairoMakie +using ComponentArrays +using Downloads +using GeometryBasics: Point3 +using JLD2 +using LinearAlgebra +using MLStyle +using OrdinaryDiffEq +Point3D = Point3{Float64}; +nothing # hide +``` + +## Specify our models + +This is [Equation 1: "The momentum conservation equation"](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-momentum-conservation-equation). This is the first formulation of mutual advection (of `v` along `V`, and `V` along `v`) that we could find in the exterior calculus. + +```@example DEC +momentum = @decapode begin + (v,V)::DualForm1 + f::Form0 + uˢ::DualForm1 + ∂tuˢ::DualForm1 + p::DualForm0 + b::DualForm0 + ĝ::DualForm1 + Fᵥ::DualForm1 + StressDivergence::DualForm1 + + ∂ₜ(v) == + -ℒ₁(v,v) + 0.5*d(ι₁₁(v,v)) - + d(ι₁₁(v,V)) + ι₁₂(v,d(V)) + ι₁₂(V,d(v)) - + (f - ∘(d,⋆)(uˢ)) ∧ᵖᵈ₀₁ v - + d(p) + + b ∧ᵈᵈ₀₁ ĝ - + StressDivergence + + ∂tuˢ + + Fᵥ +end +to_graphviz(momentum) +``` + +Why did we write "StressDivergence" instead of ∇⋅τ, as in the linked equation? According to [this docs page](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/), the user makes a selection of what model to insert in place of the term ∇⋅τ. For example, in [the isotropic case](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity), Oceananigans.jl replaces this term with: ∇⋅τ = *ν*Δv. Thus, we write StressDivergence, and replace this term with a choice of "turbulence closure" model. Using the "constant isotropic diffusivity" case, we can operate purely in terms of scalar-valued forms. + +This is [Equation 2: "The tracer conservation equation"](https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-tracer-conservation-equation). + +```@example DEC +tracer_conservation = @decapode begin + (c,C,F,FluxDivergence)::DualForm0 + (v,V)::DualForm1 + + ∂ₜ(c) == + -1*ι₁₁(v,d(c)) - + ι₁₁(V,d(c)) - + ι₁₁(v,d(C)) - + FluxDivergence + + F +end +to_graphviz(tracer_conservation) +``` + +This is [Equation 2: "Linear equation of state"](https://clima.github.io/OceananigansDocumentation/stable/physics/buoyancy_and_equations_of_state/#Linear-equation-of-state) of seawater buoyancy. + +```@example DEC +equation_of_state = @decapode begin + (b,T,S)::DualForm0 + (g,α,β)::Constant + + b == g*(α*T - β*S) +end +to_graphviz(equation_of_state) +``` + +This is [Equation 2: "Constant isotropic diffusivity"](https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity). + +```@example DEC +isotropic_diffusivity = @decapode begin + v::DualForm1 + c::DualForm0 + StressDivergence::DualForm1 + FluxDivergence::DualForm0 + (κ,nu)::Constant + + StressDivergence == nu*Δᵈ₁(v) + FluxDivergence == κ*Δᵈ₀(c) +end +to_graphviz(isotropic_diffusivity) +``` + +## Compatibility Guarantees via Operadic Composition + +Decapodes composition is formally known as an "operad algebra". That means that we don't have to encode our composition in a single undirected wiring diagram (UWD) and then apply it. Rather, we can define several UWDs, compose those, and then apply those. Of course, since the output of oapply is another Decapode, we could perform an intermediate oapply, if that is convenient. + +Besides it being convenient to break apart large UWDs into component UWDs, this hierarchical composition can enforce rules on our physical quantities. + +For example: + +1. We want all the tracers (salinity, temperature, etc.) in our physics to obey the same conservation equation. +2. We want them to obey the same "turbulence closure", which affects their flux-divergence term. +3. At the same time, a choice of turbulence closure doesn't just affect (each of) the flux-divergence terms, it also constrains which stress-divergence is physically valid in the momentum equation. + +We will use our operad algebra to guarantee model compatibility and physical consistency, guarantees that would be burdensome to fit into a one-off type system. + +Here, we specify the equations that any tracer obeys: + +```@example DEC +tracer_composition = @relation () begin + # "The turbulence closure selected by the user determines the form of ... diffusive flux divergence" + turbulence(FD,v,c) + + continuity(FD,v,c) +end +draw_composition(tracer_composition) +``` + +Let's "lock in" isotropic diffusivity by doing an intermediate oapply. + +```@example DEC +isotropic_tracer = apex(oapply(tracer_composition, [ + Open(isotropic_diffusivity, [:FluxDivergence, :v, :c]), + Open(tracer_conservation, [:FluxDivergence, :v, :c])])) +to_graphviz(isotropic_tracer) +``` + +Let's use this building-block tracer physics at the next level. The quotes that appear in this composition diagram appear directly in the Oceananigans.jl docs. + +```@example DEC +nonhydrostatic_composition = @relation () begin + # "The turbulence closure selected by the user determines the form of stress divergence" + # => Note that the StressDivergence term, SD, is shared by momentum and all the tracers. + momentum(V, v, b, SD) + + # "Both T and S obey the tracer conservation equation" + # => Temperature and Salinity both receive a copy of the tracer physics. + temperature(V, v, T, SD, nu) + salinity(V, v, S, SD, nu) + + # "Buoyancy is determined from a linear equation of state" + # => The b term in momentum is that described by the equation of state here. + eos(b, T, S) +end +draw_composition(nonhydrostatic_composition) +``` + +```@example DEC +isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ + Open(momentum, [:V, :v, :b, :StressDivergence]), + Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_StressDivergence, :turbulence_nu]), + Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_StressDivergence, :turbulence_nu]), + Open(equation_of_state, [:b, :T, :S])])); +to_graphviz(isotropic_nonhydrostatic_buoyancy) +``` + +## Our Mesh + +We execute these dynamics on the torus explicitly, instead of using a square with periodic boundary conditions. + +```@example DEC +# This is a torus with resolution of its dual mesh similar to that +# used by Oceananigans (explicitly represented as a torus, not as a +# square with periodic boundary conditions!) +Downloads.download("https://cise.ufl.edu/~luke.morris/torus.obj", "torus.obj") +s = EmbeddedDeltaSet2D("torus.obj") +sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s) +subdivide_duals!(sd, Barycenter()) +fig = Figure() # hide +ax = CairoMakie.Axis(fig[1,1], aspect=1) # hide +wf = wireframe!(ax, s; linewidth=1) # hide +save("NHS_mesh.png", fig) # hide +nothing # hide +``` + +!["NHS_torus"](NHS_mesh.png) + +## Results + +In the DEC, vorticity is encoded with `d⋆`, and speed can be encoded with `norm ♯`. We can use our operators from CombinatorialSpaces.jl to create our GIFs. + +![Vorticity](vorticity.gif) + +![Speed](speed.gif) + +```@example INFO +DocInfo.get_report(info) # hide +``` diff --git a/docs/src/nhs/speed.gif b/docs/src/nhs/speed.gif new file mode 100644 index 00000000..025802c9 Binary files /dev/null and b/docs/src/nhs/speed.gif differ diff --git a/docs/src/nhs/vorticity.gif b/docs/src/nhs/vorticity.gif new file mode 100644 index 00000000..ec784f6f Binary files /dev/null and b/docs/src/nhs/vorticity.gif differ diff --git a/docs/src/overview.md b/docs/src/overview/overview.md similarity index 73% rename from docs/src/overview.md rename to docs/src/overview/overview.md index 9865331d..abb44eb7 100644 --- a/docs/src/overview.md +++ b/docs/src/overview/overview.md @@ -1,32 +1,34 @@ # Introduction to Decapodes +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + Discrete Exterior Calculus Applied to Partial and Ordinary Differential Equations (Decapodes) is a diagrammatic language used to express systems of -ordinary and partial differential equations. The Decapode provides a visual +ordinary and partial differential equations. Decapodes provides a visual framework for understanding the coupling between variables within a PDE or ODE system, and a combinatorial data structure for working with them. Below, we provide a high-level overview of how Decapodes can be generated and interpreted. ## Your First Decapode -We begin with the most basic Decapode, one which only includes a single -variable. In the Decapode graphical paradigm, nodes represent variables and -arrows represent operators which relate variables to each other. Since the -Decapode applies this diagrammatic language specifically to the Discrete +In the Decapodes graphical paradigm, nodes represent variables and +arrows represent operators which relate variables to each other. Since +Decapodes applies this diagrammatic language specifically to the Discrete Exterior Calculus (DEC), variables are typed by the dimension and orientation of the information they contain. So a variable of type `Form0` will be the -0-dimensional data points on some space, or in a discrete context, the -values defined on points of a mesh. Similarly, `Form1` will be values -stored on edges of the mesh, and `Form2` will be values stored on the -surfaces of the mesh. Below, we provide a very simple Decapode with just a -single variable `C`. -and define a convenience function for visualization in later examples. +0-dimensional data points defined the vertices of a mesh. Similarly, `Form1` will be values +stored on edges of the mesh and `Form2` will be values stored on the +surfaces of the mesh. + +Below, we provide a Decapode with just a single variable `C` and display it. ```@example DEC -using DiagrammaticEquations -using DiagrammaticEquations.Deca +using Catlab using Decapodes -using Catlab, Catlab.Graphics +using DiagrammaticEquations Variable = @decapode begin C::Form0 @@ -36,7 +38,7 @@ to_graphviz(Variable) ``` The resulting diagram contains a single node, showing the single variable in -this system. We can then add a second variable: +this system. We can add a second variable: ```@example DEC TwoVariables = @decapode begin @@ -47,24 +49,24 @@ end; to_graphviz(TwoVariables) ``` -And then can add some relationship between them. In this case, we make an -equation which states that `dC` is the discrete derivative of `C`: +We can also add a relationship between them. In this case, we make an +equation which states that `dC` is the derivative of `C`: ```@example DEC Equation = @decapode begin C::Form0 dC::Form1 - dC == d₀(C) + dC == d(C) end; to_graphviz(Equation) ``` Here, the two nodes represent the two variables, and the arrow between them -shows how they are related by the discrete derivative. +shows how they are related by the derivative. -## A Little More Complicated +## A Little More Complicated Now that we've seen how to construct a simple equation, it's time to move on to some actual PDE systems! One classic PDE example is the diffusion equation. @@ -78,36 +80,35 @@ Diffusion = @decapode begin # Fick's first law ϕ == k(d₀(C)) + # Diffusion equation Ċ == ⋆₀⁻¹(dual_d₁(⋆₁(ϕ))) ∂ₜ(C) == Ċ + end; to_graphviz(Diffusion) ``` The resulting Decapode shows the relationships between the three variables with -the triangle diagram. Note that these diagrams are automatically layed-out by Graphviz. +the triangle diagram. Note that these diagrams are automatically layed-out by [Graphviz](https://graphviz.org/). ## Bring in the Dynamics Now that we have a reasonably complex PDE, we can demonstrate some of the developed tooling for actually solving the PDE. Currently, the tooling will automatically generate an explicit method for solving the system (using -DifferentialEquations.jl to handle time-stepping and instability detection). +[DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl?tab=readme-ov-file) to handle time-stepping and instability detection). -We begin this process by importing a mesh. The mesh has been pre-generated -within CombinatorialSpaces, and is generated such that it has periodic boundary -conditions. We will also upload a non-periodic mesh for the sake of +`Torus_30x10` is a default mesh that is downloaded via `Artifacts.jl` when a user installs [CombinatorialSpaces.jl](https://github.com/AlgebraicJulia/CombinatorialSpaces.jl). If we wanted, we could also instantiate any `.obj` file of triangulated faces as a simplicial set although we do not here. + +We will also upload a non-periodic mesh for the sake of visualization, as well as a mapping between the points on the periodic and non-periodic meshes. -See CombinatorialSpaces.jl for mesh construction and importing utilities. - ```@example DEC -using Catlab.CategoricalAlgebra -using CombinatorialSpaces, CombinatorialSpaces.DiscreteExteriorCalculus using CairoMakie +using CombinatorialSpaces plot_mesh = loadmesh(Rectangle_30x10()) periodic_mesh = loadmesh(Torus_30x10()) @@ -119,17 +120,12 @@ wireframe!(ax, plot_mesh) fig ``` -With the mesh uploaded, we also need to convert the Decapode into something -which can be scheduled with explicit time stepping. In order to do this, we take -every variable which is the time derivative of another variable and trace back -the operations needed to compute this. This process essentially generates a computation -graph in the form of a directed wiring diagram. +Now we create a function which links the names of functions used in the Decapode to their implementations. Note that many DEC operators are already +defined for you. -Since our diagram is already defined, we just need to define a function which implements each of -these symbolic operators and pass them to a scheduler for generating the -function. +As an example, we chose to define `k` as a function that multiplies an input by `0.05`. We could have alternately chosen to represent `k` as a `Constant` that we multiply by in the Decapode itself. -Note that we chose to define `k` as a function that multiplies by a value `k`. We could have alternately chosen to represent `k` as a Constant that we multiply by in the Decapode itself. +We then compile the simulation by using `gen_sim` and create functional simulation by calling the evaluated `sim` with the mesh and our `generate` function. ```@example DEC using MLStyle @@ -139,17 +135,16 @@ function generate(sd, my_symbol; hodge=DiagonalHodge()) :k => x -> 0.05*x x => error("Unmatched operator $my_symbol") end - return (args...) -> op(args...) + return op end -``` - -Next, we generate the simulation function using `gen_sim` and set up our -initial conditions for this problem. -```@example DEC sim = eval(gensim(Diffusion)) fₘ = sim(periodic_mesh, generate, DiagonalHodge()) +``` + +We go ahead and set up our initial conditions for this problem. In this case we generate a Gaussian and apply it to our mesh. +```@example DEC using Distributions c_dist = MvNormal([7, 5], [1.5, 1.5]) c = [pdf(c_dist, [p[1], p[2]]) for p in periodic_mesh[:point]] @@ -160,7 +155,7 @@ mesh!(ax, plot_mesh; color=c[point_map]) fig ``` -Finally, we solve this PDE problem using the `Tsit5()` solver and generate an animation of the result! +Finally, we solve this PDE problem using the `Tsit5()` solver provided by DifferentialEquations.jl. ```@example DEC using LinearAlgebra @@ -171,7 +166,12 @@ u₀ = ComponentArray(C=c) prob = ODEProblem(fₘ, u₀, (0.0, 100.0)) sol = solve(prob, Tsit5()); +sol.retcode +``` +Now that the simulation has succeeded we can plot out our results with [CairoMakie.jl](https://github.com/MakieOrg/Makie.jl). + +```@example DEC # Plot the result times = range(0.0, 100.0, length=150) colors = [sol(t).C[point_map] for t in times] @@ -189,12 +189,12 @@ record(fig, "diffusion.gif", range(0.0, 100.0; length=150); framerate = 30) do t end ``` -![](diffusion.gif) +![Your first Decapode!](diffusion.gif) ## Merging Multiple Physics Now that we've seen the basic pipeline, it's time for a more complex example -that demonstrates some of the benefits reaped from using Catlab.jl as the +that demonstrates some of the benefits reaped from using [Catlab.jl](https://github.com/AlgebraicJulia/Catlab.jl) as the backend to our data structures. In this example, we will take two separate physics (diffusion and advection), and combine them together using a higher-level composition pattern. @@ -217,6 +217,7 @@ Advection = @decapode begin C::Form0 ϕ::Form1 V::Form1 + ϕ == ∧₀₁(C,V) end @@ -228,19 +229,25 @@ Superposition = @decapode begin Ċ == ⋆₀⁻¹(dual_d₁(⋆₁(ϕ))) ∂ₜ(C) == Ċ end -true # hide +nothing # hide ``` +The diffusion Decapode. + ```@example DEC -to_graphviz(Diffusion) +to_graphviz(Diffusion) # hide ``` +The advection Decapode. + ```@example DEC -to_graphviz(Advection) +to_graphviz(Advection) # hide ``` +And the superposition Decapode. + ```@example DEC -to_graphviz(Superposition) +to_graphviz(Superposition) # hide ``` Next, we define the pattern of composition which we want to compose these @@ -255,7 +262,7 @@ compose_diff_adv = @relation (C, V) begin superposition(ϕ₁, ϕ₂, ϕ, C) end -to_graphviz(compose_diff_adv, box_labels=:name, junction_labels=:variable, prog="circo") +draw_composition(compose_diff_adv) ``` After this, the physics can be composed as follows: @@ -316,9 +323,6 @@ using MLStyle function generate(sd, my_symbol; hodge=DiagonalHodge()) op = @match my_symbol begin :k => x -> 0.05*x - :∧₀₁ => (x,y) -> begin - ∧(Tuple{0,1}, sd, x,y) - end x => error("Unmatched operator $my_symbol") end return (args...) -> op(args...) @@ -352,5 +356,8 @@ record(fig, "diff_adv.gif", range(0.0, 100.0; length=150); framerate = 30) do t end ``` -![](diff_adv.gif) +![Your first composed Decapode!](diff_adv.gif) + +```@example INFO +DocInfo.get_report(info) # hide ``` diff --git a/docs/src/poiseuille.md b/docs/src/poiseuille/poiseuille.md similarity index 78% rename from docs/src/poiseuille.md rename to docs/src/poiseuille/poiseuille.md index 53997083..2e2e1396 100644 --- a/docs/src/poiseuille.md +++ b/docs/src/poiseuille/poiseuille.md @@ -1,31 +1,34 @@ # Poissuille Flow for Fluid Mechanics -When modeling a fluid flowing in a pipe, one can ignore the multidimensional structure of the pipe and approximate the system as a 1 dimensional flow along the pipe. The noslip boundary condition and the geometry of the pipe enter a 1D equation in the form of a resistance term. +```@setup INFO +include(joinpath(Base.@__DIR__, "..", "..", "docinfo.jl")) +info = DocInfo.Info() +``` + +When modeling a fluid flowing in a pipe, one can ignore the multidimensional structure of the pipe and approximate the system as a 1 dimensional flow along the pipe. The no-slip boundary condition and the geometry of the pipe enter a 1D equation in the form of a resistance term. ```@example Poiseuille using Catlab using CombinatorialSpaces -using CombinatorialSpaces.ExteriorCalculus -using CombinatorialSpaces.DiscreteExteriorCalculus: ∧ using DiagrammaticEquations -using DiagrammaticEquations.Deca using Decapodes # Julia community libraries using CairoMakie -using GeometryBasics +using GeometryBasics: Point3 using LinearAlgebra using OrdinaryDiffEq +Point3D = Point3{Float64} +nothing # hide ``` -# Creating the Poiseuille Equations +## Creating the Poiseuille Equations The `@decapode` macro creates the data structure representing the equations of Poiseuille flow. The first block declares variables, the second block defines intermediate terms and the last block is the core equation. -```@example Poiseuille -# μ̃ = negative viscosity per unit area -# R = drag of pipe boundary +For these physics, `μ̃` represents the negative viscosity per unit area while `R` represents the drag of the pipe boundary. +```@example Poiseuille Poise = @decapode begin P::Form0 q::Form1 @@ -49,13 +52,13 @@ resolve_overloads!(Poise, op1_res_rules_1D, op2_res_rules_1D) to_graphviz(Poise) ``` -# Defining the Semantics +## Defining the Semantics -In order to solve our equations, we will need numerical linear operators that give meaning to our symbolic operators. The `generate` function below assigns the necessary matrices as definitions for the symbols. In order to define the viscosity effect correctly we have to identify boundary edges and apply a mask. This is because the DEC has discrete dual cells at the boundaries that need to be handled specially for the viscosity term. We found empirically that if you allow nonzero viscosity at the boundary edges, the flows at the boundaries will be incorrect. +In order to solve our equations, we will need numerical linear operators that give meaning to our symbolic operators. The `generate` function below assigns the necessary matrices as definitions for the symbols. In order to define the viscosity effect correctly we have to identify boundary edges and apply a mask. This is because the DEC has discrete dual cells at the boundaries that need to be handled specially for the viscosity term. We found empirically that if you allow nonzero viscosity at the boundary edges, the flows at the boundaries will be incorrect. You can find the file for boundary conditions [here](../boundary_helpers.jl). ```@example Poiseuille using MLStyle -include("../../examples/boundary_helpers.jl") +include("../boundary_helpers.jl") function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin @@ -63,9 +66,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) x[boundary_edges(sd)] .= 0 x end - :∧₀₁ => (x,y) -> begin - ∧(0,1, sd, x,y) - end :∂ρ => ρ -> begin ρ[1] = 0 ρ[end] = 0 @@ -73,16 +73,15 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) end x => error("Unmatched operator $my_symbol") end - return (args...) -> op(args...) + return op end ``` ## A Single Pipe Segment -We create a mesh with one pipe segment to see if we get the right answer. This simulation can be validated with the Poiseuille equation for a single pipe. First we create the mesh. +This simulation can be validated with the Poiseuille equation for a single pipe. First we create a mesh with one pipe segment. ```@example Poiseuille -Point3D = Point3{Float64} s = EmbeddedDeltaSet1D{Bool,Point3D}() add_vertices!(s, 2, point=[Point3D(-1, 0, 0), Point3D(+1, 0, 0)]) add_edge!(s, 1, 2, edge_orientation=true) @@ -103,8 +102,8 @@ P = [10.0, 5.0] u = ComponentArray(q=q,P=P) params = (k = -0.01, μ̃ = 0.5, R=0.005) prob = ODEProblem(fₘ, u, (0.0, 10000.0), params) -sol = solve(prob, Tsit5()) -sol.u +soln = solve(prob, Tsit5()) +soln.u ``` ## A Linear Pipe with Multiple Segments @@ -122,12 +121,12 @@ function linear_pipe(n::Int) end sd = linear_pipe(10) -true # hide +nothing # hide ``` Then we solve the equation. Notice that the equilibrium flow is constant down the length of the pipe. This must be true because of conservation of mass. The segments are all the same length and the total flow in must equal the total flow out of each segment. -Note that we do not generate new simulation code for Poiseuille flow with `gensim` again. We provide our new mesh so that our discrete differential operators can be instantiated. +Note that we do not generate new simulation code for Poiseuille flow with `gensim` again. We simply need to provide our new mesh so that our discrete differential operators can be re-instantiated. ```@example Poiseuille fₘ = sim(sd, generate) @@ -142,7 +141,7 @@ sol.u ## A Distribution Network -To model a distribution network such as residential drinking water or natural gas, we will build a binary tree of pipes that at each junction have a bifurcation into two pipes. We expect that the flow will be divided by two at each level of the tree. First we make the mesh. +To model a distribution network, such as residential drinking water system, we will build a binary tree of pipes that at each junction have a bifurcation into two pipes. We expect that the flow will be divided by two at each level of the tree. First we make the mesh. ```@example Poiseuille function binary_pipe(depth::Int) @@ -163,7 +162,7 @@ function binary_pipe(depth::Int) sd end sd = binary_pipe(2) -true # hide +nothing # hide ``` Then we solve the equations. @@ -183,10 +182,9 @@ sol.u Decapodes really shines when you want to extend or refine your physics. We will change our physics by adding in a term for density of the material and the corresponding changes in pressure. This is not the only formulation for including a dynamic pressure effect into this system. If you can think of a better way to include this effect, we invite you to try it as an exercise! -Because the pressure is no longer being supplied as a parameter of the system controlled by the operators, we need to introduce a density term and a boundary condition for that density. In this system you can think of forcing a prescribed amount of material per unit time through the openings of the pipe and allowing the flow (q) and the pressure (P) to fluctuate. Before we were enforcing a fixed pressure gradient and and letting the flow fluctuate to achieve equilibrium. In the prior model, we were not accounting for the amount of material that had to flow in order to achieve that (flow, pressure) combination. +Because the pressure is no longer being supplied as a parameter of the system controlled by the operators, we need to introduce a density term and a boundary condition for that density. In this system you can think of forcing a prescribed amount of material per unit time through the openings of the pipe and allowing the flow `q` and the pressure `P` to fluctuate. Before we were enforcing a fixed pressure gradient and and letting the flow fluctuate to achieve equilibrium. In the prior model, we were not accounting for the amount of material that had to flow in order to achieve that (flow, pressure) combination. - -The Decapode can be visualized with graphviz, note that the boundary conditions are explicitly represented in the Decapode as operators that implement a masking operation. This is not consistent with the Diagrammatic Equations in Physics paper [PBHF22]. This approach is more directly tied to the computational method and will eventually be replaced with one based on morphisms of diagrams. +The Decapode can be visualized with [Graphviz](https://graphviz.org/), note that the boundary conditions are explicitly represented in the Decapode as operators that implement a masking operation. This is not consistent with the [Diagrammatic Equations in Physics](https://www.aimspress.com/article/id/62989382ba35de155149669f) paper. This approach is more directly tied to the computational method and will eventually be replaced with one based on morphisms of diagrams. ```@example Poiseuille # μ̃ = negative viscosity per unit area @@ -236,3 +234,7 @@ sol.u ``` Notice that the solution contains both a vector of flows and a vector of pressures. + +```@example INFO +DocInfo.get_report(info) # hide +``` 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/src/canon/Chemistry.jl b/src/canon/Chemistry.jl index 8e9dc6ad..95cc9875 100644 --- a/src/canon/Chemistry.jl +++ b/src/canon/Chemistry.jl @@ -6,7 +6,7 @@ using ..Canon using Markdown @docapode("Brusselator" - ,"https://en.wikipedia.org/wiki/brusselator" + ,"https://en.wikipedia.org/wiki/Brusselator" ,"A model of reaction-diffusion for an oscillatory chemical system." ,brusselator ,begin diff --git a/src/canon/Physics.jl b/src/canon/Physics.jl index 44afb8b3..2a633481 100644 --- a/src/canon/Physics.jl +++ b/src/canon/Physics.jl @@ -35,7 +35,7 @@ end) f∧v - ∘(⋆,d,⋆)(uˢ)∧v - d(p) + b∧g - ∘(⋆,d,⋆)(τ) + uˢ̇ + Fᵥ uˢ̇ == force(U) - + end ) @@ -50,7 +50,7 @@ end) # Fick's first law ϕ == k(d₀(C)) end -) +) @docapode("Advection" ,"https://en.wikipedia.org/wiki/Advection" @@ -115,14 +115,14 @@ end) ) @docapode("Navier-Stokes" - ,"https://en.wikipedia.org/wiki/Navier_Stokes_equation" + ,"https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations" ,"Partial differential equations which describe the motion of viscous fluid surfaces." ,navier_stokes ,begin (V, V̇, G)::Form1{X} (ρ, ṗ, p)::Form0{X} - V̇ == neg₁(L₁′(V, V)) + + V̇ == neg₁(L₁′(V, V)) + div₁(kᵥ(Δ₁(V) + third(d₀(δ₁(V)))), avg₀₁(ρ)) + d₀(half(i₁′(V, V))) + neg₁(div₁(d₀(p),avg₀₁(ρ))) + @@ -157,7 +157,7 @@ end) ) @docapode("Poiseuille Density" - ,"https://en.wikipedia.org/wiki/hagen-poiseuille_density" + ,"https://en.wikipedia.org/wiki/Hagen%E2%80%93Poiseuille_equation" ,"" ,poiseuille_density ,begin @@ -206,7 +206,7 @@ end) ) @docapode(Lie - ,"https://en.wikipedia.org/wiki/lie_derivative" + ,"https://en.wikipedia.org/wiki/Lie_derivative" ,"" ,lie ,begin @@ -219,7 +219,7 @@ end) ) @docapode("Superposition" - ,"https://en.wikipedia.org/wiki/superposition" + ,"https://en.wikipedia.org/wiki/Superposition" ,"" ,superposition ,begin