diff --git a/docs/make.jl b/docs/make.jl index a265ade0..8b0c9be7 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -44,6 +44,7 @@ makedocs( "Decapodes.jl" => "index.md", "Overview" => "overview.md", "Equations" => "equations.md", + "BCs & ICs" => "collages.md", "ASCII Operators" => "ascii.md", "Misc Features" => "bc_debug.md", "Pipe Flow" => "poiseuille.md", diff --git a/docs/src/collages.md b/docs/src/collages.md new file mode 100644 index 00000000..08100c2f --- /dev/null +++ b/docs/src/collages.md @@ -0,0 +1,29 @@ +# Collages: Formalizing Boundary and Initial Conditions + +Representing physical equations via formal diagrams (i.e. a Decapode) solves many challenges in automatic simulation generation. In practice, however, setting initial and boundary conditions is a major part of the workflow, and this process itself must be formalized. We do so by implementing the notion of a *collage* of Decapodes. This is directly inspired by the [Diagrammatic view of differential equations in physics](https://arxiv.org/abs/2204.01843). See §4 for relevant theoretical discussion. + +## Layman's Definition: Semantics + +We can quickly introduce a *collage* as simply taking two Decapodes, specifying a *morphism* between the two (which point from quantities in one to quantities in the other), and interpreting this entire collection of arrows as itself a single large Decapode diagram. In contrast to the composition patterns used to specify composition between different physical components, we will not merge any quantities here. Rather, we interpret these new arrows as providing information on boundary and initial conditions. This information is passed directly to our simulation generator to help handle loading initial conditions and boundary conditions. + +## The Workflow + +### Goals +We have a few goals for this collage-based workflow: +- Make it quick and easy to try out different initial conditions +- Make it quick and easy to try out different boundary conditions +- Make it quick and easy to swap out different meshes +- Separate the specification of the "location" of boundaries from the "values" those boundaries should take +- Separate the specification of distributions from the properties of the mesh +- Formalize these representations for the sake of future-proofness and consistency + +And of course, if this workflow is not ergonomic for your particular circumstances, you can hand-craft your initial and boundary conditions, while still using the other simulation automation that the Decapodes framework offers. + +### Higher-Order Functions +To enforce a consistent workflow, we employ the programming technique of using higher-order functions. In particular, we will create a collage, and call `simulation_helper` on it. The simulation helper returns 3 things: a boundary condition **loader**, an initial condition **loader**, and a simulation **generator**. + +The boundary condition **loader** accepts a dictionary which maps the names of physical quantities to functions which evaluate them on an arbitrary mesh. This loader then returns a **generator** which takes in the mesh that you want to run your simulation on, and any other regular constants and parameters that you want passed to your simulation. When you execute the generator on a particular mesh, you get a named tuple of parameters that is suitable to be used with the `ODEProblem` interface. + +Similarly, the initial condition **loader** accepts a dictionary which maps the names of physical quantities to functions which evaluate them on an arbitrary mesh. This loader then returns a **generator** which takes in the mesh that you want to run your simulation on. When you execute the generator on a particular mesh, you get a state vector that is suitable to be used with the `ODEProblem` interface. + +The simulation **generator** is the output of `gensim`, with operators that apply your boundary conditions included for you. Executing the output of `gensim` on a mesh returns a simulation function that is suitable to be used with the `ODEProblem` interface. diff --git a/examples/chemistry/brusselator_bounded.jl b/examples/chemistry/brusselator_bounded.jl new file mode 100644 index 00000000..e650296c --- /dev/null +++ b/examples/chemistry/brusselator_bounded.jl @@ -0,0 +1,163 @@ +# Demonstrate how to encode boundary conditions in Decapodes using the "collage" +# technique. +using Catlab +using Catlab.Graphics +using CombinatorialSpaces +using CombinatorialSpaces.ExteriorCalculus +using Decapodes +using MultiScaleArrays +using MLStyle +using OrdinaryDiffEq +using LinearAlgebra +using GLMakie +using Logging +using JLD2 +using Printf + +using GeometryBasics: Point2, Point3 +Point2D = Point2{Float64} +Point3D = Point3{Float64} + +BrusselatorDynamics = @decapode begin + (U, V)::Form0 + U2V::Form0 + α::Constant + F::Parameter + U2V == (U .* U) .* V + ∂ₜ(U) == 1 + U2V - (4.4 * U) + (α * Δ(U)) + F + ∂ₜ(V) == (3.4 * U) - U2V + (α * Δ(U)) +end +# Visualize. You must have graphviz installed. +to_graphviz(BrusselatorDynamics) + +# This is a "discrete" Decapode, with no morphisms. +BrusselatorBoundaries = @decapode begin + L::Constant +end + +# Specify the BC morphism between Decapodes. +BrusselatorBCMorphism = BCMorphism(ACSetTransformation( + BrusselatorBoundaries, BrusselatorDynamics, + Var = [1])) + +# This is a "discrete" Decapode, with no morphisms. +BrusselatorInitialConditions = @decapode begin + (U₀, V₀)::Form0 +end + +# Specify the IC morphism between Decapodes. +BrusselatorICMorphism = ICMorphism(ACSetTransformation( + BrusselatorInitialConditions, BrusselatorDynamics, + Var = [1,2])) + +# Wrap these morphisms into a single collage. +BrusselatorCollage = Collage( + BrusselatorBCMorphism, + BrusselatorICMorphism) + +# Create the BC loader, IC loader and simulation generator. +bc_loader, ic_loader, sim = + simulation_helper(BrusselatorCollage) + +# Specify a function that finds boundaries and values to assign, generic to any +# mesh. +function left_wall_constant_1(sd) + min_x = minimum(x -> x[1], sd[:point]) + return ( + findall(x -> x[1] == min_x, sd[:point]), + ones(count(x -> x[1] == min_x, sd[:point]))) +end + +# Specify a function that assigns initial values, generic to any mesh. +# This could be an explicit function, as below, or a function that +# interpolates data from e.g. a netcdf file. +function parabola_in_y(sd) + map(sd[:point]) do (_,y) + 22 * (y *(1-y))^(3/2) + end +end +function parabola_in_x(sd) + map(sd[:point]) do (x,_) + 27 * (x *(1-x))^(3/2) + end +end + +# Specify some mesh. +s = loadmesh(Rectangle_30x10()) +scaling_mat = Diagonal([1/maximum(x->x[1], s[:point]), + 1/maximum(x->x[2], s[:point]), + 1.0]) +s[:point] = map(x -> scaling_mat*x, s[:point]) +s[:edge_orientation] = false +orient!(s) +# Visualize the mesh. +GLMakie.wireframe(s) +sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s) +subdivide_duals!(sd, Circumcenter()) + +# Demonstrate the boundary condition loader. +# This parameter data will be passed normally for demonstration. +F₁ = map(sd[:point]) do (x,y) + (x-0.3)^2 + (y-0.6)^2 ≤ (0.1)^2 ? 5.0 : 0.0 +end +F₂ = zeros(nv(sd)) + +bc_generator = bc_loader(Dict( + :L => left_wall_constant_1)) + +p = bc_generator(sd, + (α = 0.001, + F = t -> t ≥ 1.1 ? F₂ : F₁)) + +# Demonstrate the initial condition loader. +ic_generator = ic_loader(Dict( + :U₀ => parabola_in_x, + :V₀ => parabola_in_y)) +u₀ = ic_generator(sd) + +# Generate the simulation. +fₘ = eval(sim)(sd, default_dec_generate) + +# Create problem and run sim for t ∈ [0,tₑ). +tₑ = 11.5 + +@info("Precompiling Solver") +prob = ODEProblem(fₘ, u₀, (0, 1e-4), p) +soln = solve(prob, Tsit5()) +soln.retcode != :Unstable || error("Solver was not stable") +@info("Solving") +prob = ODEProblem(fₘ, u₀, (0, tₑ), p) +soln = solve(prob, Tsit5()) +@info("Done") + +@save "brusselator_bounded.jld2" soln + +# Visualize the final conditions. +GLMakie.mesh(s, color=findnode(soln(tₑ), :U), colormap=:jet) + +begin # BEGIN Gif creation +frames = 100 +# Initial frame +fig = GLMakie.Figure(resolution = (1200, 800)) +p1 = GLMakie.mesh(fig[1,2], s, color=findnode(soln(0), :U), colormap=:jet, colorrange=extrema(findnode(soln(0), :U))) +p2 = GLMakie.mesh(fig[1,4], s, color=findnode(soln(0), :V), colormap=:jet, colorrange=extrema(findnode(soln(0), :V))) +ax1 = Axis(fig[1,2], width = 400, height = 400) +ax2 = Axis(fig[1,4], width = 400, height = 400) +hidedecorations!(ax1) +hidedecorations!(ax2) +hidespines!(ax1) +hidespines!(ax2) +Colorbar(fig[1,1]) +Colorbar(fig[1,5]) +Label(fig[1,2,Top()], "U") +Label(fig[1,4,Top()], "V") +lab1 = Label(fig[1,3], "") + +# Animation +record(fig, "brusselator_bounded.gif", range(0.0, tₑ; length=frames); framerate = 15) do t + p1.plot.color = findnode(soln(t), :U) + p2.plot.color = findnode(soln(t), :V) + lab1.text = @sprintf("%.2f",t) +end + +end # END Gif creation diff --git a/src/Decapodes.jl b/src/Decapodes.jl index 62f09057..7d7e9cee 100644 --- a/src/Decapodes.jl +++ b/src/Decapodes.jl @@ -22,7 +22,10 @@ export normalize_unicode, DerivOp, append_dot, unicode!, vec_to_dec!, AbstractMeshKey, loadmesh, Icosphere, Rectangle_30x10, Torus_30x10, Point_Map, Open, OpenSummationDecapodeOb, OpenSummationDecapode, unique_by, unique_by!, oapply, CartesianPoint, SpherePoint, r, theta, phi, TangentBasis, θhat, ϕhat, - #average_rewrite, + average_rewrite, recursive_delete_parents, contract_operators, + default_dec_matrix_generate, default_dec_generate, + AbstractDecapodeMorphism, BCMorphism, ICMorphism, + AbstractCollage, Collage, simulation_helper, recursive_delete_parents, contract_operators, default_dec_matrix_generate, default_dec_generate, Term, @decapode @@ -35,6 +38,7 @@ append_dot(s::Symbol) = Symbol(string(s)*'\U0307') include("decapodeacset.jl") include("language.jl") include("composition.jl") +include("collages.jl") include("coordinates.jl") include("visualization.jl") include("simulation.jl") diff --git a/src/collages.jl b/src/collages.jl new file mode 100644 index 00000000..e152b0e9 --- /dev/null +++ b/src/collages.jl @@ -0,0 +1,180 @@ + +abstract type AbstractDecapodeMorphism end + +struct BCMorphism <: AbstractDecapodeMorphism + morphism::ACSetTransformation +end + +struct ICMorphism <: AbstractDecapodeMorphism + morphism::ACSetTransformation +end + +abstract type AbstractCollage end + +struct Collage <: AbstractCollage + bc::BCMorphism + ic::ICMorphism +end + +""" function transfer_sources!(d::SummationDecapode, from::Int, to::Int) + +Transfer any operations of which the Var `from` is the source of to `to`, excluding the special `∂ₜ` operator. +""" +function transfer_sources!(d::SummationDecapode, from::Int, to::Int) + # All Σ of which `from` is a summand: + summands = incident(d, from, :summand) + # All Op1s of which `from` is the source, excluding ∂ₜ: + op1s = filter(i -> d[i, :name] != :∂ₜ, incident(d, from, :src)) + # All Op2s of which `from` is the proj1: + proj1s = incident(d, from, :proj1) + proj2s = incident(d, from, :proj2) + + # Transfer summands: + # TODO: Change this fill idiom when ACSet assignment syntax changes. + d[summands, :summand] = fill(to, length(summands)) + # Transfer op1s: + d[op1s, :src] = fill(to, length(op1s)) + # Transfer proj1s: + d[proj1s, :proj1] = fill(to, length(proj1s)) + # Transfer proj2s: + d[proj2s, :proj2] = fill(to, length(proj2s)) + + d +end + +""" function transfer_targets!(d::SummationDecapode, from::Int, to::Int) + +Transfer any operations of which the Var `from` is the target of to `to`, excluding the special `∂ₜ` operator. +""" +function transfer_targets!(d::SummationDecapode, from::Int, to::Int) + # All Σ of which `from` is the target: + summations = incident(d, from, :sum) + # All Op1s of which `from` is the target, excluding ∂ₜ: + op1s = filter(i -> d[i, :name] != :∂ₜ, incident(d, from, :tgt)) + # All Op2s of which `from` is the target: + op2s = incident(d, from, :res) + + # Transfer sums: + # TODO: Change this fill idiom when ACSet assignment syntax changes. + d[summations, :sum] = fill(to, length(summations)) + # Transfer op1s: + d[op1s, :tgt] = fill(to, length(op1s)) + # Transfer op2s: + d[op2s, :res] = fill(to, length(op2s)) + + d +end + +""" function collate(dm::BCMorphism) + +"Compile" a collage of Decapodes to a simulatable one. +``` +""" +function collate(dm::BCMorphism) + dm = dm.morphism + d = SummationDecapode{Any, Any, Symbol}() + copy_parts!(d, dm.codom, (:Var, :TVar, :Op1, :Op2, :Σ, :Summand)) + + for (i,x) in enumerate(dm.components.Var.func) + op_name = Symbol("∂_mask") + mask_var = add_part!(d, :Var, type=dm.dom[i, :type], name=dm.dom[i, :name]) + tgt_name = dm.codom[x, :name] + tgt_idx = only(incident(d, tgt_name, :name)) + d[tgt_idx, :name] = Symbol(string(d[tgt_idx, :name]) * string(i)) + res_var = add_part!(d, :Var, type=dm.codom[x, :type], name=tgt_name) + is_tvar = !isempty(incident(d, tgt_idx, :incl)) + if is_tvar + transfer_sources!(d, tgt_idx, res_var) + add_part!(d, :Op2, proj1=tgt_idx, proj2=mask_var, res=res_var, op2=op_name) + else + transfer_targets!(d, tgt_idx, res_var) + add_part!(d, :Op2, proj1=res_var, proj2=mask_var, res=tgt_idx, op2=op_name) + end + tan_dir = is_tvar ? :tgt : :src + tangent_op1s = filter(x -> d[x, :op1]==:∂ₜ, incident(d, tgt_idx, tan_dir)) + isempty(tangent_op1s) && continue + d[only(tangent_op1s), tan_dir] = res_var + d[incident(d, tgt_idx, :incl), :incl] = res_var + end + d +end + +""" function make_bc_loader(dm::Collage, dimension) + +Given a collage, return a function that accepts a mapping of Vars to masking functions. + +This returned function will take in a mesh and non-boundary constants-and-parameters, and return a named tuple containing the parameters as evaluated on the mesh, and the regular parameters. (This final named tuple is suitable to be passed to an ODEProblem as `p`.) +""" +function make_bc_loader(dm::Collage, dimension) + function loader(mask_funcs::Dict{Symbol, N}) where {N <: Function} + function generator(sd, cs_ps::NamedTuple) + vars = keys(mask_funcs) + for (_, bc_var) in enumerate(dm.bc.morphism.dom[:name]) + bc_var ∉ vars && error("BC Variable $(string(bc_var)) is not given a generating function.") + # TODO: Also perform error checking on the return type of the + # corresponding mask_func. i.e. Base.return_types + # Parameter return types should be <: Function, and Constants should be + # tuples. Base.return_types behavior is not guaranteed to not return + # Any, though. + end + mask_pairs = map(values(mask_funcs)) do func + func(sd) + end + merge( + NamedTuple{Tuple(vars)}(mask_pairs), + cs_ps) + end + end + loader +end + +""" function make_ic_loader(dm::Collage, dimension) + +Given a collage, return a function that accepts a mapping of Vars to initial conditions. + +This returned function will take in a mesh, and return a MultScaleArray as evaluated on the mesh. (This final MultiScaleArray is suitable to be passed to an ODEProblem as `u`.) +""" +function make_ic_loader(dm::Collage, dimension) + function loader(ic_funcs::Dict{Symbol, N}) where {N <: Function} + function generator(sd) + vars = keys(ic_funcs) + for ic_var in dm.ic.morphism.dom[:name] + ic_var ∉ vars && error("IC Variable $(string(ic_var)) is not given a generating function.") + end + ics = map(values(ic_funcs)) do func + func(sd) + end + codom_names = map(collect(vars)) do var + var_idx = incident(dm.ic.morphism.dom, var, :name) + only(dm.ic.morphism.codom[dm.ic.morphism.components.Var.func[var_idx], :name]) + end + for (var, ic) in zip(vars,ics) + var_idx = incident(dm.ic.morphism.dom, var, :name) + type = only(dm.ic.morphism.dom[var_idx, :type]) + simplex = Decapodes.form_simplex(type, dimension) + if simplex != :AllocVecCall_Error && + (dm.ic.morphism.dom[var_idx, :type] ∉ [:Constant, :Parameter, :infer]) && + length(ic) != nparts(sd, simplex) + # TODO: This error-catch may not work if say, the number of edges is + # equal to the number of vertices, (i.e. in a circle). + error("IC Variable $(string(var)) was declared to be a $(string(type)), but is of length $(length(ic)), not $(nparts(sd, simplex)).") + end + end + construct(PhysicsState, + VectorForm.(ics), + Float64[], + collect(codom_names)) + end + end + loader +end + +""" function simulation_helper(dm::Collage, dimension=2) + +Given a collage, return functions to load boundary conditions, initial conditions, and a simulation. +``` +""" +simulation_helper(dm::Collage; dimension=2) = ( + make_bc_loader(dm, dimension), + make_ic_loader(dm, dimension), + gensim( dm.bc, dimension=dimension)) diff --git a/src/simulation.jl b/src/simulation.jl index 79e19fd0..0b7cb576 100644 --- a/src/simulation.jl +++ b/src/simulation.jl @@ -71,23 +71,25 @@ struct AllocVecCallError <: Exception c::AllocVecCall end +form_simplex(form, dim) = @match (form, dim) begin + (:Form0, 2) => :V + (:Form1, 2) => :E + (:Form2, 2) => :Tri + (:DualForm0, 2) => :Tri + (:DualForm1, 2) => :E + (:DualForm2, 2) => :V + + (:Form0, 1) => :V + (:Form1, 1) => :E + (:DualForm0, 1) => :E + (:DualForm1, 1) => :V + _ => return :AllocVecCall_Error +end + # TODO: There are likely better ways of dispatching on dimension instead of # storing it inside an AllocVecCall. Base.Expr(c::AllocVecCall) = begin - resolved_form = @match (c.form, c.dimension) begin - (:Form0, 2) => :V - (:Form1, 2) => :E - (:Form2, 2) => :Tri - (:DualForm0, 2) => :Tri - (:DualForm1, 2) => :E - (:DualForm2, 2) => :V - - (:Form0, 1) => :V - (:Form1, 1) => :E - (:DualForm0, 1) => :E - (:DualForm1, 1) => :V - _ => throw(AllocVecCallError(c)) - end + resolved_form = form_simplex(c.form, c.dimension) :($(Symbol(:__,c.name)) = Decapodes.FixedSizeDiffCache(Vector{$(c.T)}(undef, nparts(mesh, $(QuoteNode(resolved_form)))))) end @@ -436,6 +438,9 @@ function gensim(user_d::AbstractNamedDecapode, input_vars; dimension::Int=2, end end +gensim(dm::BCMorphism; dimension::Int=2) = + gensim(collate(dm); 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`. @@ -535,6 +540,12 @@ function default_dec_generate(sd, my_symbol, hodge=GeometricHodge()) # Lie Derivative 0 :L₀ => dec_lie_derivative_zero(sd, hodge) + # Mask application + :∂_mask => (x,y) -> begin + x[y[1]] .= y[2] + x + end + x => error("Unmatched operator $my_symbol") end diff --git a/test/collages.jl b/test/collages.jl new file mode 100644 index 00000000..783e4ca1 --- /dev/null +++ b/test/collages.jl @@ -0,0 +1,324 @@ +using Test +using Decapodes +using Catlab +using CombinatorialSpaces +using GeometryBasics: Point2 +Point2D = Point2{Float64} +using MultiScaleArrays + +# TODO: General boundaries i.e. arbitrary functions of solutions. (More complex +# relationships between what the boundary "is" and how to interpret state to +# provide values there.) + +# Test `transfer_sources!` transfers summands, op1s, proj1s, and proj2s. +Test1 = @decapode begin + To::Form0 + A == B+From+D + X == f(From) + Y == g(From,G) + Z == h(F,From) +end +Decapodes.transfer_sources!(Test1, + only(incident(Test1, :From, :name)), only(incident(Test1, :To, :name))) + +@test Test1 == @acset SummationDecapode{Any, Any, Symbol} begin + Var = 10 + Op1 = 1 + Op2 = 2 + Σ = 1 + Summand = 3 + src = [1] + tgt = [7] + proj1 = [1, 10] + proj2 = [6, 1] + res = [9, 8] + incl = Int64[] + summand = [3, 1, 5] + summation = [1, 1, 1] + sum = [2] + op1 = [:f] + op2 = [:g, :h] + type = [:Form0, :infer, :infer, :infer, :infer, :infer, :infer, :infer, :infer, :infer] + name = [:To, :A, :B, :From, :D, :G, :X, :Z, :Y, :F] +end + +# Test `transfer_targets!` transfers sums, op1s, and op2s. +Test2 = @decapode begin + To::Form0 + From == B+C+D + From == f(E) + From == g(F,G) +end +Decapodes.transfer_targets!(Test2, + only(incident(Test2, :From, :name)), only(incident(Test2, :To, :name))) + +@test Test2 == @acset SummationDecapode{Any, Any, Symbol} begin + Var = 8 + Op1 = 1 + Op2 = 1 + Σ = 1 + Summand = 3 + src = [8] + tgt = [1] + proj1 = [7] + proj2 = [6] + res = [1] + incl = Int64[] + summand = [3, 4, 5] + summation = [1, 1, 1] + sum = [1] + op1 = [:f] + op2 = [:g] + type = [:Form0, :infer, :infer, :infer, :infer, :infer, :infer, :infer] + name = [:To, :From, :B, :C, :D, :G, :F, :E] +end + +# Test simple boundary masks. +DiffusionDynamics = @decapode begin + K::Form0 + ∂ₜ(K) == ∘(d,⋆,d,⋆)(K) +end +DiffusionBoundaries = @decapode begin + (Kb1, Kb2)::Constant + Null::Parameter +end + +# Test that simple boundary masks work on state variables. +StateMorphism = BCMorphism(ACSetTransformation( + DiffusionBoundaries, DiffusionDynamics, + Var = [1,1,1])) + +Diffusion = Decapodes.collate(StateMorphism) + +@test Diffusion == @acset SummationDecapode{Any, Any, Symbol} begin + Var = 8 + TVar = 1 + Op1 = 2 + Op2 = 3 + src = [8, 1] + tgt = [2, 2] + proj1 = [4, 6, 8] + proj2 = [3, 5, 7] + res = [1, 4, 6] + incl = [2] + op1 = Any[:∂ₜ, [:d, :⋆, :d, :⋆]] + op2 = [:∂_mask, :∂_mask, :∂_mask] + type = [:Form0, :infer, :Constant, :Form0, :Constant, :Form0, :Parameter, :Form0] + name = [:K1, :K̇, :Kb1, :K2, :Kb2, :K3, :Null, :K] +end + +# Test that simple boundary masks work on state and tangent variables. +StateTangentMorphism = BCMorphism(ACSetTransformation( + DiffusionBoundaries, DiffusionDynamics, + Var = [1,1,2])) + +Diffusion = Decapodes.collate(StateTangentMorphism) + +@test Diffusion == @acset SummationDecapode{Any, Any, Symbol} begin + Var = 8 + TVar = 1 + Op1 = 2 + Op2 = 3 + src = [6, 1] + tgt = [8, 2] + proj1 = [4, 6, 2] + proj2 = [3, 5, 7] + res = [1, 4, 8] + incl = [8] + op1 = Any[:∂ₜ, [:d, :⋆, :d, :⋆]] + op2 = [:∂_mask, :∂_mask, :∂_mask] + type = [:Form0, :infer, :Constant, :Form0, :Constant, :Form0, :Parameter, :infer] + name = [:K1, :K̇3, :Kb1, :K2, :Kb2, :K, :Null, :K̇] +end + +# Test that simple boundary masks work on second order tangent variables. +SecondDerivatives = @decapode begin + (X,Y)::Form0 + ∂ₜ(X) == Z + Z == Δ(X) + ∂ₜ(Y) == Δ(Z) +end +SingleBC = @decapode begin + M::Form0 +end +to_graphviz(SecondDerivatives) +TangentTangentMorphism = BCMorphism(ACSetTransformation( + SingleBC, SecondDerivatives, + Var = [4])) + +SecondDerivativesBounded = Decapodes.collate(TangentTangentMorphism) + +@test SecondDerivativesBounded == @acset SummationDecapode{Any, Any, Symbol} begin + Var = 6 + TVar = 2 + Op1 = 4 + Op2 = 1 + src = [1, 1, 2, 6] + tgt = [6, 4, 3, 3] + proj1 = [4] + proj2 = [5] + res = [6] + incl = [6, 3] + summand = Int64[] + summation = Int64[] + sum = Int64[] + op1 = [:∂ₜ, :Δ, :∂ₜ, :Δ] + op2 = [:∂_mask] + type = [:Form0, :Form0, :infer, :infer, :Form0, :infer] + name = [:X, :Y, :Ẏ, :Z1, :M, :Z] +end + +# Test gensim on a collage. +@test gensim(StateTangentMorphism) == gensim(Diffusion) + +# Test that simple boundary masks work on intermediate variables. +DiffusionExpanded = expand_operators(DiffusionDynamics) +IntermediateMorphism = BCMorphism(ACSetTransformation( + DiffusionBoundaries, DiffusionExpanded, + Var = [3,3,3])) + +Diffusion = Decapodes.collate(IntermediateMorphism) +@test Diffusion == @acset SummationDecapode{Any, Any, Symbol} begin + Var = 11 + TVar = 1 + Op1 = 5 + Op2 = 3 + src = [1, 1, 3, 4, 5] + tgt = [2, 11, 4, 5, 2] + proj1 = [7, 9, 11] + proj2 = [6, 8, 10] + res = [3, 7, 9] + incl = [2] + op1 = [:∂ₜ, :d, :⋆, :d, :⋆] + op2 = [:∂_mask, :∂_mask, :∂_mask] + type = [:Form0, :infer, :infer, :infer, :infer, :Constant, :infer, :Constant, :infer, :Parameter, :infer] + name = [:K, :K̇, Symbol("•_2_11"), Symbol("•_2_2"), +Symbol("•_2_3"), :Kb1, Symbol("•_2_12"), :Kb2, Symbol("•_2_13"), :Null, Symbol("•_2_1")] +end + +# Test collate with empty boundaries. +EmptyBoundaries = @decapode begin +end + +EmptyMorphism = BCMorphism(ACSetTransformation( + EmptyBoundaries, DiffusionDynamics)) + +Diffusion = Decapodes.collate(EmptyMorphism) +@test Diffusion == DiffusionDynamics + +# Test initial conditions via collage. +DiffusionICs = @decapode begin + (Kic)::Form0 +end +SingleICMorphism = ICMorphism(ACSetTransformation( + DiffusionICs, DiffusionDynamics, + Var = [1])) + +# Test the collage workflow with IC and BC morphisms, step-by-step. +DiffusionCollage = Collage( + BCMorphism(ACSetTransformation( + DiffusionBoundaries, DiffusionDynamics, + Var = [1,1,2])), + ICMorphism(ACSetTransformation( + DiffusionICs, DiffusionDynamics, + Var = [1]))) + +# We will use this mesh to test that Boundary and Initial Condition generators +# generate appropriate data structures. +# The fifth subdivision is chosen so that numeric results of applying the +# boundary condition computations can be compared with analytic methods up to +# some number of digits. +s = loadmesh(Icosphere(5)) +sd = EmbeddedDeltaDualComplex2D{Bool,Float64,Point2D}(s) +subdivide_duals!(sd, Barycenter()) + +# Boundary Condition: +# Set all points with y-coordinate above cos(t) to sin(t). +# e.g. The polar ice cap is a source term, which expands and contracts +# seasonally, and so does its albedo. +time_varying_bc_function(sd) = + t -> (map(x -> x[2] > cos(t), sd[:point]), + fill(sin(t), count(x -> x[2] > 0.9, sd[:point]))) +# Boundary Condition: +# Set all points with y-coordinate below -0.9 to 0.0. +constant_bc_function(sd) = ( + map(x -> x[2] > 0.9, sd[:point]), + fill(0.0, count(x -> x[2] > 0.9, sd[:point]))) + +bc_loader = Decapodes.make_bc_loader(DiffusionCollage, 2) + +bc_generator = bc_loader(Dict( + :Kb1 => constant_bc_function, + :Kb2 => constant_bc_function, + :Null => time_varying_bc_function)) + +# Test that the BC generator returns a tuple of the correct type. +p = bc_generator(sd, (a = 1.0, b = 2.0)) +@test typeof(p) <: + NamedTuple{(:Kb2, :Kb1, :Null, :a, :b), Tuple{ + # Constant: a bit mask and a vector of floats + Tuple{Vector{Bool}, Vector{Float64}}, + Tuple{Vector{Bool}, Vector{Float64}}, + # Parameter: a function that returns '' + N, + # "Regular" constants and parameters: which happen to be floats here + Float64, Float64}} where {N <: Function} + +# Test that the p functions are passed and evaluate as expected. +@test typeof(p.Null(π/4)) == Tuple{Vector{Bool}, Vector{Float64}} +# Some basic trigonometry to check the correct number of vertices are selected. +@test count(p.Null(π/4)[1]) / nv(sd) - + abs((1/2)*(sqrt(1- 1^2)*1 + asin(1 )) - + (1/2)*(sqrt(1-cos(π/4)^2)*cos(π/4) + asin(cos(π/4)))) < 1e-3 +@test all(p.Null(π/4)[2] .== sin(π/4)) + +# Test that p does not "blow up" memory by making a copy of the mesh. +# This 1% mark is a heuristic. +@test sizeof(p) / Base.summarysize(p) < 0.01 +@test sizeof(p) / Base.summarysize(sd) < 0.01 + +# Test that the boundary condition generator will catch if initial conditions +# are not provided. +bc_generator = bc_loader(Dict( + :Kb1 => constant_bc_function, + #:Kb2 => constant_bc_function, + :Null => time_varying_bc_function)) +# Test that the BC generator returns a tuple of the correct type. +@test_throws "BC Variable Kb2 is not given a generating function." bc_generator(sd, (a = 1.0, b = 2.0)) + +# Initial Condition: +# Assign each vertex the sin of their y-coordinate. +ic_function(sd) = map(x -> sin(x[2]), sd[:point]) +# This function will return a vector of length the number of edges in the mesh, +# and we will erroneously assign it to a Form0 so as to test error-checking. +false_ic_function(sd) = ones(ne(sd)) + +ic_loader = Decapodes.make_ic_loader(DiffusionCollage, 2) +ic_generator = ic_loader(Dict(:Kic => ic_function)) +u = ic_generator(sd) + +# Test that the initial condition functions are passed and evaluate as expected. +@test typeof(u) == PhysicsState{VectorForm{Float64}, Float64} +@test findnode(u, :K) == ic_function(sd) +# Test that the generator is equivalent to setting initial conditions manually. +# Note: Testing equality with `==` fails for MultiScaleArrays. +# i.e. This test will fail: +#@test construct(PhysicsState,[VectorForm(ic_function(sd))], Float64[], [:K]) == +# construct(PhysicsState,[VectorForm(ic_function(sd))], Float64[], [:K]) +# So, we check for equality between the accessed components. +v = construct(PhysicsState, [VectorForm(ic_function(sd))], Float64[], [:K]) +@test all(findnode(u, :K) .- findnode(v, :K) .== 0) + +# Test the error-checking on the length of initial condition data. +ic_generator = ic_loader(Dict(:Kic => false_ic_function)) +@test_throws "IC Variable Kic was declared to be a Form0, but is of length 7680, not 2562." ic_generator(sd) + +# Test the `simulation_helper` function, which wraps the boundary conditions +# loader, initial conditions loader, and the simulation generator. +sim = gensim(Decapodes.collate(DiffusionCollage.bc)) + +auto_bc_loader, auto_ic_loader, auto_sim = + simulation_helper(DiffusionCollage, dimension=2) +@test auto_bc_loader == bc_loader +@test auto_ic_loader == ic_loader +@test auto_sim == sim diff --git a/test/runtests.jl b/test/runtests.jl index 8261dc46..e386fa9d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -14,6 +14,10 @@ end include("composition.jl") end +@testset "Collages" begin + include("collages.jl") +end + @testset "Coordinates" begin include("coordinates.jl") end