Refactors for Sharp/Flat (#269)

Project.toml

@@ -25,6 +25,7 @@ DecapodesCUDAExt = "CUDA"
 ACSets = "0.2"
 Aqua = "0.8"
 CUDA = "5.2"
+Catlab = "0.16.19"
 CombinatorialSpaces = "0.7"
 ComponentArrays = "0.15"
 DiagrammaticEquations = "0.1.7"
@@ -45,6 +46,7 @@ julia = "1.9"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
@@ -55,4 +57,4 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "CUDA", "Distributions", "GeometryBasics", "JSON", "OrdinaryDiffEq", "Random", "Statistics", "Test"]
+test = ["Aqua", "Catlab", "CUDA", "Distributions", "GeometryBasics", "JSON", "OrdinaryDiffEq", "Random", "Statistics", "Test"]

docs/src/bsh/budyko_sellers_halfar.md

@@ -32,16 +32,16 @@ We have defined the [Halfar ice model](../ice_dynamics/ice_dynamics.md) in other
 ``` @example DEC
 halfar_eq2 = @decapode begin
   h::Form0
-  Γ::Form1
+  Γ::Form0
   n::Constant
 
   ḣ == ∂ₜ(h)
-  ḣ == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n-1)) * avg₀₁(h^(n+2)))
+  ḣ == Γ * ∘(⋆, d, ⋆)(d(h) * avg₀₁(mag(♯(d(h)))^(n-1)) * avg₀₁(h^(n+2)))
 end
 
 glens_law = @decapode begin
-  Γ::Form1
-  A::Form1
+  Γ::Form0
+  A::Form0
   (ρ,g,n)::Constant
   
   Γ == (2/(n+2))*A*(ρ*g)^n
@@ -147,9 +147,9 @@ We need to specify physically what it means for these two terms to interact. We
 ``` @example DEC
 warming = @decapode begin
   Tₛ::Form0
-  A::Form1
-
-  A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e7)
+  A::Form0
+  
+  A == 5.8282*10^(-0.236 * Tₛ)*1.65e7
 end
 
 to_graphviz(warming)

docs/src/cism/cism.md

@@ -53,11 +53,11 @@ Halfar's equation looks a little disjoint. It seems that the front most terms ar
 # translated into the exterior calculus.
 halfar_eq2 = @decapode begin
   h::Form0
-  Γ::Form1
+  Γ::Form0
   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
 
 to_graphviz(halfar_eq2)
@@ -72,7 +72,7 @@ Here, we recognize that Gamma is in fact what glaciologists call "Glen's Flow La
 # assumptions made in glacier theory, their experimental foundations and
 # consequences. (1958)
 glens_law = @decapode begin
-  Γ::Form1
+  Γ::Form0
   (A,ρ,g,n)::Constant
   
   Γ == (2/(n+2))*A*(ρ*g)^n
@@ -154,7 +154,7 @@ g = 9.8101
 alpha = 1/9
 beta = 1/18
 flwa = 1e-16
-A = fill(1e-16, ne(sd))
+A = 1e-16
 
 Gamma = 2.0/(n+2) * flwa * (ρ * g)^n
 t0 = (beta/Gamma) * (7.0/4.0)^3 * (R₀^4 / H^7)

docs/src/grigoriev/grigoriev.md

@@ -94,7 +94,7 @@ The coordinates of a vertex are stored in `sd[:point]`. Let's use our interpolat
 n = 3
 ρ = 910
 g = 9.8101
-A = fill(1e-16, ne(sd))
+A = 1e-16
 
 h₀ = map(sd[:point]) do (x,y,_)
   tif_val = ice_interp(x,y)
