From ef97e960f845fbe54579a2387bea7ea0acd34abd Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Sun, 11 Feb 2024 19:48:34 -0500 Subject: [PATCH 01/11] Revisit NHS momentum eqn and add turbulence closure --- examples/climate/nonhydrostatic_buoyancy.jl | 76 +++++++++++++++------ 1 file changed, 55 insertions(+), 21 deletions(-) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index 91740df9..81d98085 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -1,4 +1,4 @@ -# This is a small implementation of Oceananigans.jl's NonhydrostaticModel. +# This is a discrete exterior calculus implementation of Oceananigans.jl's `NonhydrostaticModel`. ####################### # Import Dependencies # @@ -11,7 +11,7 @@ using Decapodes # External Dependencies using GeometryBasics: Point3 -using GLMakie +using CairoMakie using JLD2 using LinearAlgebra using OrdinaryDiffEq @@ -23,22 +23,38 @@ Point3D = Point3{Float64} #################### # Equation 1: "The momentum conservation equation" from https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-momentum-conservation-equation -momentum = @decapode begin - (f,b)::Form0 - (v,V,g,Fᵥ,uˢ,v_up)::Form1 - τ::Form2 - U::Parameter - - uˢ̇ == ∂ₜ(uˢ) - - v_up == -1 * L(v,v) - L(V,v) - L(v,V) - - f∧v - ∘(⋆,d,⋆)(uˢ)∧v - d(p) + b∧g - ∘(⋆,d,⋆)(τ) + uˢ̇ + Fᵥ - - uˢ̇ == force(U) +@decapode momentum begin + (v,V)::DualForm1 + f::PrimalForm0 + uˢ::DualForm1 + p::DualForm0 + b::Constant + ĝ::DualForm1 + Fᵥ::DualForm1 + + ∂ₜ(v) == + -ℒ(v,v) + (1/2)dι(v,v) - + dι(v,V) + ι(v,dV) + ι(V,dv) - + (f - ∘(d,⋆)(uˢ)) ∧ v - + d(p) + + b∧ĝ - + StressDivergence + + ∂ₜ(uˢ) + + Fᵥ end momentum = expand_operators(momentum) to_graphviz(momentum, graph_attrs=Dict(:rankdir => "LR")) +# Why write "StressDivergence" instead of ∇⋅τ? +# According to this docs page: +# https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/ +# , the user simply selects what model to insert in place of the term ∇⋅τ. +# For example, in the isotropic case, Oceananigans.jl rewrites 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. + # Equation 2: "The tracer conservation equation" from https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-tracer-conservation-equation tracer = @decapode begin (c,C,F,c_up)::Form0 @@ -57,6 +73,15 @@ equation_of_state = @decapode begin end to_graphviz(equation_of_state) +# Equation 2: "Constant isotropic diffusivity" from https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity +@decapode isotropic_diffusivity begin + StressDivergence::DualForm1 + nu::Constant + + StressDivergence == nu*Δ(v) +end +to_graphviz(equation_of_state) + boundary_conditions = @decapode begin (S,T)::Form0 (Ṡ,T_up)::Form0 @@ -72,7 +97,10 @@ end to_graphviz(boundary_conditions) buoyancy_composition_diagram = @relation () begin - momentum(V, v, v_up, b) + momentum(V, v, v_up, b, SD) + + # "The turbulence closure selected by the user determines the form of stress divergence" + turbulence(SD) # "Both T and S obey the tracer conservation equation" temperature(V, v, T, T_up) @@ -86,11 +114,12 @@ end to_graphviz(buoyancy_composition_diagram, box_labels=:name, junction_labels=:variable, prog="fdp", graph_attrs=Dict(["sep" => "1.5"])) buoyancy_cospan = oapply(buoyancy_composition_diagram, [ - Open(momentum, [:V, :v, :v_up, :b]), - Open(tracer, [:V, :v, :c, :c_up]), - Open(tracer, [:V, :v, :c, :c_up]), - Open(equation_of_state, [:b, :T, :S]), - Open(boundary_conditions, [:v, :S, :T, :v_up, :Ṡ, :T_up])]) + Open(momentum, [:V, :v, :v_up, :b, :StressDivergence]), + Open(isotropic_diffusivity, [:StressDivergence]), + Open(tracer, [:V, :v, :c, :c_up]), + Open(tracer, [:V, :v, :c, :c_up]), + Open(equation_of_state, [:b, :T, :S]), + Open(boundary_conditions, [:v, :S, :T, :v_up, :Ṡ, :T_up])]) buoyancy = apex(buoyancy_cospan) to_graphviz(buoyancy) @@ -100,7 +129,6 @@ infer_types!(buoyancy) resolve_overloads!(buoyancy) to_graphviz(buoyancy) - #s′ = loadmesh(Rectangle_30x10()) s′ = triangulated_grid(80,80, 10, 10, Point3D) s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s′) @@ -115,6 +143,11 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin :L₁ => (x,y) -> L₁′(x,y,sd,hodge) :L₂ᵈ => (x,y) -> lie_derivative_flat(2, sd, x, y) + :ℒ => (x,y) -> begin + # ℒ_y(x) := ⋆(⋆(d(x)∧y)) + # TODO: Check the order of arguments to ∧. + star2_mat * ∧(Tuple{1,1}, sd, (star1_mat \ (duald0_mat * x)), y) + end :force => x -> x :∂_spatial => x -> begin left = findall(x -> x[1] ≈ 0.0, point(sd)) @@ -271,3 +304,4 @@ triple_tracer = expand_operators(triple_tracer) infer_types!(triple_tracer) resolve_overloads!(triple_tracer) to_graphviz(triple_tracer) + From 7a721cc47c28133600925e00e413e11544ff3e1b Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Sun, 25 Feb 2024 09:11:02 -0500 Subject: [PATCH 02/11] Type and bind operators --- examples/climate/nonhydrostatic_buoyancy.jl | 207 +++++++++++--------- 1 file changed, 109 insertions(+), 98 deletions(-) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index 81d98085..56afa68a 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -7,15 +7,18 @@ # AlgebraicJulia Dependencies using Catlab using CombinatorialSpaces +using CombinatorialSpaces: interior_product_dd, ℒ_dd using Decapodes +using DiagrammaticEquations # External Dependencies -using GeometryBasics: Point3 using CairoMakie +using ComponentArrays +using GeometryBasics: Point3 using JLD2 using LinearAlgebra +using MLStyle using OrdinaryDiffEq -using ComponentArrays Point3D = Point3{Float64} #################### @@ -23,27 +26,26 @@ Point3D = Point3{Float64} #################### # Equation 1: "The momentum conservation equation" from https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-momentum-conservation-equation -@decapode momentum begin +momentum = @decapode begin (v,V)::DualForm1 - f::PrimalForm0 + f::Form0 uˢ::DualForm1 p::DualForm0 - b::Constant + b::DualForm0 ĝ::DualForm1 Fᵥ::DualForm1 + StressDivergence::DualForm1 ∂ₜ(v) == - -ℒ(v,v) + (1/2)dι(v,v) - - dι(v,V) + ι(v,dV) + ι(V,dv) - - (f - ∘(d,⋆)(uˢ)) ∧ v - + -ℒ(v,v) + 0.5*d(ι₁₁(v,v)) - + d(ι₁₁(v,V)) + ι₁₂(v,d(V)) + ι₁₂(V,d(v)) - + (f - ∘(d,⋆)(uˢ)) ∧ᵖ v - d(p) + - b∧ĝ - + b ∧ᵈ ĝ - StressDivergence + ∂ₜ(uˢ) + Fᵥ end -momentum = expand_operators(momentum) -to_graphviz(momentum, graph_attrs=Dict(:rankdir => "LR")) # Why write "StressDivergence" instead of ∇⋅τ? # According to this docs page: @@ -56,127 +58,136 @@ to_graphviz(momentum, graph_attrs=Dict(:rankdir => "LR")) # we can operate purely in terms of scalar-valued forms. # Equation 2: "The tracer conservation equation" from https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-tracer-conservation-equation -tracer = @decapode begin - (c,C,F,c_up)::Form0 - (v,V,q)::Form1 +tracer_conservation = @decapode begin + (c,C,F,c_up,FluxDivergence)::DualForm0 + (v,V)::DualForm1 - c_up == -1*⋆(L(v,⋆(c))) - ⋆(L(V,⋆(c))) - ⋆(L(v,⋆(C))) - ∘(⋆,d,⋆)(q) + F + ∂ₜ(c) == + -1*ι₁₁(v,d(c)) - + ι₁₁(V,d(c)) - + ι₁₁(v,d(C)) - + FluxDivergence + + F end -to_graphviz(tracer, graph_attrs=Dict(:rankdir => "LR")) # Equation 2: "Linear equation of state" of seawater buoyancy from https://clima.github.io/OceananigansDocumentation/stable/physics/buoyancy_and_equations_of_state/#Linear-equation-of-state equation_of_state = @decapode begin - (b,T,S)::Form0 + (b,T,S)::DualForm0 (g,α,β)::Constant b == g*(α*T - β*S) end -to_graphviz(equation_of_state) # Equation 2: "Constant isotropic diffusivity" from https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity -@decapode isotropic_diffusivity begin +isotropic_diffusivity = @decapode begin + v::DualForm1 + c::DualForm0 StressDivergence::DualForm1 - nu::Constant + FluxDivergence::DualForm0 + (κ,nu)::Constant StressDivergence == nu*Δ(v) + FluxDivergence == κ*Δ(c) end -to_graphviz(equation_of_state) - -boundary_conditions = @decapode begin - (S,T)::Form0 - (Ṡ,T_up)::Form0 - v::Form1 - v_up::Form1 - Ṫ == ∂ₜ(T) - Ṡ == ∂ₜ(S) - v̇ == ∂ₜ(v) - - Ṫ == ∂_spatial(T_up) - v̇ == ∂_noslip(v_up) + +################## +# Compose Models # +################## + +""" +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. + +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. +""" + +# Specify the equations that a tracer obeys: +tracer_composition = @relation () begin + # "The turbulence closure selected by the user determines the form of ... diffusive flux divergence" + turbulence(FD,v) + + continuity(FD,v) end -to_graphviz(boundary_conditions) -buoyancy_composition_diagram = @relation () begin - momentum(V, v, v_up, b, SD) +# Let's "lock in" isotropic diffusivity by doing an intermediate oapply. +isotropic_tracer = apex(oapply(tracer_composition, [ + Open(isotropic_diffusivity, [:FluxDivergence, :v]), + Open(tracer_conservation, [:FluxDivergence, :v])])) +# Use this building-block tracer physics at the next level: + +nonhydrostatic_composition = @relation () begin # "The turbulence closure selected by the user determines the form of stress divergence" - turbulence(SD) + # => 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(V, v, T, T_up) - salinity(V, v, S, S_up) + # => Temperature and Salinity both receive a copy of the tracer physics. + temperature(V, v, T, SD) + salinity(V, v, S, SD) # "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) - - bcs(v, S, T, v_up, S_up, T_up) end -to_graphviz(buoyancy_composition_diagram, box_labels=:name, junction_labels=:variable, prog="fdp", graph_attrs=Dict(["sep" => "1.5"])) - -buoyancy_cospan = oapply(buoyancy_composition_diagram, [ - Open(momentum, [:V, :v, :v_up, :b, :StressDivergence]), - Open(isotropic_diffusivity, [:StressDivergence]), - Open(tracer, [:V, :v, :c, :c_up]), - Open(tracer, [:V, :v, :c, :c_up]), - Open(equation_of_state, [:b, :T, :S]), - Open(boundary_conditions, [:v, :S, :T, :v_up, :Ṡ, :T_up])]) - -buoyancy = apex(buoyancy_cospan) -to_graphviz(buoyancy) - -buoyancy = expand_operators(buoyancy) -infer_types!(buoyancy) -resolve_overloads!(buoyancy) -to_graphviz(buoyancy) - -#s′ = loadmesh(Rectangle_30x10()) -s′ = triangulated_grid(80,80, 10, 10, Point3D) -s = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s′) -subdivide_duals!(s, Barycenter()) + +isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ + Open(momentum, [:V, :v, :b, :StressDivergence]), + Open(isotropic_tracer, [:continuity_V, :v, :continuity_c, :turbulence_StressDivergence]), + Open(isotropic_tracer, [:continuity_V, :v, :continuity_c, :turbulence_StressDivergence]), + Open(equation_of_state, [:b, :T, :S])])) + +################# +# Define a Mesh # +################# + +# This is the sphere with the same surface area as the Oceananigans +# example torus: +#s = loadmesh(Icosphere(5, (3*π)^(1/3))) + +# 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()) xmax = maximum(x -> x[1], point(s)) zmax = maximum(x -> x[2], point(s)) -wireframe(s) - -# Some extra operators that haven't been up-streamed into CombinatorialSpaces yet. -include("dec_operators.jl") +#wireframe(s) + +# TODO: Provide these functions by default. +i11 = interior_product_dd(Tuple{1,1}, sd); +i12 = interior_product_dd(Tuple{1,2}, sd); +lie11 = ℒ_dd(Tuple{1,1}, sd); +Λᵖ = dec_wedge_product_pd(Tuple{0,1}, sd); +Λᵈ = dec_wedge_product_dd(Tuple{0,1}, sd); +# TODO: Upstream the dual 0 Laplacian. +dd0 = dec_dual_derivative(0, sd); +ihs1 = dec_inv_hodge_star(1, sd, GeometricHodge()); +d1 = dec_differential(1,sd); +hs2 = dec_hodge_star(2, sd, GeometricHodge()); function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin - :L₁ => (x,y) -> L₁′(x,y,sd,hodge) - :L₂ᵈ => (x,y) -> lie_derivative_flat(2, sd, x, y) - :ℒ => (x,y) -> begin - # ℒ_y(x) := ⋆(⋆(d(x)∧y)) - # TODO: Check the order of arguments to ∧. - star2_mat * ∧(Tuple{1,1}, sd, (star1_mat \ (duald0_mat * x)), y) - end - :force => x -> x - :∂_spatial => x -> begin - left = findall(x -> x[1] ≈ 0.0, point(sd)) - right = findall(x -> x[1] ≈ xmax, point(sd)) - x[left] .= 0.0 - x[right] .= 0.0 - x - end - :∂_noslip => x -> begin - right = findall(x -> x[1] ≈ xmax, point(sd)) - top = findall(x -> x[2] ≈ zmax, point(sd)) - bottom = findall(x -> x[2] ≈ 0.0, point(sd)) - x[findall(x -> x[1] ≈ 0.0 , point(sd,sd[:∂v0]))] .= 0 - x[findall(x -> x[1] ≈ 0.0 , point(sd,sd[:∂v1]))] .= 0 - x[findall(x -> x[1] ≈ xmax, point(sd,sd[:∂v0]))] .= 0 - x[findall(x -> x[1] ≈ xmax, point(sd,sd[:∂v1]))] .= 0 - x[findall(x -> x[2] ≈ 0.0 , point(sd,sd[:∂v0]))] .= 0 - x[findall(x -> x[2] ≈ 0.0 , point(sd,sd[:∂v1]))] .= 0 - x[findall(x -> x[2] ≈ zmax, point(sd,sd[:∂v0]))] .= 0 - x[findall(x -> x[2] ≈ zmax, point(sd,sd[:∂v1]))] .= 0 - x - end + :ℒ => lie11 + :ι₁₁ => i11 + :ι₁₂ => i12 + :∧ᵖ => Λᵖ + :∧ᵈ => Λᵈ + :Δ => x -> hs2 * d1 * ihs1(dd0 * x) _ => default_dec_generate(sd, my_symbol, hodge) end return (args...) -> op(args...) end -sim = eval(gensim(buoyancy)) -fₘ = sim(s, generate) +sim = eval(gensim(isotropic_nonhydrostatic_buoyancy)) +fₘ = sim(sd, generate) S = map(point(s)) do (_,_,_) 35.0 From 0d4d91541daf58c28944bee199f3be912cf9aee7 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Sun, 25 Feb 2024 09:40:40 -0500 Subject: [PATCH 03/11] Compose on c --- examples/climate/nonhydrostatic_buoyancy.jl | 106 +++++++++++--------- 1 file changed, 59 insertions(+), 47 deletions(-) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index 56afa68a..02993cb1 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -59,7 +59,7 @@ end # Equation 2: "The tracer conservation equation" from https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-tracer-conservation-equation tracer_conservation = @decapode begin - (c,C,F,c_up,FluxDivergence)::DualForm0 + (c,C,F,FluxDivergence)::DualForm0 (v,V)::DualForm1 ∂ₜ(c) == @@ -110,15 +110,15 @@ We will use our operad algebra to guarantee model compatibility and physical con # Specify the equations that a tracer obeys: tracer_composition = @relation () begin # "The turbulence closure selected by the user determines the form of ... diffusive flux divergence" - turbulence(FD,v) + turbulence(FD,v,c) - continuity(FD,v) + continuity(FD,v,c) end # Let's "lock in" isotropic diffusivity by doing an intermediate oapply. isotropic_tracer = apex(oapply(tracer_composition, [ - Open(isotropic_diffusivity, [:FluxDivergence, :v]), - Open(tracer_conservation, [:FluxDivergence, :v])])) + Open(isotropic_diffusivity, [:FluxDivergence, :v, :c]), + Open(tracer_conservation, [:FluxDivergence, :v, :c])])) # Use this building-block tracer physics at the next level: @@ -139,8 +139,8 @@ end isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ Open(momentum, [:V, :v, :b, :StressDivergence]), - Open(isotropic_tracer, [:continuity_V, :v, :continuity_c, :turbulence_StressDivergence]), - Open(isotropic_tracer, [:continuity_V, :v, :continuity_c, :turbulence_StressDivergence]), + Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_StressDivergence]), + Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_StressDivergence]), Open(equation_of_state, [:b, :T, :S])])) ################# @@ -162,6 +162,10 @@ xmax = maximum(x -> x[1], point(s)) zmax = maximum(x -> x[2], point(s)) #wireframe(s) +################################# +# Define Differential Operators # +################################# + # TODO: Provide these functions by default. i11 = interior_product_dd(Tuple{1,1}, sd); i12 = interior_product_dd(Tuple{1,2}, sd); @@ -186,59 +190,67 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) return (args...) -> op(args...) end +####################### +# Generate Simulation # +####################### + sim = eval(gensim(isotropic_nonhydrostatic_buoyancy)) fₘ = sim(sd, generate) -S = map(point(s)) do (_,_,_) - 35.0 +################################# +# Define Differential Operators # +################################# + +S = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) + 0.0 end -T = map(point(s)) do (x,z,_) - #273.15 + 4 + ((zmax-z)^2 + (xmax-x)^2)^(1/2)/(1e2) - 273.15 + 4 +T = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) + 0.0 end -left = findall(x -> x[1] ≈ 0.0, point(sd)) -right = findall(x -> x[1] ≈ xmax, point(sd)) -T[left] .= 276 -T[right] .= 279 -extrema(T) -mesh(s′, color=T, colormap=:jet) -p = map(point(s)) do (x,z,_) - (zmax-z) +p = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) + 0.0 end -extrema(p) -f = zeros(nv(s)) -Fₛ = zeros(nv(s)) -Fₜ = zeros(nv(s)) -f = zeros(nv(s)) -Cₛ = zeros(nv(s)) -Cₜ = zeros(nv(s)) - -V = zeros(ne(s)) -v = zeros(ne(s)) -g = ♭(sd, DualVectorField(fill(Point3D(0,1,0), ntriangles(s)))).data -Fᵥ = zeros(ne(s)) -qₛ = zeros(ne(s)) -qₜ = zeros(ne(s)) -uˢ = zeros(ne(s)) -dtuˢ = zeros(ne(s)) - -τ = zeros(ntriangles(s)) - -u₀ = ComponentArrays(momentum_f=f,v=v,V=V,momentum_g=g, - momentum_Fᵥ=Fᵥ,momentum_uˢ=uˢ,momentum_τ=τ, - momentum_dtuˢ=dtuˢ,momentum_p=p,T=T, - temperature_F=Fₜ,temperature_q=qₜ, - temperature_C=Cₜ,S=S,salinity_F=Fₛ, - salinity_q=qₛ,salinity_C=Cₛ) +f = zeros(nv(sd)) +Fₛ = zeros(ntriangles(sd)) +Fₜ = zeros(ntriangles(sd)) +f = zeros(ntriangles(sd)) +Cₛ = zeros(ntriangles(sd)) +Cₜ = zeros(ntriangles(sd)) +V = zeros(ne(sd)) +v = zeros(ne(sd)) +ĝ = ♭(sd, DualVectorField(fill(Point3D(0,1,0), ntriangles(sd)))).data +Fᵥ = zeros(ne(sd)) +qₛ = zeros(ne(sd)) +qₜ = zeros(ne(sd)) +uˢ = zeros(ne(sd)) + +u₀ = ComponentArray( + v = v, + V = V, + momentum_f = f, + momentum_uˢ = uˢ, + 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, + temperature_turbulence_nu = 0.0, + salinity_turbulence_κ = 0.0, + salinity_turbulence_nu = 0.0, eos_g = gᶜ, eos_α = α, - eos_β = β, - momentum_U = t -> 0) + eos_β = β) + tₑ = 1.5 From 3edb89260fcaba921aa6ea42b98fd2141bb0adc6 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Mon, 26 Feb 2024 23:57:00 -0500 Subject: [PATCH 04/11] Define dual vector laplacian --- examples/climate/nonhydrostatic_buoyancy.jl | 27 ++++++++++++++++----- 1 file changed, 21 insertions(+), 6 deletions(-) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index 02993cb1..f44bd8b8 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -30,6 +30,7 @@ momentum = @decapode begin (v,V)::DualForm1 f::Form0 uˢ::DualForm1 + ∂tuˢ::DualForm1 p::DualForm0 b::DualForm0 ĝ::DualForm1 @@ -43,7 +44,7 @@ momentum = @decapode begin d(p) + b ∧ᵈ ĝ - StressDivergence + - ∂ₜ(uˢ) + + ∂tuˢ + Fᵥ end @@ -86,8 +87,8 @@ isotropic_diffusivity = @decapode begin FluxDivergence::DualForm0 (κ,nu)::Constant - StressDivergence == nu*Δ(v) - FluxDivergence == κ*Δ(c) + StressDivergence == nu*Δᵈ₁(v) + FluxDivergence == κ*Δᵈ₀(c) end ################## @@ -177,6 +178,11 @@ dd0 = dec_dual_derivative(0, sd); ihs1 = dec_inv_hodge_star(1, sd, GeometricHodge()); d1 = dec_differential(1,sd); hs2 = dec_hodge_star(2, sd, GeometricHodge()); +# TODO: Upstream the dual 1 Laplacian. +dd1 = dec_dual_derivative(1, sd); +ihs0 = dec_inv_hodge_star(0, sd, GeometricHodge()); +d0 = dec_differential(0,sd); +hs1 = dec_hodge_star(1, sd, GeometricHodge()); function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin :ℒ => lie11 @@ -184,7 +190,11 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) :ι₁₂ => i12 :∧ᵖ => Λᵖ :∧ᵈ => Λᵈ - :Δ => x -> hs2 * d1 * ihs1(dd0 * x) + :Δᵈ₀ => x -> hs2 * d1 * ihs1(dd0 * x) + :Δᵈ₁ => x -> begin + hs1 * d0 * ihs0 * dd1 * x + + dd0 * hs2 * d1 * ihs1(x) + end _ => default_dec_generate(sd, my_symbol, hodge) end return (args...) -> op(args...) @@ -194,7 +204,11 @@ end # Generate Simulation # ####################### -sim = eval(gensim(isotropic_nonhydrostatic_buoyancy)) +open("nhs.jl", "w") do f + write(f, string(gensim(expand_operators(isotropic_nonhydrostatic_buoyancy)))) +end +sim = include("nhs.jl") +#sim = eval(gensim(isotropic_nonhydrostatic_buoyancy)) fₘ = sim(sd, generate) ################################# @@ -213,7 +227,6 @@ end f = zeros(nv(sd)) Fₛ = zeros(ntriangles(sd)) Fₜ = zeros(ntriangles(sd)) -f = zeros(ntriangles(sd)) Cₛ = zeros(ntriangles(sd)) Cₜ = zeros(ntriangles(sd)) V = zeros(ne(sd)) @@ -223,12 +236,14 @@ 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ᵥ, From 139478d42563f16455f7a9976608978c640cab49 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 27 Feb 2024 00:35:21 -0500 Subject: [PATCH 05/11] Execute on random v --- examples/climate/nonhydrostatic_buoyancy.jl | 23 ++++++++++++++------- 1 file changed, 16 insertions(+), 7 deletions(-) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index f44bd8b8..db35e953 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -17,8 +17,11 @@ using ComponentArrays using GeometryBasics: Point3 using JLD2 using LinearAlgebra +using Logging: global_logger using MLStyle using OrdinaryDiffEq +using TerminalLoggers: TerminalLogger +global_logger(TerminalLogger()) Point3D = Point3{Float64} #################### @@ -230,8 +233,11 @@ Fₜ = zeros(ntriangles(sd)) Cₛ = zeros(ntriangles(sd)) Cₜ = zeros(ntriangles(sd)) V = zeros(ne(sd)) -v = zeros(ne(sd)) -ĝ = ♭(sd, DualVectorField(fill(Point3D(0,1,0), ntriangles(sd)))).data +#v = zeros(ne(sd)) +#v = rand(ne(sd)) *1e-5 +v = rand(ne(sd)) *1e-8 +#ĝ = ♭(sd, DualVectorField(fill(Point3D(0,1,0), ntriangles(sd)))).data +ĝ = ♭(sd, DualVectorField(fill(Point3D(0,0,0), ntriangles(sd)))).data Fᵥ = zeros(ne(sd)) qₛ = zeros(ne(sd)) qₜ = zeros(ne(sd)) @@ -267,7 +273,7 @@ constants_and_parameters = ( eos_β = β) -tₑ = 1.5 +tₑ = 5e4 # Julia will pre-compile the generated simulation the first time it is run. @info("Precompiling Solver") @@ -277,13 +283,16 @@ soln.retcode != :Unstable || error("Solver was not stable") @info("Solving") prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) -soln = solve(prob, Tsit5(), dtmin=1e-3, force_dtmin=true) +soln = solve(prob, Tsit5(), dtmin=1e-16, force_dtmin=true, progress=true) @show soln.retcode @info("Done") -mesh(s′, color=soln(1.5).T, colormap=:jet) -extrema(soln(0).T) -extrema(soln(1.5).T) +# Vorticity: +mesh(s, + color=ihs0*dd1*soln(tₑ).v, + colormap=:jet) + +# Speed: # Create a gif begin From 66c29d33c1138326455554567eed7c9db120e28e Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 27 Feb 2024 09:48:52 -0500 Subject: [PATCH 06/11] Share nu between turbulence closures --- examples/climate/nonhydrostatic_buoyancy.jl | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index db35e953..bb4d63e3 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -133,8 +133,8 @@ nonhydrostatic_composition = @relation () begin # "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) - salinity(V, v, S, SD) + 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. @@ -143,8 +143,8 @@ end isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ Open(momentum, [:V, :v, :b, :StressDivergence]), - Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_StressDivergence]), - Open(isotropic_tracer, [:continuity_V, :v, :c, :turbulence_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])])) ################# @@ -211,7 +211,6 @@ open("nhs.jl", "w") do f write(f, string(gensim(expand_operators(isotropic_nonhydrostatic_buoyancy)))) end sim = include("nhs.jl") -#sim = eval(gensim(isotropic_nonhydrostatic_buoyancy)) fₘ = sim(sd, generate) ################################# @@ -265,9 +264,8 @@ gᶜ = 9.81 β = 5e-4 constants_and_parameters = ( temperature_turbulence_κ = 0.0, - temperature_turbulence_nu = 0.0, salinity_turbulence_κ = 0.0, - salinity_turbulence_nu = 0.0, + nu = 0.0, eos_g = gᶜ, eos_α = α, eos_β = β) From 40a242974269dbda6e4da4cba6ba90d75bdd19c0 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 27 Feb 2024 12:38:43 -0500 Subject: [PATCH 07/11] Upstream lie, interior product, sharp --- docs/src/cism.md | 8 +--- examples/climate/nonhydrostatic_buoyancy.jl | 41 +++------------------ src/operators.jl | 36 ++++++++++++++++++ 3 files changed, 43 insertions(+), 42 deletions(-) diff --git a/docs/src/cism.md b/docs/src/cism.md index e2402a82..14f7b2db 100644 --- a/docs/src/cism.md +++ b/docs/src/cism.md @@ -194,22 +194,16 @@ We provide here the mapping from symbols to differential operators. As more of t # Define how symbols map to Julia functions # ############################################# -# This sharp operator, ♯, is scheduled to be upstreamed. -include("sharp_op.jl") function generate(sd, my_symbol; hodge=GeometricHodge()) # We pre-allocate matrices that encode differential operators. - ♯_m = ♯_mat(sd) op = @match my_symbol begin - :♯ => x -> begin - ♯_m * x - end :mag => x -> begin norm.(x) end :^ => (x,y) -> begin x .^ y end - x => error("Unmatched operator $my_symbol") + _ => default_dec_matrix_generate(sd, my_symbol, hodge) end return (args...) -> op(args...) end diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index bb4d63e3..892be5b3 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -7,7 +7,6 @@ # AlgebraicJulia Dependencies using Catlab using CombinatorialSpaces -using CombinatorialSpaces: interior_product_dd, ℒ_dd using Decapodes using DiagrammaticEquations @@ -41,11 +40,11 @@ momentum = @decapode begin StressDivergence::DualForm1 ∂ₜ(v) == - -ℒ(v,v) + 0.5*d(ι₁₁(v,v)) - + -ℒ₁(v,v) + 0.5*d(ι₁₁(v,v)) - d(ι₁₁(v,V)) + ι₁₂(v,d(V)) + ι₁₂(V,d(v)) - - (f - ∘(d,⋆)(uˢ)) ∧ᵖ v - + (f - ∘(d,⋆)(uˢ)) ∧ᵖᵈ₀₁ v - d(p) + - b ∧ᵈ ĝ - + b ∧ᵈᵈ₀₁ ĝ - StressDivergence + ∂tuˢ + Fᵥ @@ -162,43 +161,14 @@ isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ s = EmbeddedDeltaSet2D("torus.obj") sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s) subdivide_duals!(sd, Barycenter()) -xmax = maximum(x -> x[1], point(s)) -zmax = maximum(x -> x[2], point(s)) -#wireframe(s) ################################# # Define Differential Operators # ################################# -# TODO: Provide these functions by default. -i11 = interior_product_dd(Tuple{1,1}, sd); -i12 = interior_product_dd(Tuple{1,2}, sd); -lie11 = ℒ_dd(Tuple{1,1}, sd); -Λᵖ = dec_wedge_product_pd(Tuple{0,1}, sd); -Λᵈ = dec_wedge_product_dd(Tuple{0,1}, sd); -# TODO: Upstream the dual 0 Laplacian. -dd0 = dec_dual_derivative(0, sd); -ihs1 = dec_inv_hodge_star(1, sd, GeometricHodge()); -d1 = dec_differential(1,sd); -hs2 = dec_hodge_star(2, sd, GeometricHodge()); -# TODO: Upstream the dual 1 Laplacian. -dd1 = dec_dual_derivative(1, sd); -ihs0 = dec_inv_hodge_star(0, sd, GeometricHodge()); -d0 = dec_differential(0,sd); -hs1 = dec_hodge_star(1, sd, GeometricHodge()); function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin - :ℒ => lie11 - :ι₁₁ => i11 - :ι₁₂ => i12 - :∧ᵖ => Λᵖ - :∧ᵈ => Λᵈ - :Δᵈ₀ => x -> hs2 * d1 * ihs1(dd0 * x) - :Δᵈ₁ => x -> begin - hs1 * d0 * ihs0 * dd1 * x + - dd0 * hs2 * d1 * ihs1(x) - end - _ => default_dec_generate(sd, my_symbol, hodge) + _ => default_dec_matrix_generate(sd, my_symbol, hodge) end return (args...) -> op(args...) end @@ -281,7 +251,8 @@ soln.retcode != :Unstable || error("Solver was not stable") @info("Solving") prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) -soln = solve(prob, Tsit5(), dtmin=1e-16, force_dtmin=true, progress=true) +#soln = solve(prob, Tsit5(), dtmin=1e-16, force_dtmin=true, progress=true) +soln = solve(prob, Tsit5(), force_dtmin=true, progress=true) @show soln.retcode @info("Done") diff --git a/src/operators.jl b/src/operators.jl index 6228ac94..831b1ce5 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -31,6 +31,32 @@ function default_dec_matrix_generate(sd, my_symbol, hodge) :∧₂₀ => dec_pair_wedge_product(Tuple{2,0}, sd) :∧₁₁ => dec_pair_wedge_product(Tuple{1,1}, sd) + # Primal-Dual Wedge Products + :∧ᵖᵈ₁₁ => dec_wedge_product_pd(Tuple{1,1}, sd) + :∧ᵖᵈ₀₁ => dec_wedge_product_pd(Tuple{0,1}, sd) + :∧ᵈᵖ₁₁ => dec_wedge_product_dp(Tuple{1,1}, sd) + :∧ᵈᵖ₁₀ => dec_wedge_product_dp(Tuple{1,0}, sd) + + # Dual-Dual Wedge Products + :∧ᵈᵈ₁₁ => dec_wedge_product_dd(Tuple{1,1}, sd) + :∧ᵈᵈ₁₀ => dec_wedge_product_dd(Tuple{1,0}, sd) + :∧ᵈᵈ₀₁ => dec_wedge_product_dd(Tuple{0,1}, sd) + + # Dual-Dual Interior Products + :ι₁₁ => interior_product_dd(Tuple{1,1}, sd) + :ι₁₂ => interior_product_dd(Tuple{1,2}, sd) + + # Dual-Dual Lie Derivatives + :ℒ₁ => ℒ_dd(Tuple{1,1}, sd) + + # Dual Laplacians + :Δᵈ₀ => Δᵈ(Val{0},sd) + :Δᵈ₁ => Δᵈ(Val{1},sd) + + # Dual Laplacians + :♯ => dec_♯_p(sd) + :♯ᵈ => dec_♯_d(sd) + _ => error("Unmatched operator $my_symbol") end @@ -86,6 +112,16 @@ function dec_pair_wedge_product(::Type{Tuple{1,1}}, sd::HasDeltaSet2D) (α, β) -> dec_c_wedge_product(Tuple{1,1}, α, β, val_pack)) end +function dec_♯_p(sd::HasDeltaSet2D) + ♯_m = ♯_mat(sd, AltPPSharp()) + x -> ♯_m * x +end + +function dec_♯_d(sd::HasDeltaSet2D) + ♯_m = ♯_mat(sd, LLSDDSharp()) + x -> ♯_m * x +end + function default_dec_generate(sd, my_symbol, hodge=GeometricHodge()) op = @match my_symbol begin From cf6dc45cf737b442bbd67239cafbeb4cea806a5f Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 27 Feb 2024 16:39:48 -0500 Subject: [PATCH 08/11] Plot vorticity and speed in NHS --- examples/climate/nonhydrostatic_buoyancy.jl | 79 +++++++++++++-------- 1 file changed, 51 insertions(+), 28 deletions(-) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/examples/climate/nonhydrostatic_buoyancy.jl index 892be5b3..29c55595 100644 --- a/examples/climate/nonhydrostatic_buoyancy.jl +++ b/examples/climate/nonhydrostatic_buoyancy.jl @@ -157,7 +157,7 @@ isotropic_nonhydrostatic_buoyancy = apex(oapply(nonhydrostatic_composition, [ # 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") +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()) @@ -202,10 +202,7 @@ Fₜ = zeros(ntriangles(sd)) Cₛ = zeros(ntriangles(sd)) Cₜ = zeros(ntriangles(sd)) V = zeros(ne(sd)) -#v = zeros(ne(sd)) -#v = rand(ne(sd)) *1e-5 -v = rand(ne(sd)) *1e-8 -#ĝ = ♭(sd, DualVectorField(fill(Point3D(0,1,0), ntriangles(sd)))).data +v = rand(ne(sd)) * 1e-8 ĝ = ♭(sd, DualVectorField(fill(Point3D(0,0,0), ntriangles(sd)))).data Fᵥ = zeros(ne(sd)) qₛ = zeros(ne(sd)) @@ -235,13 +232,12 @@ gᶜ = 9.81 constants_and_parameters = ( temperature_turbulence_κ = 0.0, salinity_turbulence_κ = 0.0, - nu = 0.0, + nu = 1e-5, eos_g = gᶜ, eos_α = α, eos_β = β) - -tₑ = 5e4 +tₑ = 50 # Julia will pre-compile the generated simulation the first time it is run. @info("Precompiling Solver") @@ -251,31 +247,58 @@ soln.retcode != :Unstable || error("Solver was not stable") @info("Solving") prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) -#soln = solve(prob, Tsit5(), dtmin=1e-16, force_dtmin=true, progress=true) -soln = solve(prob, Tsit5(), force_dtmin=true, progress=true) +soln = solve(prob, Vern7(), force_dtmin=true, dtmax=0.2, progress=true,progress_steps=1) @show soln.retcode @info("Done") -# Vorticity: -mesh(s, - color=ihs0*dd1*soln(tₑ).v, - colormap=:jet) - -# Speed: - -# Create a gif -begin - frames = 100 - #fig, ax, ob = GLMakie.mesh(s′, color=soln(0).T, colormap=:jet, colorrange=extrema(soln(tₑ).h)) - fig, ax, ob = GLMakie.mesh(s′, color=soln(0).T, colormap=:jet, colorrange=extrema(soln(1.5).T)) - Colorbar(fig[1,2], ob) - #record(fig, "oceananigans.gif", range(0.0, tₑ; length=frames); framerate = 30) do t - record(fig, "oceananigans.gif", range(0.0, 1.5; length=frames); framerate = 30) do t - ob.color = soln(t).T - end +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 -begin 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(true) + +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(true) # Track a single tracer. single_tracer_composition_diagram = @relation () begin From b83cfeba891773fbc8e77dd4caae8bec88f32179 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 27 Feb 2024 16:41:18 -0500 Subject: [PATCH 09/11] Turn NHS into docs page --- examples/climate/nonhydrostatic_buoyancy.jl => docs/src/nhs.md | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename examples/climate/nonhydrostatic_buoyancy.jl => docs/src/nhs.md (100%) diff --git a/examples/climate/nonhydrostatic_buoyancy.jl b/docs/src/nhs.md similarity index 100% rename from examples/climate/nonhydrostatic_buoyancy.jl rename to docs/src/nhs.md From 9909b61306a245563f9328b907cd3cd94b89041a Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 27 Feb 2024 17:13:15 -0500 Subject: [PATCH 10/11] Use Documenter format for NHS page --- docs/src/nhs.md | 146 +++++++++++++++++------------------------------- 1 file changed, 51 insertions(+), 95 deletions(-) diff --git a/docs/src/nhs.md b/docs/src/nhs.md index 29c55595..0f39b8ed 100644 --- a/docs/src/nhs.md +++ b/docs/src/nhs.md @@ -1,9 +1,8 @@ -# This is a discrete exterior calculus implementation of Oceananigans.jl's `NonhydrostaticModel`. +# Implement Oceananigans.jl's NonhydrostaticModel in the Discrete Exterior Calculus -####################### -# Import Dependencies # -####################### +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. +```@example DEC # AlgebraicJulia Dependencies using Catlab using CombinatorialSpaces @@ -16,18 +15,16 @@ using ComponentArrays using GeometryBasics: Point3 using JLD2 using LinearAlgebra -using Logging: global_logger using MLStyle using OrdinaryDiffEq -using TerminalLoggers: TerminalLogger -global_logger(TerminalLogger()) -Point3D = Point3{Float64} +Point3D = Point3{Float64}; +``` -#################### -# Define the model # -#################### +## Specify our models. -# Equation 1: "The momentum conservation equation" from https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-momentum-conservation-equation +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 @@ -49,18 +46,12 @@ momentum = @decapode begin ∂tuˢ + Fᵥ end +``` + +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 write "StressDivergence" instead of ∇⋅τ? -# According to this docs page: -# https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/ -# , the user simply selects what model to insert in place of the term ∇⋅τ. -# For example, in the isotropic case, Oceananigans.jl rewrites 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. - -# Equation 2: "The tracer conservation equation" from https://clima.github.io/OceananigansDocumentation/stable/physics/nonhydrostatic_model/#The-tracer-conservation-equation +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 @@ -72,16 +63,20 @@ tracer_conservation = @decapode begin FluxDivergence + F end +``` -# Equation 2: "Linear equation of state" of seawater buoyancy from https://clima.github.io/OceananigansDocumentation/stable/physics/buoyancy_and_equations_of_state/#Linear-equation-of-state +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 +``` -# Equation 2: "Constant isotropic diffusivity" from https://clima.github.io/OceananigansDocumentation/stable/physics/turbulence_closures/#Constant-isotropic-diffusivity +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 @@ -92,12 +87,10 @@ isotropic_diffusivity = @decapode begin StressDivergence == nu*Δᵈ₁(v) FluxDivergence == κ*Δᵈ₀(c) end +``` -################## -# Compose Models # -################## +## 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. Besides it being convenient to break apart large UWDs into component UWDs, this hierarchical composition can enforce rules on our physical quantities. @@ -108,23 +101,27 @@ For example: 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. -""" -# Specify the equations that a tracer obeys: +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 +``` -# 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]), Open(tracer_conservation, [:FluxDivergence, :v, :c])])) +``` -# Use this building-block tracer physics at the next level: +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. @@ -144,16 +141,14 @@ 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])])) + Open(equation_of_state, [:b, :T, :S])])); +``` -################# -# Define a Mesh # -################# +## Meshes and Initial Conditions -# This is the sphere with the same surface area as the Oceananigans -# example torus: -#s = loadmesh(Icosphere(5, (3*π)^(1/3))) +Let's 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!) @@ -162,10 +157,6 @@ s = EmbeddedDeltaSet2D("torus.obj") sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point3D}(s) subdivide_duals!(sd, Barycenter()) -################################# -# Define Differential Operators # -################################# - function generate(sd, my_symbol; hodge=GeometricHodge()) op = @match my_symbol begin _ => default_dec_matrix_generate(sd, my_symbol, hodge) @@ -173,20 +164,12 @@ function generate(sd, my_symbol; hodge=GeometricHodge()) return (args...) -> op(args...) end -####################### -# Generate Simulation # -####################### - open("nhs.jl", "w") do f write(f, string(gensim(expand_operators(isotropic_nonhydrostatic_buoyancy)))) end sim = include("nhs.jl") fₘ = sim(sd, generate) -################################# -# Define Differential Operators # -################################# - S = map(sd[sd[:tri_center], :dual_point]) do (_,_,_) 0.0 end @@ -235,8 +218,14 @@ constants_and_parameters = ( nu = 1e-5, eos_g = gᶜ, eos_α = α, - eos_β = β) + eos_β = β); +``` +## Execute the Simulation + +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. @@ -250,7 +239,13 @@ prob = ODEProblem(fₘ, u₀, (0, tₑ), constants_and_parameters) soln = solve(prob, Vern7(), force_dtmin=true, dtmax=0.2, progress=true,progress_steps=1) @show soln.retcode @info("Done") +``` + +## 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}, sd, GeometricHodge()) dd1 = dec_dual_derivative(1, sd) ♯_m = ♯_mat(sd, LLSDDSharp()) @@ -299,48 +294,9 @@ function save_speed(is_2d=false) frames = 200 end end save_speed(true) +``` -# Track a single tracer. -single_tracer_composition_diagram = @relation () begin - momentum(V, v) - tracer(V, v) -end -to_graphviz(single_tracer_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") - -single_tracer_cospan = oapply(single_tracer_composition_diagram, - [Open(momentum, [:V, :v]), - Open(tracer, [:V, :v])]) - -single_tracer = apex(single_tracer_cospan) -to_graphviz(single_tracer) - -single_tracer = expand_operators(single_tracer) -infer_types!(single_tracer) -resolve_overloads!(single_tracer) -to_graphviz(single_tracer) - - -# Track multiple tracers. -triple_tracer_composition_diagram = @relation () begin - momentum(V, v) - - mercury(V, v) - phosphate(V, v) - oil(V, v) -end -to_graphviz(triple_tracer_composition_diagram, box_labels=:name, junction_labels=:variable, prog="circo") - -triple_tracer_cospan = oapply(triple_tracer_composition_diagram, - [Open(momentum, [:V, :v]), - Open(tracer, [:V, :v]), - Open(tracer, [:V, :v]), - Open(tracer, [:V, :v])]) - -triple_tracer = apex(triple_tracer_cospan) -to_graphviz(triple_tracer) +![Vorticity](vorticity.gif) -triple_tracer = expand_operators(triple_tracer) -infer_types!(triple_tracer) -resolve_overloads!(triple_tracer) -to_graphviz(triple_tracer) +![Speed](speed.gif) From be97b11731c69debd3bd2409888c2bb71e089659 Mon Sep 17 00:00:00 2001 From: Luke Morris Date: Tue, 27 Feb 2024 17:44:46 -0500 Subject: [PATCH 11/11] Plot torus in 3D and precompile solver --- docs/src/nhs.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/nhs.md b/docs/src/nhs.md index 0f39b8ed..30ff8cc6 100644 --- a/docs/src/nhs.md +++ b/docs/src/nhs.md @@ -231,7 +231,7 @@ 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, Tsit5(), dtmin=1e-3, force_dtmin=true) +soln = solve(prob, Vern7()) soln.retcode != :Unstable || error("Solver was not stable") @info("Solving") @@ -274,7 +274,7 @@ function save_vorticity(is_2d=false) time[] = t end end -save_vorticity(true) +save_vorticity(false) function save_speed(is_2d=false) frames = 200 time = Observable(0.0) @@ -293,7 +293,7 @@ function save_speed(is_2d=false) frames = 200 time[] = t end end -save_speed(true) +save_speed(false) ``` ![Vorticity](vorticity.gif)