diff --git a/docs/src/bsh/budyko_sellers_halfar.md b/docs/src/bsh/budyko_sellers_halfar.md index 30d8b612..4d2dcfba 100644 --- a/docs/src/bsh/budyko_sellers_halfar.md +++ b/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 diff --git a/docs/src/cism/cism.md b/docs/src/cism/cism.md index b4d624cf..e10b192c 100644 --- a/docs/src/cism/cism.md +++ b/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) diff --git a/docs/src/grigoriev/grigoriev.md b/docs/src/grigoriev/grigoriev.md index 7bebba6b..076c5485 100644 --- a/docs/src/grigoriev/grigoriev.md +++ b/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 diff --git a/docs/src/halmo/halmo.md b/docs/src/halmo/halmo.md index 55c8522d..4b24e1ed 100644 --- a/docs/src/halmo/halmo.md +++ b/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 diff --git a/docs/src/ice_dynamics/ice_dynamics.md b/docs/src/ice_dynamics/ice_dynamics.md index 0cf499c1..e4780404 100644 --- a/docs/src/ice_dynamics/ice_dynamics.md +++ b/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 diff --git a/src/operators.jl b/src/operators.jl index 80935186..89b736f0 100644 --- a/src/operators.jl +++ b/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) diff --git a/src/simulation.jl b/src/simulation.jl index bb592e09..03d32f7f 100644 --- a/src/simulation.jl +++ b/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) diff --git a/test/simulation.jl b/test/simulation.jl index 1b6e4b5a..de3cbed5 100644 --- a/test/simulation.jl +++ b/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 -