@@ -118,15 +118,15 @@ For exposition on this Halfar Decapode, see our [Glacial Flow](../ice_dynamics/i
 ``` @example DEC
 halfar_eq2 = @decapode begin
   h::Form0
-  Γ::Form1
+  Γ::Form0
   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
-  Γ::Form1
+  Γ::Form0
   (A,ρ,g,n)::Constant
   
   Γ == (2/(n+2))*A*(ρ*g)^n
@@ -151,10 +151,6 @@ to_graphviz(ice_dynamics)
 function generate(sd, my_symbol; hodge=GeometricHodge())
   op = @match my_symbol begin
     :mag => x -> norm.(x)
-    :♯ => begin
-      sharp_mat = ♯_mat(sd, AltPPSharp())
-      x -> sharp_mat * x
-    end
     x => error("Unmatched operator $my_symbol")
   end
   return op

docs/src/halmo/halmo.md

@@ -52,14 +52,14 @@ Halfar's equation and Glen's law are composed like so:
 ```@example DEC_halmo
 halfar_eq2 = @decapode begin
   h::Form0
-  Γ::Form1
+  Γ::Form0
   n::Constant
 
-  ∂ₜ(h) == ∘(⋆, d, ⋆)(Γ * d(h) * avg₀₁(mag(♯(d(h)))^(n-1)) * avg₀₁(h^(n+2)))
+  ∂ₜ(h) == Γ * ∘(⋆, d, ⋆)(d(h) ∧₁₀ ((mag(♯ᵖᵖ(d(h)))^(n-1)) ∧₀₀ h^(n+2)))
 end
 
 glens_law = @decapode begin
-  Γ::Form1
+  Γ::Form0
   (A,ρ,g,n)::Constant
   
   Γ == (2/(n+2))*A*(ρ*g)^n
@@ -147,7 +147,7 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
     :σ => sigmoid
     :mag => x -> norm.(x)
     # Remaining operations (such as our differential operators) are built-in.
-    _ => default_dec_matrix_generate(sd, my_symbol, hodge)
+    _ => error("Unmatched operator $my_symbol")
   end
   return op
 end;
@@ -181,7 +181,7 @@ u₀ = ComponentArray(
 
 constants_and_parameters = (
   glacier_dynamics_n = 3,
-  glacier_dynamics_stress_A = fill(1e-16, ne(sd)),
+  glacier_dynamics_stress_A = 1e-16,
   glacier_dynamics_stress_ρ = 910,
   glacier_dynamics_stress_g = 9.8101,
   water_dynamics_μ = 0.01);

docs/src/ice_dynamics/ice_dynamics.md

@@ -43,19 +43,17 @@ We'll change the term out front to Γ so we can demonstrate composition in a mom
 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})).
