Use upstreamed averaging operator (#243)

docs/src/ascii.md

@@ -17,6 +17,7 @@ Some users may have trouble entering unicode characters like ⋆ or ∂ in their
 | ∧        | wedge      | wedge product, generalizing the cross product |
 | ⋆⁻¹      | star\_inv  | Hodge star, generalizing transpose            |
 | ∘(⋆,d,⋆) | div        | divergence, a.k.a. ∇⋅                         |
+| avg₀₁    | avg_01     | average values stored on endpoints of edges   |
 
 ## Vector Calculus Equivalents
 

docs/src/bsh/budyko_sellers_halfar.md

@@ -271,18 +271,6 @@ constants_and_parameters = (
 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
@@ -304,10 +292,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
       end
     end
     :mag => x -> norm.(x)
-    :avg₀₁ => begin
-      avg_mat = create_average_matrix(sd)
-      x -> avg_mat * x
-    end
     x => error("Unmatched operator $my_symbol")
   end
   return (args...) -> op(args...)

docs/src/ice_dynamics/ice_dynamics.md

@@ -166,18 +166,6 @@ constants_and_parameters = (
 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
@@ -199,10 +187,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
       end
     end
     :mag => x -> norm.(x)
-    :avg₀₁ => begin
-      avg_mat = create_average_matrix(sd)
-      x -> avg_mat * x
-    end
     x => error("Unmatched operator $my_symbol")
   end
   return (args...) -> op(args...)
@@ -360,10 +344,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
       x -> sharp_mat * x
     end
     :mag => x -> norm.(x)
-    :avg₀₁ => begin
-      avg_mat = create_average_matrix(sd)
-      x -> avg_mat * x
-    end
     x => error("Unmatched operator $my_symbol")
   end
   return (args...) -> op(args...)

examples/climate/budyko_sellers_halfar.jl

@@ -39,7 +39,6 @@ warming = @decapode begin
   (Tₛ)::Form0
   (A)::Form1
 
-  ## A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e-17)
   A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e7)
 
 end
@@ -160,18 +159,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
     :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

examples/diff_adv/cht.jl

@@ -269,19 +269,6 @@ begin
     mul!(x, matrices[:cross], matrices[:β_cache])
     #x .= matrices[:cross] * ((matrices[:t2c]*α).*(matrices[:e2c]*β))
   end
-
-
-  function avg_mat(::Type{Val{(0,1)}},s)
-    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
-    sparse(I,J,V)
-  end
 end
 
 # TODO: Produce the implementations of these DEC operators
@@ -449,7 +436,6 @@ begin
   g = ♭(sd, DualVectorField(gravity.(sd[triangle_center(sd),:dual_point]))).data;
   p = [density for p in s[:point]] * (288.15 * R₀)
 end
-m_avg = avg_mat(Val{(0,1)}, sd)
 # sim = eval(gensim(HeatXFer))
 
 wedge_cache = init_wedge_ops(sd)
@@ -463,7 +449,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
             x[cyl_inner] .= 0
             x
           end
-    :avg₀₁ => x -> m_avg * x
     :∂ₜₐ => (x) -> begin
       x[cyl_inner] .= 0
       x

examples/diff_adv/cht_cuda.jl

@@ -290,19 +290,6 @@ begin
     mul!(x, matrices[:cross], matrices[:β_cache])
     #x .= matrices[:cross] * ((matrices[:t2c]*α).*(matrices[:e2c]*β))
   end
-
-
-  function avg_mat(::Type{Val{(0,1)}},s)
-    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
-    CuSparseMatrixCSC(sparse(I,J,V))
-  end
 end
 
 ## Constants
@@ -435,7 +422,6 @@ begin
   g = CuArray(♭(sd, DualVectorField(gravity.(sd[triangle_center(sd),:dual_point]))).data);
   p = CuArray([density for p in s[:point]] * (288.15 * R₀))
 end
-m_avg = avg_mat(Val{(0,1)}, sd)
 # sim = eval(gensim(HeatXFer, code_target=CUDATarget()))
 
 wedge_cache = init_wedge_ops(sd)
@@ -453,7 +439,6 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
             x
       end
     end
-    :avg₀₁ => x -> m_avg * x
     :∂ₜₐ => (x) -> begin
       x[cyl_inner] .= 0
       x
@@ -515,6 +500,7 @@ function simulate(mesh, operators, hodge = GeometricHodge())
       ∂ᵥ = operators(mesh, :∂ᵥ)
       ∂ᵣ = operators(mesh, :∂ᵣ)
       ∂ₚ = operators(mesh, :∂ₚ)
+      # XXX: This auto-generated output was created before the averaging operator was upstreamed in CombinatorialSpaces PR 97.
       # avg₀₁ = operators(mesh, :avg₀₁)
       avg₀₁ = avg_mat(Val{(0,1)}, sd)
       ∂ₜₐ = operators(mesh, :∂ₜₐ)

ext/DecapodesCUDAExt.jl

@@ -36,6 +36,9 @@ function default_dec_cu_matrix_generate(sd, my_symbol, hodge)
     :∧₂₀ => dec_cu_pair_wedge_product(Tuple{2,0}, sd)
     :∧₁₁ => dec_cu_pair_wedge_product(Tuple{1,1}, sd)
 
+    # Averaging Operator
+    :avg₀₁ => dec_cu_avg₀₁(sd)
+
     _ => error("Unmatched operator $my_symbol")
   end
 
@@ -127,4 +130,11 @@ function dec_cu_flat(sd::HasDeltaSet2D)
   ♭_m = ♭_mat(sd)
   x -> ♭_m * x
 end
+
+function dec_cu_avg₀₁(sd::HasDeltaSet)
+  # TODO: Cast this in the CombinatorialSpaces CUDA extension directly.
+  avg_mat = CuSparseMatrixCSC(avg₀₁_mat(sd))
+  (avg_mat, x -> avg_mat * x)
+end
+
 end

src/operators.jl

@@ -62,6 +62,9 @@ function default_dec_matrix_generate(sd, my_symbol, hodge)
 
     :♭ => dec_♭(sd)
 
+    # Averaging Operator
+    :avg₀₁ => dec_avg₀₁(sd)
+
     :neg => x -> -1 .* x
      _ => error("Unmatched operator $my_symbol")
   end
@@ -137,6 +140,11 @@ function dec_♭(sd::HasDeltaSet2D)
   x -> ♭_m * x
 end
 
+function dec_avg₀₁(sd::HasDeltaSet)
+  avg_mat = avg₀₁_mat(sd)
+  (avg_mat, x -> avg_mat * x)
+end
+
 function default_dec_generate(sd, my_symbol, hodge=GeometricHodge())
 
   op = @match my_symbol begin

src/simulation.jl

@@ -544,7 +544,8 @@ function gensim(user_d::AbstractNamedDecapode, input_vars; dimension::Int=2, sta
 
   # This will generate all of the fundemental DEC operators present
   optimizable_dec_operators = Set([:⋆₀, :⋆₁, :⋆₂, :⋆₀⁻¹, :⋆₂⁻¹,
-                                  :d₀, :d₁, :dual_d₀, :d̃₀, :dual_d₁, :d̃₁])
+                                  :d₀, :d₁, :dual_d₀, :d̃₀, :dual_d₁, :d̃₁,
+                                  :avg₀₁])
   extra_dec_operators = Set([:⋆₁⁻¹, :∧₀₁, :∧₁₀, :∧₁₁, :∧₀₂, :∧₂₀])
 
   init_dec_matrices!(d′, dec_matrices, union(optimizable_dec_operators, extra_dec_operators))