Generate Musical Isomorphisms by Default

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
@@ -140,9 +140,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::Form0
 
-  A == avg₀₁(5.8282*10^(-0.236 * Tₛ)*1.65e7)
+  A == 5.8282*10^(-0.236 * Tₛ)*1.65e7
 
 end
 

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

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

src/operators.jl

@@ -5,6 +5,31 @@ using Krylov
 using LinearAlgebra
 using SparseArrays
 
+# Define mappings for default DEC operations that are not optimizable.
+# --------------------------------------------------------------------
+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))
+
+    # Musical Isomorphisms
+    :♯ᵖᵖ => dec_♯_p(sd)
+    :♯ᵈᵈ => dec_♯_d(sd)
+    :♭ᵈᵖ => dec_♭(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(sd::HasDeltaSet, my_symbol::Symbol, hodge::DiscreteHodge)
@@ -57,12 +82,6 @@ function default_dec_matrix_generate(sd::HasDeltaSet, my_symbol::Symbol, hodge::
     :Δᵈ₀ => Δᵈ(Val{0},sd)
     :Δᵈ₁ => Δᵈ(Val{1},sd)
 
-    # Musical Isomorphisms
-    :♯ => dec_♯_p(sd)
-    :♯ᵈ => dec_♯_d(sd)
-
-    :♭ => dec_♭(sd)
-
     # Averaging Operator
     :avg₀₁ => dec_avg₀₁(sd)
 
@@ -146,20 +165,6 @@ function dec_avg₀₁(sd::HasDeltaSet)
   (avg_mat, x -> avg_mat * 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; dimension::Int=2)
   e = deepcopy(d)
   open_operators!(e, dimension=dimension)

src/simulation.jl

@@ -223,6 +223,14 @@ end
 
 Base.showerror(io::IO, e::InvalidCodeTargetException) = print(io, "Provided code target $(e.code_target) is not yet supported in simulations")
 
+opt_generator_function(code_target::AbstractGenerationTarget) = throw(InvalidCodeTargetException(code_target))
+opt_generator_function(::CPUBackend) = :default_dec_matrix_generate
+opt_generator_function(::CUDABackend) = :default_dec_cu_matrix_generate
+
+generator_function(code_target::AbstractGenerationTarget) = throw(InvalidCodeTargetException(code_target))
+generator_function(::CPUBackend) = :default_dec_generate
+generator_function(::CUDABackend) = :default_dec_cu_generate
+
 """
     compile_env(d::SummationDecapode, dec_matrices::Vector{Symbol}, con_dec_operators::Set{Symbol}, code_target::AbstractGenerationTarget)
 
@@ -234,27 +242,30 @@ function compile_env(d::SummationDecapode, dec_matrices::Vector{Symbol}, con_dec
 
   defs = quote end
 
+  # These are optimizable default DEC functions.
   for op in dec_matrices
     op in defined_ops && continue
 
-    quote_op = QuoteNode(op)
-    mat_op = add_stub(GENSIM_INPLACE_STUB, op)
+    def = :(($(add_inplace_stub(op)), $op) = $(opt_generator_function(code_target))(mesh, $(QuoteNode(op)), hodge))
+    push!(defs.args, def)
 
-    # TODO: Add support for user-defined code targets
-    default_generation = @match code_target begin
-      ::CPUBackend => :default_dec_matrix_generate
-      ::CUDABackend => :default_dec_cu_matrix_generate
-      _ => throw(InvalidCodeTargetException(code_target))
-    end
+    push!(defined_ops, op)
+  end
+
+  operators = vcat(d[:op1], d[:op2])
+
+  # These are nonoptimizable default DEC functions.
+  for op in intersect(operators, non_optimizable(code_target))
+    op in defined_ops && continue
 
-    def = :(($mat_op, $op) = $(default_generation)(mesh, $quote_op, hodge))
+    def = :($op = $(generator_function(code_target))(mesh, $(QuoteNode(op)), hodge))
     push!(defs.args, def)
 
     push!(defined_ops, op)
   end
 
   # Add in user-defined operations
-  for op in vcat(d[:op1], d[:op2])
+  for op in operators
     if op == DerivOp || op in defined_ops || op in ARITHMETIC_OPS
       continue
     end
@@ -418,7 +429,7 @@ function compile(d::SummationDecapode, inputs::Vector{Symbol}, alloc_vectors::Ve
         if preallocate && is_form(d, t)
           if operator in optimizable_dec_operators
             equality = PROMOTE_ARITHMETIC_MAP[equality]
-            operator = add_stub(GENSIM_INPLACE_STUB, operator)
+            operator = add_inplace_stub(operator)
             push!(alloc_vectors, AllocVecCall(tname, d[t, :type], dimension, stateeltype, code_target))
 
           elseif operator == :(-) || operator == :.-
@@ -463,7 +474,7 @@ function compile(d::SummationDecapode, inputs::Vector{Symbol}, alloc_vectors::Ve
               push!(alloc_vectors, AllocVecCall(rname, d[r, :type], dimension, stateeltype, code_target))
             end
           elseif operator in optimizable_dec_operators
-            operator = add_stub(GENSIM_INPLACE_STUB, operator)
+            operator = add_inplace_stub(operator)
             equality = PROMOTE_ARITHMETIC_MAP[equality]
             push!(alloc_vectors, AllocVecCall(rname, d[r, :type], dimension, stateeltype, code_target))
           end
@@ -678,15 +689,28 @@ end
 
 Base.showerror(io::IO, e::UnsupportedStateeltypeException) = print(io, "Decapodes does not support state element types as $(e.type), only Float32 or Float64")
 
+const MATRIX_OPTIMIZABLE_DEC_OPERATORS = Set([:⋆₀, :⋆₁, :⋆₂, :⋆₀⁻¹, :⋆₂⁻¹,
+                                              :d₀, :d₁, :dual_d₀, :d̃₀, :dual_d₁, :d̃₁,
+                                              :avg₀₁])
+
+const NONMATRIX_OPTIMIZABLE_DEC_OPERATORS = Set([:⋆₁⁻¹, :∧₀₁, :∧₁₀, :∧₁₁, :∧₀₂, :∧₂₀])
+
+const NON_OPTIMIZABLE_CPU_OPERATORS = Set([:♯ᵖᵖ, :♯ᵈᵈ, :♭ᵈᵖ])
+const NON_OPTIMIZABLE_CUDA_OPERATORS = Set{Symbol}()
+
+non_optimizable(::AbstractGenerationTarget) = throw(InvalidCodeTargetException(code_target))
+non_optimizable(::CPUBackend) = NON_OPTIMIZABLE_CPU_OPERATORS
+non_optimizable(::CUDABackend) = NON_OPTIMIZABLE_CUDA_OPERATORS
+
 """
     gensim(user_d::SummationDecapode, input_vars::Vector{Symbol}; dimension::Int=2, stateeltype::DataType = Float64, code_target::AbstractGenerationTarget = CPUTarget(), preallocate::Bool = true)
 
-Generates the entire code body for the simulation function. The returned simulation function can then be combined with a mesh, provided by `CombinatorialSpaces`, and a function describing symbol 
+Generates the entire code body for the simulation function. The returned simulation function can then be combined with a mesh, provided by `CombinatorialSpaces`, and a function describing symbol
 to operator mappings to return a simulator that can be used to solve the represented equations given initial conditions.
-  
+
 **Arguments:**
-  
-`user_d`: The user passed Decapode for which simulation code will be generated. (This is not modified) 
+
+`user_d`: The user passed Decapode for which simulation code will be generated. (This is not modified)
 
 `input_vars` is the collection of variables whose values are known at the beginning of the simulation. (Defaults to all state variables and literals in the Decapode)
 
@@ -732,20 +756,14 @@ function gensim(user_d::SummationDecapode, input_vars::Vector{Symbol}; dimension
   open_operators!(gen_d, dimension = dimension)
   infer_overload_compiler!(gen_d, dimension)
 
-  # This will generate all of the fundemental DEC operators present
-  optimizable_dec_operators = Set([:⋆₀, :⋆₁, :⋆₂, :⋆₀⁻¹, :⋆₂⁻¹,
-                                  :d₀, :d₁, :dual_d₀, :d̃₀, :dual_d₁, :d̃₁,
-                                  :avg₀₁])
-  extra_dec_operators = Set([:⋆₁⁻¹, :∧₀₁, :∧₁₀, :∧₁₁, :∧₀₂, :∧₂₀])
-
-  init_dec_matrices!(gen_d, dec_matrices, union(optimizable_dec_operators, extra_dec_operators))
+  init_dec_matrices!(gen_d, dec_matrices, union(MATRIX_OPTIMIZABLE_DEC_OPERATORS, NONMATRIX_OPTIMIZABLE_DEC_OPERATORS))
 
   # This contracts matrices together into a single matrix
   contracted_dec_operators = Set{Symbol}()
-  contract_operators!(gen_d, white_list = optimizable_dec_operators)
+  contract_operators!(gen_d, white_list = MATRIX_OPTIMIZABLE_DEC_OPERATORS)
   cont_defs = link_contract_operators(gen_d, contracted_dec_operators, stateeltype, code_target)
 
-  union!(optimizable_dec_operators, contracted_dec_operators, extra_dec_operators)
+  optimizable_dec_operators = union(MATRIX_OPTIMIZABLE_DEC_OPERATORS, contracted_dec_operators, NONMATRIX_OPTIMIZABLE_DEC_OPERATORS)
 
   # Compilation of the simulation
   equations = compile(gen_d, input_vars, alloc_vectors, optimizable_dec_operators, dimension, stateeltype, code_target, preallocate)

test/simulation.jl

@@ -641,6 +641,8 @@ end
 
 end
 
+filter_lnn(arr::AbstractVector) = filter(x -> !(x isa LineNumberNode), arr)
+
 @testset "1-D Mat Generation" begin
   Point2D = Point2{Float64}
   function generate_dual_mesh(s::HasDeltaSet1D)
@@ -655,19 +657,20 @@ end
   add_edges!(primal_line, [1,2], [2,3])
   line = generate_dual_mesh(primal_line)
 
-    # Testing Diagonal inverse hodge 1
-    DiagonalInvHodge1 = @decapode begin
-      A::DualForm1
+  # Testing Diagonal inverse hodge 1
+  DiagonalInvHodge1 = @decapode begin
+    A::DualForm1
 
-      B == ∂ₜ(A)
-      B == ⋆(A)
-    end
-    g = gensim(DiagonalInvHodge1)
-    @test gensim(DiagonalInvHodge1).args[2].args[2].args[3].args[2].args[2].args[3].value == :⋆₁⁻¹
-    sim = eval(g)
+    B == ∂ₜ(A)
+    B == ⋆(A)
+  end
+  g = gensim(DiagonalInvHodge1)
+  @test g.args[2].args[2].args[3].args[2].args[2].args[3].value == :⋆₁⁻¹
+  @test length(filter_lnn(g.args[2].args[2].args[3].args)) == 1
+  sim = eval(g)
 
-    # Test that no error is thrown here
-    f = sim(line, default_dec_generate, DiagonalHodge())
+  # TODO: Error is being thrown here
+  @test f = sim(line, default_dec_generate, DiagonalHodge()) isa Any
 end
 
 @testset "GenSim Compilation" begin
@@ -829,4 +832,3 @@ haystack = string(gensim(LargeSum))
 @test occursin(needle, haystack)
 
 end
-