+\partial_t(h) = \Gamma\quad \circ(\star, d, \star)(d(h)\quad \wedge \quad|d(h)^\sharp|^{n-1} \quad \wedge \quad (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
-  Γ::Form1
+  Γ::Form0
   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
 
 to_graphviz(halfar_eq2)
@@ -65,8 +63,7 @@ And here, a formulation of Glen's law from J.W. Glen's 1958 ["The flow law of ic
 
 ``` @example DEC
 glens_law = @decapode begin
-  #Γ::Form0
-  Γ::Form1
+  Γ::Form0
   (A,ρ,g,n)::Constant
   
   Γ == (2/(n+2))*A*(ρ*g)^n

docs/src/navier_stokes/ns.jl

@@ -221,7 +221,7 @@ function generate(s, my_symbol; hodge=GeometricHodge())
     :dinv => x -> fd0 \ x
     :♭♯ => x -> ♭♯_m * x
     :dᵦ => x -> dᵦ * x
-    _ => default_dec_matrix_generate(s, my_symbol, hodge)
+    _ => error("Unmatched operator $my_symbol")
   end
   return (args...) -> op(args...)
 end;
@@ -536,10 +536,10 @@ if PLOT_FCS
     ax = CairoMakie.Axis(fig[1,1],
       title="Vorticity along θ=0.4",
       xlabel="Azimuth Angle, φ")
-    lns_ic = CairoMakie.lines!(ax, 
+    lns_ic = CairoMakie.lines!(ax,
                                range(0.0,2π; length=length(in_lat_range)),
                                form_ic.data[in_lat_range][by_phi])
-    lns_fc = CairoMakie.lines!(ax, 
+    lns_fc = CairoMakie.lines!(ax,
                                range(0.0,2π; length=length(in_lat_range)),
                                form_fc.data[in_lat_range][by_phi])
     Legend(fig[1,2], [lns_ic, lns_fc], ["T=0.0", "T=$(tₑ)"])
@@ -570,7 +570,7 @@ if PLOT_FCS
     plot_vort_fc())
 
   azimuth_name = "azimuth.png"
-  save(joinpath(DATA_DR, azimuth_name), 
+  save(joinpath(DATA_DR, azimuth_name),
     azimuth_form(0.4, 0.01, VForm(s0inv*soln(0).d𝐮), VForm(s0inv*soln(tₑ).d𝐮)))
 
   diff_name = "relsolchng.png"
@@ -586,4 +586,3 @@ if PLOT_FCS
   end
 end
 end # Plot FCs
-

docs/src/nhs/nhs_full.md

@@ -175,7 +175,7 @@ subdivide_duals!(sd, Barycenter())
 
 function generate(sd, my_symbol; hodge=GeometricHodge())
   op = @match my_symbol begin
-    _ => default_dec_matrix_generate(sd, my_symbol, hodge) # TODO: Don't do this!
+    _ => error("Unmatched operator $my_symbol")
   end
   return op
 end

src/operators.jl

@@ -5,43 +5,28 @@ using Krylov
 using LinearAlgebra
 using SparseArrays
 
-function default_dec_cu_matrix_generate() end;
-
-function default_dec_matrix_generate(fs::PrimalGeometricMapSeries, my_symbol::Symbol, hodge::DiscreteHodge)
+# Define mappings for default DEC operations that are not optimizable.
+# --------------------------------------------------------------------
+function default_dec_generate(fs::PrimalGeometricMapSeries, my_symbol::Symbol, hodge::DiscreteHodge)
   op = @match my_symbol begin
     # Inverse Laplacians
     :Δ₀⁻¹ => dec_Δ⁻¹(Val{0}, fs)
-    _ => default_dec_matrix_generate(finest_mesh(fs), my_symbol, hodge)
+    _ => default_dec_generate(finest_mesh(fs), my_symbol, hodge)
   end
 end
 
-function default_dec_matrix_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::DiscreteHodge)
-  op = @match my_symbol begin
-
-    # Regular Hodge Stars
-    :⋆₀ => dec_mat_hodge(0, sd, hodge)
-    :⋆₁ => dec_mat_hodge(1, sd, hodge)
-    :⋆₂ => dec_mat_hodge(2, sd, hodge)
-
-    # Inverse Hodge Stars
-    :⋆₀⁻¹ => dec_mat_inverse_hodge(0, sd, hodge)
-    :⋆₁⁻¹ => dec_pair_inv_hodge(Val{1}, sd, hodge) # Special since Geo is a solver
-    :⋆₂⁻¹ => dec_mat_inverse_hodge(1, sd, hodge)
+function default_dec_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::DiscreteHodge=GeometricHodge())
 
-    # Differentials
-    :d₀ => dec_mat_differential(0, sd)
-    :d₁ => dec_mat_differential(1, sd)
+  op = @match my_symbol begin
 
-    # Dual Differentials
-    :dual_d₀ || :d̃₀ => dec_mat_dual_differential(0, sd)
-    :dual_d₁ || :d̃₁ => dec_mat_dual_differential(1, sd)
+    # :plus => (+)
+    :(-) || :neg => x -> -1 .* x
+    :ln => (x -> log.(x))
 
-    # Wedge Products
-    :∧₀₁ => dec_pair_wedge_product(Tuple{0,1}, sd)
-    :∧₁₀ => dec_pair_wedge_product(Tuple{1,0}, sd)
-    :∧₀₂ => dec_pair_wedge_product(Tuple{0,2}, sd)
-    :∧₂₀ => dec_pair_wedge_product(Tuple{2,0}, sd)
-    :∧₁₁ => dec_pair_wedge_product(Tuple{1,1}, sd)
+    # Musical Isomorphisms
+    :♯ᵖᵖ => dec_♯_p(sd)
+    :♯ᵈᵈ => dec_♯_d(sd)
+    :♭ᵈᵖ => dec_♭(sd)
 
     # Primal-Dual Wedge Products
     :∧ᵖᵈ₁₁ => dec_wedge_product_pd(Tuple{1,1}, sd)
@@ -68,16 +53,55 @@ function default_dec_matrix_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::
     # Inverse Laplacians
     :Δ₀⁻¹ => dec_inv_lap_solver(Val{0}, sd)
 
-    # Musical Isomorphisms
-    :♯ => dec_♯_p(sd)
-    :♯ᵈ => dec_♯_d(sd)
+    _ => error("Unmatched operator $my_symbol")
+  end
+
+  return (args...) -> op(args...)
+end
+
+function default_dec_cu_generate() end;
+
+# Define mappings for default DEC operations that are optimizable.
+# ----------------------------------------------------------------
+function default_dec_cu_matrix_generate() end;
+
+function default_dec_matrix_generate(fs::PrimalGeometricMapSeries, my_symbol::Symbol, hodge::DiscreteHodge)
+  op = @match my_symbol begin
+    _ => default_dec_matrix_generate(finest_mesh(fs), my_symbol, hodge)
+  end
+end
+
+function default_dec_matrix_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::DiscreteHodge)
+  op = @match my_symbol begin
+
+    # Regular Hodge Stars
+    :⋆₀ => dec_mat_hodge(0, sd, hodge)
+    :⋆₁ => dec_mat_hodge(1, sd, hodge)
+    :⋆₂ => dec_mat_hodge(2, sd, hodge)
+
+    # Inverse Hodge Stars
+    :⋆₀⁻¹ => dec_mat_inverse_hodge(0, sd, hodge)
+    :⋆₁⁻¹ => dec_pair_inv_hodge(Val{1}, sd, hodge) # Special since Geo is a solver
+    :⋆₂⁻¹ => dec_mat_inverse_hodge(1, sd, hodge)
+
+    # Differentials
+    :d₀ => dec_mat_differential(0, sd)
+    :d₁ => dec_mat_differential(1, sd)
+
+    # Dual Differentials
+    :dual_d₀ || :d̃₀ => dec_mat_dual_differential(0, sd)
+    :dual_d₁ || :d̃₁ => dec_mat_dual_differential(1, sd)
 
-    :♭ => dec_♭(sd)
+    # Wedge Products
+    :∧₀₁ => dec_pair_wedge_product(Tuple{0,1}, sd)
+    :∧₁₀ => dec_pair_wedge_product(Tuple{1,0}, sd)
+    :∧₀₂ => dec_pair_wedge_product(Tuple{0,2}, sd)
+    :∧₂₀ => dec_pair_wedge_product(Tuple{2,0}, sd)
+    :∧₁₁ => dec_pair_wedge_product(Tuple{1,1}, sd)
 
     # Averaging Operator
     :avg₀₁ => dec_avg₀₁(sd)
 
-    :neg => x -> -1 .* x
      _ => error("Unmatched operator $my_symbol")
   end
 
@@ -162,20 +186,6 @@ function dec_inv_lap_solver(::Type{Val{0}}, sd::HasDeltaSet)
   x -> inv_lap \ x
 end
 
-function default_dec_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::DiscreteHodge=GeometricHodge())
-
-  op = @match my_symbol begin
-
-    :plus => (+)
-    :(-) || :neg => x -> -1 .* x
-    :ln => (x -> log.(x))
-
-    _ => error("Unmatched operator $my_symbol")
-  end
-
-  return (args...) -> op(args...)
-end
-
 function open_operators(d::SummationDecapode; kwargs...)
   e = deepcopy(d)
   open_operators!(e, kwargs...)