diff --git a/Project.toml b/Project.toml
index 649cd9b4..c9f6d7eb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,26 +5,15 @@ version = "0.5.3"
 
 [deps]
 ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8"
-Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
-DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiagrammaticEquations = "6f00c28b-6bed-4403-80fa-30e0dc12f317"
-Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
-FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
-JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
-Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-Unicode = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
 
 [weakdeps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -34,36 +23,36 @@ DecapodesCUDAExt = "CUDA"
 
 [compat]
 ACSets = "0.2"
-Catlab = "0.15, 0.16"
+Aqua = "0.8"
+CUDA = "5.2"
 CombinatorialSpaces = "0.6.3"
 ComponentArrays = "0.15"
-CUDA = "5.2"
-DataStructures = "0.18.13"
 DiagrammaticEquations = "0.1"
 Distributions = "0.25"
-FileIO = "1.16"
 GeometryBasics = "0.4.2"
-JLD2 = "0.4"
 JSON = "0.21"
 Krylov = "0.9.6"
 LinearAlgebra = "1.9"
 MLStyle = "0.4.17"
+Markdown = "1.9"
 OrdinaryDiffEq = "6.47"
 PreallocationTools = "0.4"
-Requires = "1.3"
+Random = "1.9"
 SparseArrays = "1.9"
-StaticArrays = "1.7"
-Unicode = "1.9"
+Statistics = "1.9"
+Test = "1.9"
 julia = "1.9"
 
 [extras]
-AlgebraicPetri = "4f99eebe-17bf-4e98-b6a1-2c4f205a959b"
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["AlgebraicPetri", "CUDA", "Distributions", "Test", "Random", "GeometryBasics", "Statistics"]
+test = ["Aqua", "CUDA", "Distributions", "GeometryBasics", "JSON", "OrdinaryDiffEq", "Random", "Statistics", "Test"]
diff --git a/docs/make.jl b/docs/make.jl
index 5968bc13..94f8e1dd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -42,7 +42,8 @@ makedocs(
   checkdocs = :none,
   pagesonly = true,
   linkcheck = true,
-  linkcheck_ignore = [r"agupubs.onlinelibrary.wiley.com"], # This gives a 403 Forbidden
+  linkcheck_ignore = [r"agupubs\.onlinelibrary\.wiley\.com", # This gives a 403 Forbidden
+                      r"Decapodes\.jl/dev"], # 404, probably due to bad self-rerference
   pages     = Any[
     "Decapodes.jl" => "index.md",
     "Overview" => "overview/overview.md",
diff --git a/src/coordinates.jl b/examples/sw/coordinates.jl
similarity index 100%
rename from src/coordinates.jl
rename to examples/sw/coordinates.jl
diff --git a/test/coordinates.jl b/examples/sw/coordinates_test.jl
similarity index 100%
rename from test/coordinates.jl
rename to examples/sw/coordinates_test.jl
diff --git a/examples/sw/pressure.jl b/examples/sw/pressure.jl
index 806fec05..d0dc357c 100644
--- a/examples/sw/pressure.jl
+++ b/examples/sw/pressure.jl
@@ -62,7 +62,7 @@ function generate(sd, my_symbol; hodge=GeometricHodge())
   return (args...) ->  op(args...)
 end
 
-#include("coordinates.jl")
+include("coordinates.jl")
 
 radius = 6371+90
 
diff --git a/examples/sw/sw.jl b/examples/sw/sw.jl
index 9db62ced..db7edca7 100644
--- a/examples/sw/sw.jl
+++ b/examples/sw/sw.jl
@@ -43,7 +43,7 @@ ddp = SummationDecapode(diffExpr)
 gensim(expand_operators(ddp), [:C])
 f = eval(gensim(expand_operators(ddp), [:C]))
 
-#include("coordinates.jl")
+include("coordinates.jl")
 
 const RADIUS = 6371+90
 #primal_earth, npi, spi = makeSphere(0, 180, 5, 0, 360, 5, RADIUS);
diff --git a/ext/DecapodesCUDAExt.jl b/ext/DecapodesCUDAExt.jl
index f15e40e8..29f68809 100644
--- a/ext/DecapodesCUDAExt.jl
+++ b/ext/DecapodesCUDAExt.jl
@@ -2,7 +2,6 @@ module DecapodesCUDAExt
 using CombinatorialSpaces
 using LinearAlgebra
 using Base.Iterators
-using Catlab
 using Krylov
 using CUDA
 using CUDA.CUSPARSE
@@ -109,6 +108,10 @@ function dec_cu_pair_wedge_product(::Type{Tuple{1,1}}, sd::HasDeltaSet2D)
     (α, β) -> dec_c_wedge_product(Tuple{1,1}, α, β, val_pack, Val{:CUDA}))
 end
 
+function dec_pair_wedge_product(::Type{Tuple{0,0}}, sd::HasDeltaSet)
+  error("Replace me in compiled code with element-wise multiplication (.*)")
+end
+
 # TODO: These need to be converted into CuArrays/kernels
 function dec_cu_sharp_p(sd::HasDeltaSet2D)
   ♯_m = ♯_mat(sd, AltPPSharp())
diff --git a/src/Decapodes.jl b/src/Decapodes.jl
index 6b7faa81..b50fbea9 100644
--- a/src/Decapodes.jl
+++ b/src/Decapodes.jl
@@ -1,29 +1,15 @@
 module Decapodes
 
-using Catlab
-using Catlab.Theories
-using Catlab.Programs
-using Catlab.CategoricalAlgebra
-using Catlab.WiringDiagrams
-using Catlab.WiringDiagrams.DirectedWiringDiagrams
-using Catlab.ACSetInterface
-using MLStyle
-using Base.Iterators
-using SparseArrays
-using PreallocationTools
-
+using ACSets
 using DiagrammaticEquations
 using DiagrammaticEquations.Deca
+using MLStyle
 
 export
-findname, flat_op,
-gensim, evalsim, closest_point, findnode, compile, compile_env,  PhysicsState,  default_dec_matrix_generate, default_dec_cu_matrix_generate, default_dec_generate, VectorForm,
-CartesianPoint, SpherePoint, r, theta, phi, TangentBasis, θhat, ϕhat,
+gensim, evalsim, compile, compile_env, default_dec_matrix_generate, default_dec_cu_matrix_generate, default_dec_generate,
 CPUTarget, CUDATarget
 
-append_dot(s::Symbol) = Symbol(string(s)*'\U0307')
 
-include("coordinates.jl")
 include("operators.jl")
 include("simulation.jl")
 
diff --git a/src/operators.jl b/src/operators.jl
index 600d29c3..83ece938 100644
--- a/src/operators.jl
+++ b/src/operators.jl
@@ -1,10 +1,11 @@
+using Base.Iterators
 using CombinatorialSpaces
+using Krylov
 using LinearAlgebra
-using Base.Iterators
-using Catlab
+using SparseArrays
 
 function default_dec_cu_matrix_generate() end;
-  
+
 function default_dec_matrix_generate(sd, my_symbol, hodge)
   op = @match my_symbol begin
 
@@ -117,6 +118,10 @@ function dec_pair_wedge_product(::Type{Tuple{1,1}}, sd::HasDeltaSet2D)
     (α, β) -> dec_c_wedge_product(Tuple{1,1}, α, β, val_pack))
 end
 
+function dec_pair_wedge_product(::Type{Tuple{0,0}}, sd::HasDeltaSet)
+  error("Replace me in compiled code with element-wise multiplication (.*)")
+end
+
 function dec_♯_p(sd::HasDeltaSet2D)
   ♯_m = ♯_mat(sd, AltPPSharp())
   x -> ♯_m * x
diff --git a/src/simulation.jl b/src/simulation.jl
index b428285f..4354236c 100644
--- a/src/simulation.jl
+++ b/src/simulation.jl
@@ -1,13 +1,9 @@
-using Catlab
+using Base.Iterators
 using CombinatorialSpaces
 using ComponentArrays
-using OrdinaryDiffEq
-using GeometryBasics
 using LinearAlgebra
-using Base.Iterators
-using Catlab
 using MLStyle
-import Catlab.Programs.GenerateJuliaPrograms: compile
+using PreallocationTools
 
 const gensim_in_place_stub = Symbol("GenSim-M")
 
@@ -18,7 +14,7 @@ struct CUDATarget <: GenerationTarget end
 
 abstract type AbstractCall end
 
-struct UnaryCall <: AbstractCall 
+struct UnaryCall <: AbstractCall
   operator
   equality
   input
@@ -44,7 +40,7 @@ Base.Expr(c::UnaryCall) = begin
   end
 end
 
-struct BinaryCall <: AbstractCall 
+struct BinaryCall <: AbstractCall
   operator
   equality
   input1
@@ -65,7 +61,7 @@ Base.Expr(c::BinaryCall) = begin
   return Expr(c.equality, c.output, Expr(:call, c.operator, c.input1, c.input2))
 end
 
-struct VarargsCall <: AbstractCall 
+struct VarargsCall <: AbstractCall
   operator
   equality
   inputs
@@ -79,7 +75,7 @@ Base.Expr(c::VarargsCall) = begin
   return Expr(c.equality, c.output, Expr(:call, c.operator, c.inputs...))
 end
 
-struct AllocVecCall <: AbstractCall 
+struct AllocVecCall <: AbstractCall
   name
   form
   dimension
@@ -141,7 +137,7 @@ end =#
 
 function is_form(d::SummationDecapode, var_id::Int)
   type = d[var_id, :type]
-  return (type == :Form0 || type == :Form1 || type == :Form2 || 
+  return (type == :Form0 || type == :Form1 || type == :Form2 ||
     type == :DualForm0 || type == :DualForm1 || type == :DualForm2)
 end
 
@@ -345,10 +341,10 @@ function compile(d::SummationDecapode, inputs::Vector, alloc_vectors::Vector{All
             operator = promote_arithmetic_map[operator]
             equality = promote_arithmetic_map[equality]
             push!(alloc_vectors, AllocVecCall(rname, d[r, :type], dimension, stateeltype, code_target))
-          
+
           # TODO: Do we want to support the ability of a user to use the backslash operator?
           elseif(operator == :(*) || operator == :(/) || operator == :.* || operator == :./)
-            # WARNING: This part may break if we add more compiler types that have different 
+            # WARNING: This part may break if we add more compiler types that have different
             # operations for basic and broadcast modes, e.g. matrix multiplication vs broadcast
             if(!is_infer(d, arg1) && !is_infer(d, arg2))
               operator = promote_arithmetic_map[operator]
@@ -547,7 +543,7 @@ function gensim(user_d::AbstractNamedDecapode, input_vars; dimension::Int=2, sta
   infer_overload_compiler!(d′, dimension)
 
   # This will generate all of the fundemental DEC operators present
-  optimizable_dec_operators = Set([:⋆₀, :⋆₁, :⋆₂, :⋆₀⁻¹, :⋆₂⁻¹, 
+  optimizable_dec_operators = Set([:⋆₀, :⋆₁, :⋆₂, :⋆₀⁻¹, :⋆₂⁻¹,
                                   :d₀, :d₁, :dual_d₀, :d̃₀, :dual_d₁, :d̃₁])
   extra_dec_operators = Set([:⋆₁⁻¹, :∧₀₁, :∧₁₀, :∧₁₁, :∧₀₂, :∧₂₀])
 
@@ -580,19 +576,19 @@ function gensim(user_d::AbstractNamedDecapode, input_vars; dimension::Int=2, sta
   end
 end
 
-gensim(c::Collage; dimension::Int=2) = 
+gensim(c::Collage; dimension::Int=2) =
 gensim(collate(c); dimension=dimension)
 
 """    function gensim(d::AbstractNamedDecapode; dimension::Int=2)
 
 Generate a simulation function from the given Decapode. The returned function can then be combined with a mesh and a function describing function mappings to return a simulator to be passed to `solve`.
 """
-gensim(d::AbstractNamedDecapode; dimension::Int=2, stateeltype = Float64, code_target = CPUTarget()) = 
+gensim(d::AbstractNamedDecapode; dimension::Int=2, stateeltype = Float64, code_target = CPUTarget()) =
   gensim(d, vcat(collect(infer_state_names(d)), d[incident(d, :Literal, :type), :name]), dimension=dimension, stateeltype=stateeltype, code_target=code_target)
 
-evalsim(d::AbstractNamedDecapode; dimension::Int=2, stateeltype = Float64, code_target = CPUTarget()) = 
+evalsim(d::AbstractNamedDecapode; dimension::Int=2, stateeltype = Float64, code_target = CPUTarget()) =
   eval(gensim(d, dimension=dimension, stateeltype=stateeltype, code_target=code_target))
-evalsim(d::AbstractNamedDecapode, input_vars; dimension::Int=2, stateeltype = Float64, code_target = CPUTarget()) = 
+evalsim(d::AbstractNamedDecapode, input_vars; dimension::Int=2, stateeltype = Float64, code_target = CPUTarget()) =
   eval(gensim(d, input_vars, dimension=dimension, stateeltype=stateeltype, code_target=code_target))
 
 """
diff --git a/test/aqua.jl b/test/aqua.jl
new file mode 100644
index 00000000..25bd1ce9
--- /dev/null
+++ b/test/aqua.jl
@@ -0,0 +1,7 @@
+using Aqua, Decapodes
+
+Aqua.test_ambiguities(Decapodes)
+Aqua.test_all(
+  Decapodes, ambiguities=false,
+  deps_compat=(ignore=[:Markdown, :Random, :Test],),
+)
diff --git a/test/componentarrays.jl b/test/componentarrays.jl
index 47137b6a..4b066284 100644
--- a/test/componentarrays.jl
+++ b/test/componentarrays.jl
@@ -1,19 +1,14 @@
-using ComponentArrays
-using OrdinaryDiffEq
-using GeometryBasics
-using JSON
-using Distributions
-
-using Catlab
-using Catlab.CategoricalAlgebra
+using ACSets
 using CombinatorialSpaces
-
-using DiagrammaticEquations
-using DiagrammaticEquations.Deca
+using ComponentArrays
 using Decapodes
-using Test
-using MLStyle
+using DiagrammaticEquations
+using Distributions
+using GeometryBasics: Point3
 using LinearAlgebra
+using MLStyle
+using OrdinaryDiffEq
+using Test
 
 C = ones(Float64, 10)
 V = ones(Float64, 100)
@@ -173,4 +168,3 @@ soln = solve(prob, Tsit5())
 # record(fig, "diff_adv.gif", range(0.0, tₑ; length=150); framerate = 30) do t
 #     ob.color = findnode(soln(t), :C)[point_map]
 # end
-
diff --git a/test/cuda_sims.jl b/test/cuda_sims.jl
index dca520ad..5947aa2d 100644
--- a/test/cuda_sims.jl
+++ b/test/cuda_sims.jl
@@ -1,17 +1,15 @@
-using Catlab
+using ACSets
+using CUDA, CUDA.CUSPARSE
+using CombinatorialSpaces
+using ComponentArrays
 using Decapodes
 using DiagrammaticEquations
-using CombinatorialSpaces
-using GeometryBasics
+using Distributions
+using GeometryBasics: Point2, Point3
+using LinearAlgebra
 using MLStyle
-using ComponentArrays
 using OrdinaryDiffEq
-using CUDA
-using CUDA.CUSPARSE
-using LinearAlgebra
-using SparseArrays
 using Statistics
-using Distributions
 using Test
 Point2D = Point2{Float64}
 Point3D = Point3{Float64}
@@ -50,13 +48,13 @@ end
   # CUDA Setup and Solve
   cuda_sim = eval(gensim(Heat, code_target=CUDATarget()))
   cuda_fₘ = cuda_sim(sd, nothing, DiagonalHodge())
-  
+
   cuda_u₀ = ComponentArray(U=CuArray{Float64}(U))
   prob = ODEProblem(cuda_fₘ, cuda_u₀, (0, tₑ), constants_and_parameters)
   cuda_soln = solve(prob, Tsit5(), save_everystep=false)
 
   @test all(isapprox(cpu_soln(tₑ).U, Array(cuda_soln(tₑ).U); atol=1e-12))
-  @test RMSE(cpu_soln(tₑ).U, Array(cuda_soln(tₑ).U)) < 1e-13 
+  @test RMSE(cpu_soln(tₑ).U, Array(cuda_soln(tₑ).U)) < 1e-13
 end
 
 @testset "Heat Equation Float32" begin
@@ -88,7 +86,7 @@ end
   # CUDA Setup and Solve
   cuda_sim = eval(gensim(Heat, code_target=CUDATarget(), stateeltype=Float32))
   cuda_fₘ = cuda_sim(sd, nothing, DiagonalHodge())
-  
+
   cuda_u₀ = ComponentArray(U=CuArray{Float32}(U))
   prob = ODEProblem(cuda_fₘ, cuda_u₀, (0, tₑ), constants_and_parameters)
   cuda_soln = solve(prob, Tsit5(), save_everystep=false)
@@ -125,20 +123,20 @@ end
   U = map(sd[:point]) do (_,y,_)
     abs(y)
   end
-  
+
   V = map(sd[:point]) do (x,_,_)
     abs(x)
   end
-  
+
   # TODO: Try making this sparse.
   F₁ = map(sd[:point]) do (_,_,z)
     z ≥ 0.8 ? 5.0 : 0.0
   end
-    
+
   F₂ = zeros(Float64, nv(sd))
 
   α = 0.001
-  
+
   # CPU Setup and Solve
   cpu_sim = eval(gensim(Brusselator))
   cpu_fₘ = cpu_sim(sd, nothing, DiagonalHodge())
@@ -147,7 +145,7 @@ end
 
   cpu_constants_and_parameters = (
     α = α,
-    F = t -> t ≥ 1.1 ? F₂ : F₁)  
+    F = t -> t ≥ 1.1 ? F₂ : F₁)
 
   prob = ODEProblem(cpu_fₘ, cpu_u₀, (0, tₑ), cpu_constants_and_parameters)
   cpu_soln = solve(prob, Tsit5())
@@ -155,20 +153,20 @@ end
   # CUDA Setup and Solve
   cuda_sim = eval(gensim(Brusselator, code_target=CUDATarget()))
   cuda_fₘ = cuda_sim(sd, nothing, DiagonalHodge())
-  
+
   cuda_u₀ = ComponentArray(U=CuArray(U), V=CuArray(V))
 
   cuda_F₁ = CuArray(F₁)
-   
+
   cuda_F₂ = CuArray(F₂)
 
   cuda_constants_and_parameters = (
     α = α,
-    F = t -> t ≥ 1.1 ? cuda_F₂ : cuda_F₁)  
+    F = t -> t ≥ 1.1 ? cuda_F₂ : cuda_F₁)
 
   prob = ODEProblem(cuda_fₘ, cuda_u₀, (0, tₑ), cuda_constants_and_parameters)
   cuda_soln = solve(prob, Tsit5())
-  
+
   @test all(isapprox(cpu_soln(tₑ).U, Array(cuda_soln(tₑ).U); atol=1e-11))
   @test RMSE(cpu_soln(tₑ).U, Array(cuda_soln(tₑ).U)) < 1e-13
 end
@@ -178,11 +176,11 @@ end
     (n,w)::DualForm0
     dX::Form1
     (a,ν,m)::Constant
-  
+
     ∂ₜ(w) == a - w - w * n^2 + ν * L(dX, w)
     ∂ₜ(n) == w * n^2 - m*n + Δ(n)
   end
-  
+
   Klausmeier[9, :type] = :DualForm0
   Klausmeier[10, :type] = :DualForm0
   Klausmeier[15, :type] = :DualForm0
@@ -263,11 +261,11 @@ end
     (n,w)::DualForm0
     dX::Form1
     (a,ν,m)::Constant
-  
+
     ∂ₜ(w) == a - w - w * n^2 + ν * L(dX, w)
     ∂ₜ(n) == w * n^2 - m*n + Δ(n)
   end
-  
+
   Klausmeier[9, :type] = :DualForm0
   Klausmeier[10, :type] = :DualForm0
   Klausmeier[15, :type] = :DualForm0
@@ -341,4 +339,3 @@ end
 
   @test RMSE(cpu_soln(tₑ).n, Array(cuda_soln(tₑ).n)) < 1e-5
 end
-
diff --git a/test/runtests.jl b/test/runtests.jl
index f96d6cd2..7c52dc8b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,7 @@
 using Test
 
-@testset "Coordinates" begin
-  include("coordinates.jl")
+@testset "Code Quality (Aqua.jl)" begin
+  include("aqua.jl")
 end
 
 @testset "ComponentArrays.jl Integration" begin
diff --git a/test/simulation.jl b/test/simulation.jl
index 2e73587e..7ddd2bf7 100644
--- a/test/simulation.jl
+++ b/test/simulation.jl
@@ -1,21 +1,14 @@
-using Decapodes
-import Decapodes: compile, gensim
-
-using Catlab
-using Catlab.CategoricalAlgebra
-
-using Test
-
-using MLStyle
+using ACSets
 using CombinatorialSpaces
-using GeometryBasics: Point3
-using LinearAlgebra
-using Distributions
 using ComponentArrays
-using OrdinaryDiffEq
+using Decapodes
 using DiagrammaticEquations
-using DiagrammaticEquations.Deca
-using GeometryBasics
+using Distributions
+using GeometryBasics: Point2, Point3
+using LinearAlgebra
+using MLStyle
+using OrdinaryDiffEq
+using Test
 
 function test_hodge(k, sd::HasDeltaSet, hodge)
   hodge = ⋆(k,sd,hodge=hodge)
@@ -187,7 +180,7 @@ f = eval(gensim(expand_operators(ddp)))
 fₘₚ = f(torus, generate)
 
 @test norm(fₘₛ(du, u₀, (k=2.0,), 0)  - fₘₚ(du, u₀, (k=t->2.0,), 0)) < 1e-4
- 
+
 # to solve the ODE over a duration, use the ODEProblem from OrdinaryDiffEq
 # tₑ = 10
 # using OrdinaryDiffEq
@@ -221,7 +214,7 @@ flatten(vfield::Function, mesh) =  ♭(mesh, DualVectorField(vfield.(mesh[triang
   end
 
   Brusselator = SummationDecapode(parse_decapode(quote
-      (U, V)::Form0{X} 
+      (U, V)::Form0{X}
       (U2V, aTU)::Form0{X}
       (U̇, V̇)::Form0{X}
 
@@ -230,7 +223,7 @@ flatten(vfield::Function, mesh) =  ♭(mesh, DualVectorField(vfield.(mesh[triang
 
       U2V == (U .* U) .* V
       aTU == α * Δ(U)
-      
+
       U̇ == One + U2V - (4.4 * U) + aTU + F
       V̇ == (3.4 * U) - U2V + aTU
 
@@ -283,13 +276,13 @@ flatten(vfield::Function, mesh) =  ♭(mesh, DualVectorField(vfield.(mesh[triang
     U = map(earth[:point]) do (_,y,_)
       abs(y)
     end
-    
+
     V = map(earth[:point]) do (x,_,_)
       abs(x)
     end
 
     One = ones(nv(earth))
-        
+
     F₁ = map(earth[:point]) do (_,_,z)
       z ≥ 0.8 ? 5.0 : 0.0
     end
@@ -311,13 +304,13 @@ flatten(vfield::Function, mesh) =  ♭(mesh, DualVectorField(vfield.(mesh[triang
       U = map(earth[:point]) do (_,y,_)
         abs(y)
       end
-      
+
       V = map(earth[:point]) do (x,_,_)
         abs(x)
       end
 
       One = ones(nv(earth))
-          
+
       F₁ = map(earth[:point]) do (_,_,z)
         z ≥ 0.8 ? 5.0 : 0.0
       end
@@ -445,7 +438,7 @@ end
     D == d(d(C))
   end
   @test 4 == length(checkForContractionInGensim(single_contract))
-  
+
   sim = eval(gensim(contract_with_summation))
   f = sim(earth, default_dec_generate)
   A = 2 * ones(nv(earth))
@@ -472,7 +465,7 @@ end
     D == d(d(C))
   end
   @test 4 == length(checkForContractionInGensim(single_contract))
-  
+
   sim = eval(gensim(contract_with_op2))
   f = sim(earth, default_dec_generate)
   A = 3 * ones(nv(earth))
@@ -495,7 +488,7 @@ end
     D == ⋆(⋆(B))
   end
   @test 4 == length(checkForContractionInGensim(single_contract))
-  
+
   sim = eval(gensim(later_contraction))
   f = sim(earth, default_dec_generate)
   A = 4 * ones(nv(earth))
@@ -515,7 +508,7 @@ end
     D == d(A)
   end
   @test 0 == length(checkForContractionInGensim(no_contraction))
-  
+
   sim = eval(gensim(no_contraction))
   f = sim(earth, default_dec_generate)
   A = [i for i in 1:nv(earth)]
@@ -535,7 +528,7 @@ end
     D == d(k(A))
   end
   @test 0 == length(checkForContractionInGensim(no_unallowed))
-  
+
   sim = eval(gensim(no_unallowed))
 
   function generate_no_unallowed(sd, my_symbol; hodge=GeometricHodge())
@@ -569,7 +562,7 @@ end
 
     F == A ∧ (C ∧ B)
   end
-  
+
   sim = eval(gensim(wedges01))
 
   f = sim(earth, default_dec_generate)
@@ -583,19 +576,19 @@ end
   f(du, u, constants_and_parameters, 0)
 
   @test du.A == du.B == du.C
-  
+
   # Testing wedge 11 operators function
   wedges11 = @decapode begin
     (A, B)::Form1
     (D, E)::Form2
 
     D == ∂ₜ(A)
-    E == ∂ₜ(B)  
+    E == ∂ₜ(B)
 
     D == A ∧ B
     E == B ∧ A
   end
-  
+
   sim = eval(gensim(wedges11))
 
   f = sim(earth, default_dec_generate)
@@ -618,12 +611,12 @@ end
     (D, E)::Form2
 
     D == ∂ₜ(A)
-    E == ∂ₜ(B)  
+    E == ∂ₜ(B)
 
     D == A ∧ B
     E == B ∧ A
   end
-  
+
   sim = eval(gensim(wedges02))
 
   f = sim(earth, default_dec_generate)
@@ -646,7 +639,7 @@ end
     B == ∂ₜ(A)
     B == ⋆(⋆(A))
   end
-  
+
   sim = eval(gensim(GeoInvHodge1))
 
   f = sim(earth, default_dec_generate)
@@ -657,7 +650,7 @@ end
   constants_and_parameters = ()
   f(du, u, constants_and_parameters, 0)
 
-  @test all(isapprox.(du.A, -1 * ones(ne(earth))))  
+  @test all(isapprox.(du.A, -1 * ones(ne(earth))))
 
   # Testing Diagonal inverse hodge 1
   DiagonalInvHodge1 = @decapode begin
@@ -666,7 +659,7 @@ end
     B == ∂ₜ(A)
     B == ⋆(⋆(A))
   end
-  
+
   sim = eval(gensim(DiagonalInvHodge1))
 
   f = sim(earth, default_dec_generate, DiagonalHodge())
@@ -677,7 +670,7 @@ end
   constants_and_parameters = ()
   f(du, u, constants_and_parameters, 0)
 
-  @test all(isapprox.(du.A, -1 * ones(ne(earth))))  
+  @test all(isapprox.(du.A, -1 * ones(ne(earth))))
 
 end
 
@@ -698,16 +691,16 @@ end
     # 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)
-  
+
     # Test that no error is thrown here
-    f = sim(line, default_dec_generate, DiagonalHodge())  
+    f = sim(line, default_dec_generate, DiagonalHodge())
 end
 
 @testset "GenSim Compilation" begin
@@ -721,19 +714,19 @@ end
       end
     return op
   end
-  
+
   HeatTransfer = @decapode begin
-    (HT, Tₛ)::Form0      
-    (D, cosϕᵖ, cosϕᵈ)::Constant      
+    (HT, Tₛ)::Form0
+    (D, cosϕᵖ, cosϕᵈ)::Constant
     HT == (D ./ cosϕᵖ) .* (⋆)(d(cosϕᵈ .* (⋆)(d(Tₛ))))
   end
 
   sim_HT = evalsim(HeatTransfer)
   @test sim_HT(d_rect, generate, DiagonalHodge()) isa Any
 
-  Jordan_Kinderlehrer_Otto = @decapode begin 
-    (ρ, Ψ)::Form0       
-    β⁻¹::Constant          
+  Jordan_Kinderlehrer_Otto = @decapode begin
+    (ρ, Ψ)::Form0
+    β⁻¹::Constant
     ∂ₜ(ρ) == (∘(⋆, d, ⋆))(d(Ψ) ∧ ρ) + β⁻¹ * Δ(ρ)
   end
 
@@ -774,7 +767,7 @@ end
   @test sim_LJ(d_rect, generate, DiagonalHodge()) isa Any
 
   Tracer = @decapode begin
-    (c, C, F, c_up)::Form0   
+    (c, C, F, c_up)::Form0
     (v, V, q)::Form1
     c_up == (((-1 * (⋆)(L(v, (⋆)(c))) - (⋆)(L(V, (⋆)(c)))) - (⋆)(L(v, (⋆)(C)))) - (∘(⋆, d, ⋆))(q)) + F
   end
@@ -812,4 +805,3 @@ nallocs = @allocations f(du, u₀, p, (0,1.0))
   sizeof(p.D) +
   sizeof(Decapodes.FixedSizeDiffCache(Vector{Float64}(undef , nparts(sd, :V))))
 end
-