diff --git a/packages/algjulia-service/Project.toml b/packages/algjulia-service/Project.toml index f991b4fe..08fe0659 100644 --- a/packages/algjulia-service/Project.toml +++ b/packages/algjulia-service/Project.toml @@ -8,6 +8,7 @@ version = "0.1.0" ACSets = "227ef7b5-1206-438b-ac65-934d6da304b8" CombinatorialSpaces = "b1c52339-7909-45ad-8b6a-6e388f7c67f2" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +CoordRefSystems = "b46f11dc-f210-4604-bfba-323c1ec968cb" Decapodes = "679ab3ea-c928-4fe6-8d59-fd451142d391" DiagrammaticEquations = "6f00c28b-6bed-4403-80fa-30e0dc12f317" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" @@ -17,11 +18,13 @@ JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ACSets = "0.2.21" ComponentArrays = "0.15" +CoordRefSystems = "0.16.0" Decapodes = "0.5.6" DiagrammaticEquations = "0.1.7" Distributions = "0.25" @@ -30,4 +33,5 @@ JSON3 = "1" LinearAlgebra = "1" MLStyle = "0.4" OrdinaryDiffEq = "6" +Reexport = "1.2.2" StaticArrays = "1" diff --git a/packages/algjulia-service/src/AlgebraicJuliaService.jl b/packages/algjulia-service/src/AlgebraicJuliaService.jl index 7a2b222c..a15f9938 100644 --- a/packages/algjulia-service/src/AlgebraicJuliaService.jl +++ b/packages/algjulia-service/src/AlgebraicJuliaService.jl @@ -1,6 +1,36 @@ module AlgebraicJuliaService +using MLStyle +using Reexport + +# this code tracks integrations and allows for basic theory/model-building code to dispatch from it. +# the intent is that this is an interface for AlgebraicJulia code to interoperate with CatColab +@data AlgebraicJuliaIntegration begin + ThDecapode() +end +export ThDecapode + +# these are just objects for dispatching `to_theory` +@data ElementType begin + ObType() + HomType() +end +export ObType, HomType + +""" Functions to build a dictionary associating ids in the theory to elements in the model""" +function to_theory end; export to_theory + +struct ImplError <: Exception + name::String +end +export ImplError + +Base.showerror(io::IO, e::ImplError) = print(io, "$(e.name) not implemented") + + include("kernel_support.jl") -include("decapodes.jl") +include("decapodes-service/DecapodesService.jl") + +@reexport using .DecapodesService end diff --git a/packages/algjulia-service/src/decapodes-service/DecapodesService.jl b/packages/algjulia-service/src/decapodes-service/DecapodesService.jl new file mode 100644 index 00000000..5e9712b8 --- /dev/null +++ b/packages/algjulia-service/src/decapodes-service/DecapodesService.jl @@ -0,0 +1,180 @@ +module DecapodesService + +# algebraicjulia dependencies +using ACSets +using DiagrammaticEquations +using Decapodes +using Decapodes: dec_mat_dual_differential, dec_mat_inverse_hodge +using CombinatorialSpaces + +# dependencies +import JSON3 +using StaticArrays +using MLStyle +using LinearAlgebra +using ComponentArrays +using Distributions # for initial conditions + +# meshing +using CoordRefSystems +using GeometryBasics: Point2, Point3 +Point3D = Point3{Float64}; + +# simulation +using OrdinaryDiffEq + +using ..AlgebraicJuliaService +import ..AlgebraicJuliaService: to_theory + +# necessary to export +export infer_types!, evalsim, default_dec_generate, default_dec_matrix_generate, + DiagonalHodge, ComponentArray + +# funcitons for geometry and initial conditions +include("geometry.jl") +include("initial_conditions.jl") +include("theory.jl") ## theory-building +include("model.jl") ## model-building (is this actually a diagram of a model?) + +struct PodeSystem + pode::SummationDecapode + plotvar::Vector{Symbol} + scalars::Dict{Symbol, Any} # closures + geometry::Geometry + init::ComponentArray # TODO Is it always ComponentVector? + generate::Any + uuiddict::Dict{Symbol, String} +end +export PodeSystem + +function PodeSystem(json_string::String, args...) + json_object = JSON3.read(json_string) + PodeSystem(json_object, args...) +end + +""" +Construct a `PodeSystem` object from a JSON string. +""" +function PodeSystem(json_object::AbstractDict, hodge=GeometricHodge()) + # make a theory of the DEC, valued in Julia + theory = Theory(json_object[:model]) + + # this is a diagram in the model of the DEC. it wants to be a decapode! + diagram = json_object[:diagram] + + # any scalars? + scalars = haskey(json_object, :scalars) ? json_object[:scalars] : [] + + # pode, anons, and vars (UUID => ACSetId) + decapode, anons, vars = Decapode(diagram, theory; scalars=scalars) + dot_rename!(decapode) + uuid2symb = uuid_to_symb(decapode, vars) + + # plotting variables + plotvars = [uuid2symb[uuid] for uuid in json_object[:plotVariables]] + + # extract the domain in order to create the mesh, dualmesh + geometry = Geometry(json_object) + + # initialize operators + ♭♯_m = ♭♯_mat(geometry.dualmesh) + wedge_dp10 = dec_wedge_product_dp(Tuple{1,0}, geometry.dualmesh) + dual_d1_m = dec_mat_dual_differential(1, geometry.dualmesh) + star0_inv_m = dec_mat_inverse_hodge(0, geometry.dualmesh, hodge) + Δ0 = Δ(0,geometry.dualmesh) + #fΔ0 = factorize(Δ0); + function sys_generate(s, my_symbol, hodge=hodge) + op = @match my_symbol begin + sym && if sym ∈ keys(anons) end => anons[sym] + :♭♯ => x -> ♭♯_m * x # [1] + :dpsw => x -> wedge_dp10(x, star0_inv_m[1]*(dual_d1_m[1]*x)) + :Δ⁻¹ => x -> begin + y = Δ0 \ x + y .- minimum(y) + end + _ => default_dec_matrix_generate(s, my_symbol, hodge) + end + return (args...) -> op(args...) + end + # end initialize + + # initial conditions + u0 = initial_conditions(json_object, geometry, uuid2symb) + + # symbol => uuid. we need this to reassociate the var + symb2uuid = Dict([v => k for (k,v) in pairs(uuid2symb)]) + + return PodeSystem(decapode, plotvars, anons, geometry, u0, sys_generate, symb2uuid) +end +export PodeSystem + +points(system::PodeSystem) = collect(values(system.geometry.mesh.subparts.point.m)) +indexing_bounds(system::PodeSystem) = indexing_bounds(system.geometry.domain) + +## SIM HELPERS ## + +function run_sim(fm, u0, t0, constparam) + prob = ODEProblem(fm, u0, (0, t0), constparam) + soln = solve(prob, Tsit5(), saveat=0.01) +end +export run_sim + +struct SimResult + time::Vector{Float64} + state::Dict{String, Vector{AbstractArray{SVector{3, Float64}}}} + x::Vector{Float64} # axis + y::Vector{Float64} +end +export SimResult + +function SimResult(soln::ODESolution, system::PodeSystem) + idx_bounds = indexing_bounds(system) + state_val_dict = variables_state(soln, system) # Dict("UUID1" => VectorMatrixSVectr...) + SimResult(soln.t, state_val_dict, 0:idx_bounds.x, 0:idx_bounds.y) +end +# TODO generalize to HasDeltaSet + +""" for the variables in a system, associate them to their state values over the duration of the simulation """ +function variables_state(soln::ODESolution, system::PodeSystem) + Dict([ system.uuiddict[var] => state_entire_sim(soln, system, var) for var ∈ system.plotvar ]) +end + +""" given a simulation, a domain, and a variable, gets the state values over the duration of a simulation. +Called by `variables_state`[@ref] """ +function state_entire_sim(soln::ODESolution, system::PodeSystem, var::Symbol) + map(1:length(soln.t)) do i + state_at_time(soln, system, var, i) + end +end + +# TODO type `points` +function state_at_time(soln::ODESolution, system::PodeSystem, plotvar::Symbol, t::Int) + @match system.geometry.domain begin + # TODO check time indexing here + domain::Rectangle => state_at_time(soln, domain, plotvar, t) + domain::Sphere => state_at_time(soln, domain, plotvar, t, points(system)) + _ => throw(ImplError("state_at_time function for domain $domain")) + end +end + +function state_at_time(soln::ODESolution, domain::Rectangle, var::Symbol, t::Int) # last two args can be one Point2 + (x, y) = indexing_bounds(domain) + [SVector(i, j, getproperty(soln.u[t], var)[(x+1)*(i-1) + j]) for i in 1:x+1, j in 1:y+1] +end + +# TODO just separated this from the SimResult function and added type parameters, but need to generalize +function grid(pt3::Point3, grid_size::Vector{Int}) + pt2 = [(pt3[1]+1)/2, (pt3[2]+1)/2] + [round(Int, pt2[1]*grid_size[1]), round(Int, pt2[2]*grid_size[2])] +end + +function state_at_time(soln::ODESolution, domain::Sphere, var::Symbol, t::Int, points) + l , _ = indexing_bounds(domain) # TODO this is hardcoded to return 100, 100 + northern_indices = filter(i -> points[i][3] > 0, keys(points)) + map(northern_indices) do n + i, j = grid(points[n], [l, l]) # TODO + SVector(i, j, getproperty(soln.u[t], var)[n]) + end +end + +end diff --git a/packages/algjulia-service/src/decapodes-service/geometry.jl b/packages/algjulia-service/src/decapodes-service/geometry.jl new file mode 100644 index 00000000..329abdec --- /dev/null +++ b/packages/algjulia-service/src/decapodes-service/geometry.jl @@ -0,0 +1,117 @@ +## INTEROP + +""" Supported domains. """ +const domain_names = [:Plane, :Sphere] + +""" Mapping from supported domains to meshes for the domain. """ +const mesh_names = Dict( + :Plane => [:Rectangle, :Periodic], + :Sphere => [:Icosphere6, :Icosphere7, :Icosphere8, :UV], +) + +""" Mapping from supported domains to initial/boundary conditions. """ +const initial_condition_names = Dict( + :Plane => [:Gaussian], + :Sphere => [:TaylorVortex, :SixVortex], +) + +""" Supported geometries, in the JSON format expected by the frontend. """ +function supported_decapodes_geometries() + domains = map(domain_names) do domain + Dict( + :name => domain, + :meshes => mesh_names[domain], + :initialConditions => initial_condition_names[domain], + ) + end + Dict(:domains => domains) +end +export supported_decapodes_geometries + +## DOMAINS + +abstract type Domain end + +# meshes associated with Planes +@data Planar <: Domain begin + Rectangle(max_x::Int, max_y::Int, dx::Float64, dy::Float64) + Periodic # TODO +end + +# rectangle methods + +# TODO it is semantically better to case to Point2? +middle(r::Rectangle) = [r.max_x/2, r.max_y/2] + +function indexing_bounds(r::Rectangle) + (x=floor(Int, r.max_x/r.dx), y=floor(Int, r.max_y/r.dy)) +end + +# meshes associated with Spheres +@data Spheric <: Domain begin + Sphere(dim::Int, radius::Float64) + UV(minlat::Int, maxlat::Int, dlat::Float64, minlong::Int, maxlong::Int, dlong::Float64, radius::Float64) +end + +# default +Sphere(dim) = Sphere(dim, 1.0) + +# TODO XXX hardcoded alert! +function indexing_bounds(m::Sphere) + (x=100, y=100) +end + +""" helper function for UV """ +function makeSphere(m::UV) + makeSphere(m.minlat, m.maxlat, m.dlat, m.minlong, m.maxlong, m.dlong, m.radius) +end + +## GEOMETRY + +struct Geometry + domain::Domain + mesh::HasDeltaSet + dualmesh::HasDeltaSet +end + +function Geometry(json_object::AbstractDict) + mesh_name = Symbol(json_object[:mesh]) + domain = PREDEFINED_MESHES[mesh_name] + Geometry(domain) +end + +# function Geometry(d::Domain, args...) +# throw(ImplError("The mesh ($(d)) is")) +# end + +function Geometry(r::Rectangle, division::SimplexCenter=Circumcenter()) + s = triangulated_grid(r.max_x, r.max_y, r.dx, r.dy, Point2{Float64}) + sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point2{Float64}}(s) + subdivide_duals!(sd, division) + Geometry(r, s, sd) +end + +# function Geometry(r::Periodic, division::SimplexCenter=Circumcenter()) end + +function Geometry(m::Sphere, division::SimplexCenter=Circumcenter()) + s = loadmesh(Icosphere(m.dim, m.radius)) + sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3{Float64}}(s) + subdivide_duals!(sd, division) + Geometry(m, s, sd) +end + +function Geometry(m::UV, division::SimplexCenter=Circumcenter()) + s, _, _ = makeSphere(m) + sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3{Float64}}(s) + subdivide_duals!(sd, division) + Geometry(m, s, sd) +end + +## Prefined meshes + +const PREDEFINED_MESHES = Dict( + :Rectangle => Rectangle(100, 100, 2, 2), + :Icosphere6 => Sphere(6, 1.0), + :Icosphere7 => Sphere(7, 1.0), + :Icosphere8 => Sphere(8, 1.0), + :UV => UV(0, 180, 2.5, 0, 360, 2.5, 1.0)) diff --git a/packages/algjulia-service/src/decapodes-service/initial_conditions.jl b/packages/algjulia-service/src/decapodes-service/initial_conditions.jl new file mode 100644 index 00000000..40e8ee7a --- /dev/null +++ b/packages/algjulia-service/src/decapodes-service/initial_conditions.jl @@ -0,0 +1,109 @@ +## INITIAL CONDITIONS + +# TAYLOR VORTEX CODE +include("ns_helper.jl") +#### + +# This ADT defines the parameters for initial conditions data. +@data InitialConditionsData begin + GaussianData(μ::Vector{Float64}, Σ::Diagonal{Float64, Vector{Float64}}) + TaylorVortexData(lat::Float64, vortices::Int, p::AbstractVortexParams) +end + +function GaussianData(μ::Vector{Float64}, Σ::Vector{Float64}) + GaussianData(μ, LinearAlgebra.Diagonal(abs.(Σ))) +end + +# default method +function GaussianData(r::Rectangle) + μ = middle(r) + GaussianData(μ, μ/10) +end + +""" Normal distribution should understand GaussianData """ +Distributions.MvNormal(ξ::GaussianData) = MvNormal(ξ.μ, ξ.Σ) + +TaylorVortexData() = TaylorVortexData(0.2, 2, TaylorVortexParams(0.5, 0.1)) + +#= +This IC contains the domain and the initial conditions data. + +While these are currently tightly-interlinked with InitialConditionsData, they are formally separated to distinguish between the initial conditions schema and the data it might be parameterized over. +=# +@data InitialConditions begin + # planar + GaussianIC(r::Rectangle, ξ::GaussianData) + # spherical + TaylorVortexIC(d::Sphere, ξ::TaylorVortexData) + SixVortexIC(m::Sphere, data::Any) +end + +# DEFAULT METHOD +GaussianIC(r::Rectangle) = GaussianIC(r, GaussianData(r)) +TaylorVortexIC(d::Sphere) = TaylorVortexIC(d, TaylorVortexData()) + +function initial_conditions(json_object::AbstractDict, geometry::Geometry, uuid2symb::Dict{String, Symbol}) + ic_specs = json_object[:initialConditions] # this is "C" + dict = Dict([uuid2symb[string(uuid)] => ic_specs[string(uuid)] for uuid ∈ keys(ic_specs)]...) + initial_conditions(dict, geometry) # the resulting sim will only have (C,) as initial conditions +end + +""" Takes a string, a domain, and a mesh and returns the initial conditios object associated to it. + +Example: +``` +associate("TaylorVortex", Sphere(6, 1.0), sd) == TaylorVortexIC(Sphere(6, 1.0), sd) +``` +""" +function associate(str::String, geometry::Geometry) + @match str begin + "Gaussian" => GaussianIC(geometry.domain) + "TaylorVortex" => TaylorVortexIC(geometry.domain) + _ => error("$str is not implemented") + end +end + +""" Methods for this function implement initial conditions for their given schema. There are also helper functions.""" +function initial_conditions end; export initial_conditions + +""" associates the values in a dictionary to their initial condition flags, and passes the output to initial_conditions +""" +function initial_conditions(ics::Dict{Symbol, String}, geometry::Geometry) + ic_dict = Dict([var => associate(ics[var], geometry) for var in keys(ics)]...) + # Now we have a mapping between variables and their initial condition specs. + initial_conditions(ic_dict, geometry) +end + +""" builds a mapping between symbols and their initial conditions """ +function initial_conditions(ics::Dict{Symbol,<:InitialConditions}, geometry::Geometry) + u0 = ComponentArray(; Dict([ + var => initial_conditions(ics[var], geometry) for var ∈ keys(ics) + ])...) + return u0 +end + +function initial_conditions(x::InitialConditions, args...) + throw(ImplError("These initial conditions ($(x)) are")) +end + +function initial_conditions(ics::GaussianIC, geometry::Geometry) + c_dist = MvNormal(ics.ξ) + c = [pdf(c_dist, [p[1], p[2]]) for p ∈ geometry.dualmesh[:point]] + return c +end + +function vort_ring(ics::TaylorVortexIC, geometry::Geometry) + vort_ring(ics.d, ics.ξ.lat, ics.ξ.vortices, ics.ξ.p, geometry.dualmesh, taylor_vortex) +end + +function initial_conditions(ics::TaylorVortexIC, geometry::Geometry) + # TODO prefer not to load `s0` here but che sara sara + s0 = dec_hodge_star(0, geometry.dualmesh, GeometricHodge()) + X = vort_ring(ics, geometry) + du = s0 * X + return du +end + +function initial_conditions(ics::SixVortexIC, geometry::Geometry) + X = vort_ring(0.4, 6, PointVortexParams(3.0, 0.15), point_vortex) +end diff --git a/packages/algjulia-service/src/decapodes-service/model.jl b/packages/algjulia-service/src/decapodes-service/model.jl new file mode 100644 index 00000000..17110939 --- /dev/null +++ b/packages/algjulia-service/src/decapodes-service/model.jl @@ -0,0 +1,89 @@ +## MODEL BUILDING + +function add_to_pode! end +export add_to_pode! + +function add_to_pode!(d::SummationDecapode, + vars::Dict{String, Int}, # mapping between UUID and ACSet ID + theory::Theory, + content::AbstractDict, + nc::Vector{Int}, + ::ObType) + theory_elem = theory.data[content[:over][:content]] # indexes the theory by UUID + # checks if the cell is an anonymous (intermediate) variable. + # if so, we increment the intermediate variable counter and make an intermediate variable name. + # otherwise we use the existing name of the given content. + name = if isempty(content[:name]) + nc[1] += 1 + Symbol("•$(nc[1])") + else + Symbol(content[:name]) + end + id = add_part!(d, :Var, name=name, type=nameof(theory_elem)) + push!(vars, content[:id] => id) + return d +end + +function Base.nameof(theory::Theory, content::AbstractDict) + Symbol(theory.data[content[:over][:content]].name) +end + +# TODO we are restricted to Op1 +function add_to_pode!(d::SummationDecapode, + vars::Dict{String, Int}, # mapping between UUID and ACSet ID + theory::Theory, + content::AbstractDict, + scalars::Any, + anons::Dict{Symbol, Any}, + ::HomType) + + dom = content[:dom][:content] + cod = content[:cod][:content] + # TODO we need a safe way to fail this + if haskey(vars, dom) && haskey(vars, cod) + # get the name of the Op1 and add it to the theory + op1 = nameof(theory, content) + add_part!(d, :Op1, src=vars[dom], tgt=vars[cod], op1=op1) + # we need to add an inclusion to the TVar table + if op1 == :∂ₜ + add_part!(d, :TVar, incl=vars[cod]) + end + # if the dom is anonymous, we treat it as a something which will receive x -> k * x. + # we store its value in another array + if !isempty(scalars) && haskey(scalars, Symbol(content[:over][:content])) + scalar = scalars[Symbol(content[:over][:content])] + push!(anons, op1 => x -> scalar * x) + end + end + d +end + +""" Decapode(diagram::AbstractVector{<:AbstractDict}, theory::Theory) => (::SummationDecapode, ::Dict{Symbol, Any}, ::Dict{String, Int}) + +This returns + 1. a Decapode + 2. a dictionary of symbols mapped to anonymous functions + 3. a dictionary of JSON UUIDs mapped to symbols +""" +function Decapode(diagram::AbstractVector{<:AbstractDict}, theory::Theory; scalars=[]) + # initiatize decapode and its mapping between UUIDs and ACSet IDs + pode = SummationDecapode(parse_decapode(quote end)) + vars = Dict{String, Int}() # UUID => ACSetID + nc = [0] # array is a mutable container + anons = Dict{Symbol, Any}() + # for each cell in the notebook, add it to the diagram + foreach(diagram) do cell + @match cell begin + # TODO merge nameless_count into vars + IsObject(content) => add_to_pode!(pode, vars, theory, content, nc, ObType()) + IsMorphism(content) => add_to_pode!(pode, vars, theory, content, scalars, anons, HomType()) + _ => throw(ImplError(cell[:content][:tag])) + end + end + return pode, anons, vars +end +export Decapode + +function uuid_to_symb(decapode::SummationDecapode, vars::Dict{String, Int}) + Dict([key => (subpart(decapode, vars[key], :name)) for key ∈ keys(vars)]) +end diff --git a/packages/algjulia-service/src/decapodes-service/ns_helper.jl b/packages/algjulia-service/src/decapodes-service/ns_helper.jl new file mode 100644 index 00000000..b74e579a --- /dev/null +++ b/packages/algjulia-service/src/decapodes-service/ns_helper.jl @@ -0,0 +1,87 @@ +### +#This code was lifted from the Navier-Stokes simulation `ns.jl` in the Decapodes docs page, originally authored by Luke Morris +### +abstract type AbstractVortexParams end + +struct TaylorVortexParams <: AbstractVortexParams + G::Real + a::Real +end + +struct PointVortexParams <: AbstractVortexParams + τ::Real + a::Real +end + +""" function ring_centers(lat, n) + +Find n equispaced points at the given latitude. +""" +function ring_centers(lat, n, radius=1.0) + ϕs = range(0.0, 2π; length=n+1)[1:n] + map(ϕs) do ϕ + v_sph = Spherical(radius, lat, ϕ) + v_crt = convert(Cartesian, v_sph) + Point3D(v_crt.x.val, v_crt.y.val, v_crt.z.val) + end +end + +""" function great_circle_dist(pnt,G,a,cntr) + +Compute the length of the shortest path along a sphere, given Cartesian coordinates. +""" +function great_circle_dist(radius::Float64, pnt1::Point3D, pnt2::Point3D) + radius * acos(dot(pnt1,pnt2)) +end + +function taylor_vortex(sd::HasDeltaSet, radius::Float64, cntr::Point3D, p::TaylorVortexParams) + map(x -> taylor_vortex(x, radius, cntr, p), point(sd)) +end + +""" function taylor_vortex(pnt::Point3D, cntr::Point3D, p::TaylorVortexParams) + +Compute the value of a Taylor vortex at the given point. +""" +function taylor_vortex(pnt::Point3D, radius::Float64, cntr::Point3D, p::TaylorVortexParams) + gcd = great_circle_dist(radius, pnt, cntr) + (p.G/p.a) * (2 - (gcd/p.a)^2) * exp(0.5 * (1 - (gcd/p.a)^2)) +end + +""" function vort_ring(lat, n_vorts, p::T, formula) where {T<:AbstractVortexParams} + +Compute vorticity as primal 0-forms for a ring of vortices. + +Specify the latitude, number of vortices, and a formula for computing vortex strength centered at a point. +""" +function vort_ring(d::Sphere, lat, n_vorts, p::T, sd, formula) where {T<:AbstractVortexParams} + sum(map(x -> formula(sd, d.radius, x, p), ring_centers(lat, n_vorts, d.radius))) +end + +""" function vort_ring(lat, n_vorts, p::PointVortexParams, formula) + +Compute vorticity as primal 0-forms for a ring of vortices. + +Specify the latitude, number of vortices, and a formula for computing vortex strength centered at a point. + +Additionally, place a counter-balance vortex at the South Pole such that the integral of vorticity is 0. +""" +function vort_ring(radius, lat, n_vorts, p::PointVortexParams, formula) + Xs = sum(map(x -> formula(radius, sd, x, p), ring_centers(lat, n_vorts))) + Xsp = point_vortex(sd, Point3D(0.0, 0.0, -1.0), PointVortexParams(-1*n_vorts*p.τ, p.a)) + Xs + Xsp +end + + + + + + +""" function point_vortex(pnt::Point3D, cntr::Point3D, p::PointVortexParams) + +Compute the value of a smoothed point vortex at the given point. +""" +function point_vortex(pnt::Point3D, cntr::Point3D, p::PointVortexParams) + gcd = great_circle_dist(pnt,cntr) + p.τ / (cosh(3gcd/p.a)^2) +end + diff --git a/packages/algjulia-service/src/decapodes-service/theory.jl b/packages/algjulia-service/src/decapodes-service/theory.jl new file mode 100644 index 00000000..5e7dbb80 --- /dev/null +++ b/packages/algjulia-service/src/decapodes-service/theory.jl @@ -0,0 +1,129 @@ +""" Helper function to convert CatColab values (Obs) in Decapodes """ +function to_theory(theory::ThDecapode, type::ObType, name::String) + @match lowercase(name) begin + "0-form" => :Form0 + "1-form" => :Form1 + "2-form" => :Form2 + "primal 0-form" => :Form0 + "primal 1-form" => :Form1 + "primal 2-form" => :Form2 + "dual 0-form" => :DualForm0 + "dual 1-form" => :DualForm1 + "dual 2-form" => :DualForm2 + x => throw(ImplError(x)) + end +end + +""" Helper function to convert CatColab values (Homs) in Decapodes """ +function to_theory(theory::ThDecapode, type::HomType, name::String) + @match replace(name," " => "") begin + "∂t" || "∂ₜ" => :∂ₜ + # "∂ₜ" => :∂ₜ + "Δ" => :Δ + "Δ⁻¹" => :Δ⁻¹ + "d*" || "d̃₁" => :dual_d₁ + "⋆" || "⋆₁" || "★₁" || "★1" => :⋆₁ + "⋆⁻¹" || "⋆₀⁻¹" => :⋆₀⁻¹ + "★" || "★⁻¹" => :⋆₁ + "diffusivity" => :diffusivity + "d" || "d₀" || "d01" => :d₀ + "d12" => :d₁ + "⋆2" => :⋆₂ + "♭♯" => :♭♯ + "lamb" => :dpsw # dual-primal self-wedge + "-" => :neg + x => throw(ImplError(x)) + end +end + +# Build the theory + +#= +@active patterns are MLStyle-implementations of F# active patterns that forces us to work in the Maybe/Option pattern. +Practically, yet while a matter of opinion, they make @match statements cleaner; a statement amounts to a helpful pattern +name and the variables we intend to capture. +=# +@active IsObject(x) begin + x[:tag] == "object" ? Some(x) : nothing +end + +@active IsMorphism(x) begin + x[:tag] == "morphism" ? Some(x) : nothing +end + +export IsObject, IsMorphism + +""" Obs, Homs """ +abstract type ElementData end + +""" +Struct capturing the name of the object and its relevant information. +ElementData may be objects or homs, each of which has different data. +""" +struct TheoryElement + name::Union{Symbol, Nothing} + val::Union{ElementData, Nothing} + function TheoryElement(;name::Symbol=nothing,val::Any=nothing) + new(name, val) + end +end +export TheoryElement + +Base.nameof(t::TheoryElement) = t.name + +struct ObData <: ElementData end +# TODO not being used right now but added for completeness. + +struct HomData <: ElementData + dom::Any + cod::Any + function HomData(;dom::Any,cod::Any) + new(dom,cod) + end +end +export HomData +# TODO type dom/cod + +""" Struct wrapping a dictionary """ +struct Theory{T<:ThDecapode} + data::Dict{String, TheoryElement} + function Theory(::T) where T + new{T}(Dict{String, TheoryElement}()) + end +end +export Theory + +Base.values(theory::Theory) = values(theory.data) + +function add_to_theory! end; export add_to_theory! + +function add_to_theory!(theory::Theory{T}, content::AbstractDict, type::ObType) where T + push!(theory.data, + content[:id] => TheoryElement(;name=to_theory(T(), type, content[:name]))) +end + +function add_to_theory!(theory::Theory{T}, content::AbstractDict, type::HomType) where T + push!(theory.data, content[:id] => + TheoryElement(;name=to_theory(T(), type, content[:name]), + val=HomData(dom=content[:dom][:content], + cod=content[:cod][:content]))) +end + +# for each cell, if it is... +# ...an object, we convert its type to a symbol and add it to the theorydict +# ...a morphism, we add it to the theorydict with a field for the ids of its +# domain and codomain to its +function Theory(model::AbstractVector{<:AbstractDict}) + theory = Theory(ThDecapode()) + foreach(model) do cell + @match cell begin + IsObject(content) => add_to_theory!(theory, content, ObType()) + IsMorphism(content) => add_to_theory!(theory, content, HomType()) + x => throw(ImplError(x)) + end + end + return theory +end +export Theory +# TODO parameterize by theory + diff --git a/packages/algjulia-service/src/decapodes.jl b/packages/algjulia-service/src/decapodes.jl deleted file mode 100644 index 3ad51e78..00000000 --- a/packages/algjulia-service/src/decapodes.jl +++ /dev/null @@ -1,331 +0,0 @@ -# - parameterize by the theory. this is currently fixed to decapodes -# - switch to different meshes -# - use enum instead of val - -# NOTES: -# TODO "anonymous objects" are •n - -# algebraicjulia dependencies -using ACSets -using Decapodes -using DiagrammaticEquations -using CombinatorialSpaces - -# dependencies -import JSON3 -using StaticArrays -using MLStyle -using LinearAlgebra -using ComponentArrays -using Distributions # for initial conditions -using GeometryBasics: Point2, Point3 -using OrdinaryDiffEq - -export infer_types!, evalsim, default_dec_generate, default_dec_matrix_generate, DiagonalHodge, ComponentArray - -struct ImplError <: Exception - name::String -end - -Base.showerror(io::IO, e::ImplError) = print(io, "$(e.name) not implemented") - -## THEORY BUILDING - -function to_pode end -export to_pode - -""" Helper function to convert CatColab values (Obs) in Decapodes """ -function to_pode(::Val{:Ob}, name::String) - @match lowercase(name) begin - "0-form" => :Form0 - "1-form" => :Form1 - "2-form" => :Form2 - "dual 0-form" => :DualForm0 - "dual 1-form" => :DualForm1 - "dual 2-form" => :DualForm2 - x => throw(ImplError(x)) - end -end - -""" Helper function to convert CatColab values (Homs) in Decapodes """ -function to_pode(::Val{:Hom}, name::String) - @match name begin - "∂t" => :∂ₜ - "∂ₜ" => :∂ₜ - "Δ" => :Δ - "d" => :d₀ - "d*" => :dual_d₁ - # \star on LHS - "⋆" => :⋆₁ - "⋆⁻¹" => :⋆₀⁻¹ - # \bigstar on LHS - "★" => :⋆₁ - "★⁻¹" => :⋆₀⁻¹ - "diffusivity" => :diffusivity - x => throw(ImplError(x)) - end -end - -# Build the theory - -# @active patterns are MLStyle-implementations of F# active patterns that forces us to work in the Maybe/Option design pattern. They make @match statements cleaner. -@active IsObject(x) begin - x[:tag] == "object" ? Some(x) : nothing -end - -@active IsMorphism(x) begin - x[:tag] == "morphism" ? Some(x) : nothing -end - -export IsObject, IsMorphism - -""" Obs, Homs """ -abstract type ElementData end - -""" Struct capturing the name of the object and its relevant information. ElementData may be objects or homs, each of which has different data. -""" -struct TheoryElement - name::Union{Symbol, Nothing} - val::Union{ElementData, Nothing} - function TheoryElement(;name::Symbol=nothing,val::Any=nothing) - new(name, val) - end -end -export TheoryElement - -Base.nameof(t::TheoryElement) = t.name - -struct HomData <: ElementData - dom::Any - cod::Any - function HomData(;dom::Any,cod::Any) - new(dom,cod) - end -end -export HomData - -struct Theory - data::Dict{String, TheoryElement} - function Theory() - new(Dict{String, TheoryElement}()) - end -end -export Theory - -# TODO engooden -Base.show(io::IO, theory::Theory) = println(io, theory.data) - -Base.values(theory::Theory) = values(theory.data) - -function add_to_theory! end -export add_to_theory! - -function add_to_theory!(theory::Theory, content::Any, type::Val{:Ob}) - push!(theory.data, content[:id] => TheoryElement(;name=to_pode(type, content[:name]))) -end - -function add_to_theory!(theory::Theory, content::Any, type::Val{:Hom}) - push!(theory.data, content[:id] => - TheoryElement(;name=to_pode(type, content[:name]), - val=HomData(dom=content[:dom][:content], - cod=content[:cod][:content]))) -end - -# for each cell, if it is... -# ...an object, we convert its type to a symbol and add it to the theorydict -# ...a morphism, we add it to the theorydict with a field for the ids of its -# domain and codomain to its -function Theory(model::AbstractVector{JSON3.Object}) - theory = Theory(); - foreach(model) do cell - @match cell begin - IsObject(content) => add_to_theory!(theory, content, Val(:Ob)) - IsMorphism(content) => add_to_theory!(theory, content, Val(:Hom)) - x => throw(ImplError(x)) - end - end - return theory -end -export Theory - -## MODEL BUILDING - -function add_to_pode! end -export add_to_pode! - -function add_to_pode!(d::SummationDecapode, - vars::Dict{String, Int}, # mapping between UUID and ACSet ID - theory::Theory, - content::JSON3.Object, - nc::Vector{Int}, - ::Val{:Ob}) - theory_elem = theory.data[content[:over][:content]] # indexes the theory by UUID - name = if isempty(content[:name]) - nc[1] += 1 - Symbol("•$(nc[1])") - else - Symbol(content[:name]) - end - id = add_part!(d, :Var, name=name, type=nameof(theory_elem)) - push!(vars, content[:id] => id) - return d -end - -# TODO we are restricted to Op1 -function add_to_pode!(d::SummationDecapode, - vars::Dict{String, Int}, # mapping between UUID and ACSet ID - theory::Theory, - content::JSON3.Object, - scalars::Any, - anons::Dict{Symbol, Any}, - ::Val{:Hom}) - dom = content[:dom][:content] - cod = content[:cod][:content] - if haskey(vars, dom) && haskey(vars, cod) - op1 = Symbol(theory.data[content[:over][:content]].name) - add_part!(d, :Op1, src=vars[dom], tgt=vars[cod], op1=op1) - # we need to add an inclusion to the TVar table - if op1 == :∂ₜ - add_part!(d, :TVar, incl=vars[cod]) - end - # if the dom is anonymous, we treat it as a something which will receive x -> k * x. - # we store its value in another array - if !isempty(scalars) && haskey(scalars, Symbol(content[:over][:content])) - scalar = scalars[Symbol(content[:over][:content])] - push!(anons, op1 => x -> scalar * x) - end - # TODO if scalars were typed correctly, we could probably do away with the !isempty check - end - d -end - -""" Decapode(jsondiagram::JSON3.Object, theory::Theory) => SummationDecapode - -This returns a Decapode given a jsondiagram and a theory. -""" -function Decapode(diagram::AbstractVector{JSON3.Object}, theory::Theory; scalars=[]) - # initiatize decapode and its mapping between UUIDs and ACSet IDs - pode = SummationDecapode(parse_decapode(quote end)); - vars = Dict{String, Int}(); - nc = [0]; # array is a mutable container - anons = Dict{Symbol, Any}(); - - # for each cell in the notebook, add it to the diagram - foreach(diagram) do cell - @match cell begin - # TODO merge nameless_count into vars - IsObject(content) => add_to_pode!(pode, vars, theory, content, nc, Val(:Ob)) - IsMorphism(content) => add_to_pode!(pode, vars, theory, content, scalars, anons, Val(:Hom)) - _ => throw(ImplError(cell[:content][:tag])) - end - end - pode, anons -end -export Decapode -# the proper name for this constructor should be "SummationDecapode" - -# TODO we need to make this dynamic -function create_mesh(statevar::Symbol) - s = triangulated_grid(100,100,2,2,Point2{Float64}) - sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point2{Float64}}(s) - subdivide_duals!(sd, Circumcenter()) - - c_dist = MvNormal([50, 50], LinearAlgebra.Diagonal(map(abs2, [7.5, 7.5]))) - c = [pdf(c_dist, [p[1], p[2]]) for p in sd[:point]] - u0 = ComponentArray(; Dict(statevar=>c)...) - - return (s, sd, u0, ()) -end -export create_mesh - - -struct System - pode::SummationDecapode - statevar::Symbol - scalars::Dict{Symbol, Any} # closures # TODO rename scalars => anons - mesh::HasDeltaSet - dualmesh::HasDeltaSet - init::Any # TODO specify. Is it always ComponentVector? - generate::Any -end -export System - -function System(json_string::String) - json_object = JSON3.read(json_string); - - # converts the JSON of (the fragment of) the theory - # into theory of the DEC, valued in Julia - theory = Theory(json_object[:model]); - - # this is a diagram in the model of the DEC. it wants to be a decapode! - diagram = json_object[:diagram]; - - # any scalars? - scalars = haskey(json_object, :scalars) ? json_object[:scalars] : []; - - # pode, anons, and statevar - decapode, anons = Decapode(diagram, theory; scalars=scalars); - statevars = infer_state_names(decapode); - statevar = length(statevars) == 1 ? first(statevars) : error("$statevars must be length one") - - # mesh and initial conditions - s, sd, u0, _ = create_mesh(statevar); - - # generate - function sys_generate(s, my_symbol; hodge=GeometricHodge()) - op = @match my_symbol begin - sym && if sym ∈ keys(anons) end => anons[sym] - _ => default_dec_matrix_generate(s, my_symbol, hodge) - end - return (args...) -> op(args...) - end - - return System(decapode, statevar, anons, s, sd, u0, sys_generate) -end - - -function run_sim(fm, u0, t0, constparam) - prob = ODEProblem(fm, u0, (0, t0), constparam) - soln = solve(prob, Tsit5(), saveat=0.1) -end -export run_sim - -struct SimResult - time::Vector{Float64} - state::Vector{Matrix{SVector{3, Float64}}} - x::Vector{Float64} # axis - y::Vector{Float64} -end -export SimResult - -function SimResult(sol::ODESolution, system::System) - - points = collect(values(system.mesh.subparts.point.m)); - - xlen = 51; ylen = 51; - - function at_time(sol::ODESolution, timeidx::Int) - [SVector(i, j, getproperty(sol.u[timeidx], system.statevar)[xlen*(i-1) + j]) for i in 1:xlen, j in 1:ylen] - end - - state_vals = map(1:length(sol.t)) do i - at_time(sol, i) - end - - # TODO engooden - SimResult(sol.t, state_vals, 0:xlen-1, 0:ylen-1) -end -# TODO generalize to HasDeltaSet - -## PLOTTING CODE - -abstract type AbstractPlotType end - -struct Heatmap <: AbstractPlotType end - -# TODO make length a conditional value so we can pass it in if we want -function Base.reshape(::Heatmap, data) - l = floor(Int64, sqrt(length(data))) - reshape(data, l, l) -end - diff --git a/packages/algjulia-service/test/Project.toml b/packages/algjulia-service/test/Project.toml index c12e2820..a2ed5fe7 100644 --- a/packages/algjulia-service/test/Project.toml +++ b/packages/algjulia-service/test/Project.toml @@ -9,6 +9,5 @@ JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/packages/algjulia-service/test/runtests.jl b/packages/algjulia-service/test/runtests.jl index 86f4cb89..76f00592 100644 --- a/packages/algjulia-service/test/runtests.jl +++ b/packages/algjulia-service/test/runtests.jl @@ -1,11 +1,11 @@ using Test - +# using AlgebraicJuliaService using ACSets using CombinatorialSpaces using Decapodes using DiagrammaticEquations - +# using MLStyle using JSON3 using ComponentArrays @@ -13,34 +13,34 @@ using StaticArrays using LinearAlgebra import OrdinaryDiffEq: ReturnCode -# visualization -using Plots +const KEYS = Set([:mesh, :plotVariables, :initialConditions, :domain, :diagram, :model, :scalars]) # load data -data = open(JSON3.read, joinpath(@__DIR__, "diffusion_data.json"), "r") -diagram = data[:diagram]; -model = data[:model]; +data = open(JSON3.read, joinpath(@__DIR__, "test_jsons", "diffusion_data.json"), "r") +diagram = data[:diagram] +model = data[:model] @testset "Text-to-Pode" begin - @test to_pode(Val(:Ob), "0-form") == :Form0 - @test to_pode(Val(:Ob), "1-form") == :Form1 - @test to_pode(Val(:Ob), "2-form") == :Form2 - @test to_pode(Val(:Ob), "dual 0-form") == :DualForm0 - @test to_pode(Val(:Ob), "dual 1-form") == :DualForm1 - @test to_pode(Val(:Ob), "dual 2-form") == :DualForm2 + @test to_theory(ThDecapode(), ObType(), "0-form") == :Form0 + @test to_theory(ThDecapode(), ObType(), "1-form") == :Form1 + @test to_theory(ThDecapode(), ObType(), "2-form") == :Form2 + @test to_theory(ThDecapode(), ObType(), "dual 0-form") == :DualForm0 + @test to_theory(ThDecapode(), ObType(), "dual 1-form") == :DualForm1 + @test to_theory(ThDecapode(), ObType(), "dual 2-form") == :DualForm2 - @test_throws AlgebraicJuliaService.ImplError to_pode(Val(:Ob), "Form3") + @test_throws AlgebraicJuliaService.ImplError to_theory(ThDecapode(), ObType(), "Form3") - @test to_pode(Val(:Hom), "∂t") == :∂ₜ - @test to_pode(Val(:Hom), "Δ") == :Δ - @test_throws AlgebraicJuliaService.ImplError to_pode(Val(:Hom), "∧") + @test to_theory(ThDecapode(), HomType(), "∂t") == :∂ₜ + @test to_theory(ThDecapode(), HomType(), "Δ") == :Δ + @test_throws AlgebraicJuliaService.ImplError to_theory(ThDecapode(), HomType(), "∧") end + @testset "Parsing the Theory JSON Object" begin - @test Set(keys(data)) == Set([:diagram, :model]) + @test Set(keys(data)) == KEYS @test @match model[1] begin IsObject(_) => true @@ -52,70 +52,89 @@ end _ => false end - theory = Theory(); + theory = Theory(ThDecapode()) @match model[1] begin - IsObject(content) => add_to_theory!(theory, content, Val(:Ob)) + IsObject(content) => add_to_theory!(theory, content, ObType()) _ => nothing end - @test theory.data["019323fa-49cb-7373-8c5d-c395bae4006d"] == TheoryElement(;name=:Form0, val=nothing) + + _id = "019323fa-49cb-7373-8c5d-c395bae4006d" + @test theory.data[_id] == TheoryElement(;name=:Form0, val=nothing) end @testset "Making the Decapode" begin - theory = Theory(model); + theory = Theory(model) @test Set(nameof.(values(theory))) == Set([:Form0, :Form1, :Form2, :Δ, :∂ₜ]) - handcrafted_pode = SummationDecapode(parse_decapode(quote end)); - add_part!(handcrafted_pode, :Var, name=:C, type=:Form0); - add_part!(handcrafted_pode, :Var, name=Symbol("dC/dt"), type=:Form0); - add_part!(handcrafted_pode, :TVar, incl=2); - add_part!(handcrafted_pode, :Op1, src=1, tgt=2, op1=:∂ₜ); - add_part!(handcrafted_pode, :Op1, src=1, tgt=2, op1=:Δ); + handcrafted_pode = SummationDecapode(parse_decapode(quote end)) + add_part!(handcrafted_pode, :Var, name=:C, type=:Form0) + add_part!(handcrafted_pode, :Var, name=Symbol("dC/dt"), type=:Form0) + add_part!(handcrafted_pode, :TVar, incl=2) + add_part!(handcrafted_pode, :Op1, src=1, tgt=2, op1=:∂ₜ) + add_part!(handcrafted_pode, :Op1, src=1, tgt=2, op1=:Δ) # no scalars in second position - decapode, _ = Decapode(diagram, theory); + decapode, _, _ = Decapode(diagram, theory) @test decapode == handcrafted_pode end -@testset "Simulation" begin +@testset "Simulation - Diffusion" begin - json_string = read(joinpath(@__DIR__, "diffusion_data.json"), String); - system = System(json_string); + json_string = read(joinpath(@__DIR__, "test_jsons", "diffusion_data.json"), String) + @test Set(keys(JSON3.read(json_string))) == KEYS + system = PodeSystem(json_string) + simulator = evalsim(system.pode) - f = simulator(system.dualmesh, default_dec_generate, DiagonalHodge()); + f = simulator(system.geometry.dualmesh, default_dec_generate, DiagonalHodge()) - # time soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)); - # returns ::ODESolution - # - retcode - # - interpolation - # - t - # - u::Vector{ComponentVector} @test soln.retcode == ReturnCode.Success - result = SimResult(soln, system); + result = SimResult(soln, system) - @test typeof(result.state) == Vector{Matrix{SVector{3, Float64}}} + @test typeof(result.state) == Dict{String, Vector{AbstractArray{SVector{3, Float64}}}} - jv = JsonValue(result); + jv = JsonValue(result) end +#= XXX ERRORING +we are trying to index `soln.u[t]` by `Ċ `, but only `C` is present. I think this is because we generating initial conditions based on what was specified to the `initial_conditions` parameter, whereas in the past we were using `infer_state_names` to obtain that. +=# +# @testset "Simulation - Diffusion with Two Variables" begin -##### +# json_string = read(joinpath(@__DIR__, "test_jsons", "diffusion_data_twovars.json"), String); +# @test Set(keys(JSON3.read(json_string))) == KEYS -data = open(JSON3.read, joinpath(@__DIR__, "diffusion_long_trip.json"), "r") -diagram = data[:diagram]; -model = data[:model]; +# system = PodeSystem(json_string); -@testset "Parsing the Theory JSON Object" begin +# simulator = evalsim(system.pode) +# f = simulator(system.geometry.dualmesh, system.generate, DiagonalHodge()) + +# soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)); + +# @test soln.retcode == ReturnCode.Success + +# result = SimResult(soln, system); + +# @test typeof(result.state) == Dict{String, Vector{AbstractArray{SVector{3, Float64}}}} + +# jvs = JsonValue(result); - @test Set(keys(data)) == Set([:diagram, :model, :scalars]) +# end + +@testset "Parsing the Theory JSON Object - Diffusion Long Trip" begin + + data = open(JSON3.read, joinpath(@__DIR__, "test_jsons", "diffusion_long_trip.json"), "r") + diagram = data[:diagram] + model = data[:model] + @test Set(keys(data)) == KEYS @test @match model[1] begin IsObject(_) => true @@ -127,9 +146,9 @@ model = data[:model]; _ => false end - theory = Theory(); + theory = Theory(ThDecapode()) @match model[1] begin - IsObject(content) => add_to_theory!(theory, content, Val(:Ob)) + IsObject(content) => add_to_theory!(theory, content, ObType()) _ => nothing end @@ -137,20 +156,19 @@ model = data[:model]; end -@testset "Simulation ..." begin +# # GOOD +@testset "Simulation - Diffusion Long Trip" begin - json_string = read(joinpath(@__DIR__, "diffusion_long_trip.json"), String); - system = System(json_string); + json_string = read(joinpath(@__DIR__, "test_jsons", "diffusion_long_trip.json"), String) + @test Set(keys(JSON3.read(json_string))) == KEYS + + system = PodeSystem(json_string) simulator = evalsim(system.pode) - # open("test_sim.jl", "w") do f - # write(f, string(gensim(system.pode))) - # end - # simulator = include("test_sim.jl") - f = simulator(system.dualmesh, default_dec_matrix_generate, DiagonalHodge()); + f = simulator(system.geometry.dualmesh, default_dec_matrix_generate, DiagonalHodge()) # time - soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)); + soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)) # returns ::ODESolution # - retcode # - interpolation @@ -159,24 +177,23 @@ end @test soln.retcode == ReturnCode.Success - result = SimResult(soln, system); + result = SimResult(soln, system) - @test typeof(result.state) == Vector{Matrix{SVector{3, Float64}}} + @test typeof(result.state) == Dict{String, Vector{AbstractArray{SVector{3, Float64}}}} - jv = JsonValue(result); + jv = JsonValue(result) end -### - -data = open(JSON3.read, joinpath(@__DIR__, "diffusivity_constant.json"), "r") -diagram = data[:diagram]; -model = data[:model]; -scalars = data[:scalars]; - -@testset "Parsing the Theory JSON Object" begin +# # GOOD +@testset "Parsing the Theory JSON Object - Diffusivity Constant" begin + + data = open(JSON3.read, joinpath(@__DIR__, "test_jsons", "diffusivity_constant.json"), "r") + diagram = data[:diagram] + model = data[:model] + scalars = data[:scalars] - @test Set(keys(data)) == Set([:diagram, :model, :scalars]) + @test Set(keys(data)) == KEYS @test @match model[1] begin IsObject(_) => true @@ -188,9 +205,9 @@ scalars = data[:scalars]; _ => false end - theory = Theory(); + theory = Theory(ThDecapode()) @match model[1] begin - IsObject(content) => add_to_theory!(theory, content, Val(:Ob)) + IsObject(content) => add_to_theory!(theory, content, ObType()) _ => nothing end @@ -198,51 +215,46 @@ scalars = data[:scalars]; end -@testset "Simulation ..." begin +@testset "Simulation - Diffusivity Constant" begin - json_string = read(joinpath(@__DIR__, "diffusivity_constant.json"), String); - system = System(json_string); + json_string = read(joinpath(@__DIR__, "test_jsons", "diffusivity_constant.json"), String) + system = PodeSystem(json_string) simulator = evalsim(system.pode) - # open("test_sim.jl", "w") do f - # write(f, string(gensim(system.pode))) - # end - # simulator = include("test_sim.jl"); - - f = simulator(system.dualmesh, system.generate, DiagonalHodge()); + + f = simulator(system.geometry.dualmesh, system.generate, DiagonalHodge()) - soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)); - # returns ::ODESolution - # - retcode - # - interpolation - # - t - # - u::Vector{ComponentVector} + soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)) @test soln.retcode == ReturnCode.Success - result = SimResult(soln, system); + result = SimResult(soln, system) - @test typeof(result.state) == Vector{Matrix{SVector{3, Float64}}} + @test typeof(result.state) == Dict{String, Vector{AbstractArray{SVector{3, Float64}}}} - jv = JsonValue(result); + jv = JsonValue(result) end -# PLOTTING UTILITIES +@testset "Simulation - Navier-Stokes Vorticity" begin -# TODO size fixed -# function at_time(sr::SimResult, t::Int) -# [sr.state[t][i,j][3] for i in 1:51 for j in 1:51] -# end + json_string = read(joinpath(@__DIR__, "test_jsons", "ns_vorticity.json"), String) + @test Set(keys(JSON3.read(json_string))) == KEYS -# function show_heatmap(sr::SimResult, t::Int) -# data = at_time(result, t) -# ℓ = floor(Int64, sqrt(length(data))); -# reshaped = reshape(data, ℓ, ℓ) -# Plots.heatmap(1:51, 1:51, reshaped, clims=(minimum(data), maximum(data)); palette=:redsblues) -# end + system = PodeSystem(json_string) + + simulator = evalsim(system.pode) + + f = simulator(system.geometry.dualmesh, system.generate, DiagonalHodge()) + + soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)) -# @gif for t ∈ 1:length(result.time) -# show_heatmap(result, t) -# end every 5 + @test soln.retcode == ReturnCode.Success + + result = SimResult(soln, system); + @test typeof(result.state) == Dict{String, Vector{AbstractArray{SVector{3, Float64}}}} + + jv = JsonValue(result) + +end diff --git a/packages/algjulia-service/test/diffusion_data.json b/packages/algjulia-service/test/test_jsons/diffusion_data.json similarity index 93% rename from packages/algjulia-service/test/diffusion_data.json rename to packages/algjulia-service/test/test_jsons/diffusion_data.json index 452b2435..36ab0a1a 100644 --- a/packages/algjulia-service/test/diffusion_data.json +++ b/packages/algjulia-service/test/test_jsons/diffusion_data.json @@ -1,4 +1,11 @@ -{ +{ + "domain": "Plane", + "mesh": "Rectangle", + "initialConditions": { + "01932402-bcf5-7432-8d14-dbae9eabf907": "Gaussian" + }, + "scalars": {}, + "plotVariables": ["01932402-bcf5-7432-8d14-dbae9eabf907"], "diagram": [ { "id": "01932402-bcf5-7432-8d14-dbae9eabf907", diff --git a/packages/algjulia-service/test/test_jsons/diffusion_data_twovars.json b/packages/algjulia-service/test/test_jsons/diffusion_data_twovars.json new file mode 100644 index 00000000..56dacb6f --- /dev/null +++ b/packages/algjulia-service/test/test_jsons/diffusion_data_twovars.json @@ -0,0 +1,154 @@ +{ + "domain": "Plane", + "mesh": "Rectangle", + "initialConditions": { + "01932402-bcf5-7432-8d14-dbae9eabf907": "Gaussian" + }, + "scalars": {}, + "plotVariables": ["01932402-bcf5-7432-8d14-dbae9eabf907","01932403-5b6c-7231-90d7-d7cece275eb2"], + "diagram": [ + { + "id": "01932402-bcf5-7432-8d14-dbae9eabf907", + "name": "C", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "over": { + "content": "019323fa-49cb-7373-8c5d-c395bae4006d", + "tag": "Basic" + }, + "tag": "object" + }, + { + "id": "01932403-5b6c-7231-90d7-d7cece275eb2", + "name": "dC/dt", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "over": { + "content": "019323fa-49cb-7373-8c5d-c395bae4006d", + "tag": "Basic" + }, + "tag": "object" + }, + { + "cod": { + "content": "01932403-5b6c-7231-90d7-d7cece275eb2", + "tag": "Basic" + }, + "dom": { + "content": "01932402-bcf5-7432-8d14-dbae9eabf907", + "tag": "Basic" + }, + "id": "01932403-c4cd-7563-8ebd-080dd37a9c7e", + "morType": { + "content": { + "content": "Object", + "tag": "Basic" + }, + "tag": "Hom" + }, + "name": "", + "over": { + "content": "019323fb-3652-7e91-aee9-06187a954fc6", + "tag": "Basic" + }, + "tag": "morphism" + }, + { + "cod": { + "content": "01932403-5b6c-7231-90d7-d7cece275eb2", + "tag": "Basic" + }, + "dom": { + "content": "01932402-bcf5-7432-8d14-dbae9eabf907", + "tag": "Basic" + }, + "id": "01932404-10e5-7128-bb94-835e5d8d643f", + "morType": { + "content": { + "content": "Object", + "tag": "Basic" + }, + "tag": "Hom" + }, + "name": "", + "over": { + "content": "019323ff-1af6-79da-b776-8ee11c88a8a0", + "tag": "Basic" + }, + "tag": "morphism" + } + ], + "model": [ + { + "id": "019323fa-49cb-7373-8c5d-c395bae4006d", + "name": "0-form", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "tag": "object" + }, + { + "id": "019323fa-783b-72c8-af20-c0718fde3ac8", + "name": "1-form", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "tag": "object" + }, + { + "id": "019323fb-175b-784e-aab8-7b78fa576571", + "name": "2-form", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "tag": "object" + }, + { + "cod": { + "content": "019323fa-49cb-7373-8c5d-c395bae4006d", + "tag": "Basic" + }, + "dom": { + "content": "019323fa-49cb-7373-8c5d-c395bae4006d", + "tag": "Basic" + }, + "id": "019323fb-3652-7e91-aee9-06187a954fc6", + "morType": { + "content": { + "content": "Object", + "tag": "Basic" + }, + "tag": "Hom" + }, + "name": "∂t", + "tag": "morphism" + }, + { + "cod": { + "content": "019323fa-49cb-7373-8c5d-c395bae4006d", + "tag": "Basic" + }, + "dom": { + "content": "019323fa-49cb-7373-8c5d-c395bae4006d", + "tag": "Basic" + }, + "id": "019323ff-1af6-79da-b776-8ee11c88a8a0", + "morType": { + "content": { + "content": "Object", + "tag": "Basic" + }, + "tag": "Hom" + }, + "name": "Δ", + "tag": "morphism" + } + ] +} diff --git a/packages/algjulia-service/test/diffusion_long_trip.json b/packages/algjulia-service/test/test_jsons/diffusion_long_trip.json similarity index 97% rename from packages/algjulia-service/test/diffusion_long_trip.json rename to packages/algjulia-service/test/test_jsons/diffusion_long_trip.json index 08050f8c..6e8eaec4 100644 --- a/packages/algjulia-service/test/diffusion_long_trip.json +++ b/packages/algjulia-service/test/test_jsons/diffusion_long_trip.json @@ -1,4 +1,10 @@ -{ +{ + "domain": "Plane", + "mesh": "Rectangle", + "initialConditions": { + "01936ee4-8a1d-7bd7-a2c4-0f5291893346": "Gaussian" + }, + "plotVariables": ["01936ee4-8a1d-7bd7-a2c4-0f5291893346"], "diagram": [ { "tag": "object", diff --git a/packages/algjulia-service/test/diffusivity_constant.json b/packages/algjulia-service/test/test_jsons/diffusivity_constant.json similarity index 95% rename from packages/algjulia-service/test/diffusivity_constant.json rename to packages/algjulia-service/test/test_jsons/diffusivity_constant.json index fe132633..4fbdb543 100644 --- a/packages/algjulia-service/test/diffusivity_constant.json +++ b/packages/algjulia-service/test/test_jsons/diffusivity_constant.json @@ -1,4 +1,10 @@ -{ +{ + "domain": "Plane", + "mesh": "Rectangle", + "initialConditions": { + "01936f2e-2f39-738d-8af2-5caf75b29f3a": "Gaussian" + }, + "plotVariables": ["01936f2e-2f39-738d-8af2-5caf75b29f3a"], "diagram": [ { "tag": "object", diff --git a/packages/algjulia-service/test/example.json b/packages/algjulia-service/test/test_jsons/example.json similarity index 99% rename from packages/algjulia-service/test/example.json rename to packages/algjulia-service/test/test_jsons/example.json index 3f30feca..d74244d4 100644 --- a/packages/algjulia-service/test/example.json +++ b/packages/algjulia-service/test/test_jsons/example.json @@ -1,4 +1,4 @@ -{ +{ "plotVariables": ["C"], "diagram": [ { "tag": "object", diff --git a/packages/algjulia-service/test/test_jsons/ns_vorticity.json b/packages/algjulia-service/test/test_jsons/ns_vorticity.json new file mode 100644 index 00000000..55a46072 --- /dev/null +++ b/packages/algjulia-service/test/test_jsons/ns_vorticity.json @@ -0,0 +1,544 @@ +{ + "diagram": [ + { + "tag": "object", + "name": "", + "id": "01938869-c0ae-7b1f-9179-35d64102bc6e", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938829-293b-726a-8e86-f7462a8f6f07" + } + }, + { + "tag": "object", + "name": "uu", + "id": "01938867-e117-7be0-b32a-1e4cd0af5cd9", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-ac3d-711a-a508-35087405d2eb" + } + }, + { + "tag": "object", + "name": "", + "id": "0193886f-3fe1-740e-9a0e-2f3d8add2fa1", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-b4c3-7996-850f-2fecb9dc5343" + } + }, + { + "tag": "object", + "name": "", + "id": "0193886b-bba3-7b84-9311-4412f4c0718f", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-ac3d-711a-a508-35087405d2eb" + } + }, + { + "tag": "object", + "name": "", + "id": "01938868-dfd1-74b2-8b89-32d599ef94e6", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-972c-7fd6-a4f7-bbf1e0ee52fc" + } + }, + { + "tag": "object", + "name": "", + "id": "01938dad-a5b0-7ce3-a0bf-615e4111ce02", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-b4c3-7996-850f-2fecb9dc5343" + } + }, + { + "tag": "object", + "name": "duu", + "id": "01938fb0-449c-7614-8086-2384cf6bc784", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-b4c3-7996-850f-2fecb9dc5343" + } + }, + { + "tag": "object", + "name": "", + "id": "0193886c-f9bd-7872-8d57-1e3c88e7c31b", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-ac3d-711a-a508-35087405d2eb" + } + }, + { + "tag": "object", + "name": "", + "id": "0193886c-55b4-71a9-816c-3c30896be5d6", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938829-293b-726a-8e86-f7462a8f6f07" + } + }, + { + "tag": "object", + "name": "ψ", + "id": "01938867-f97d-7189-97e5-e5cd19500108", + "obType": { + "tag": "Basic", + "content": "Object" + }, + "over": { + "tag": "Basic", + "content": "01938827-972c-7fd6-a4f7-bbf1e0ee52fc" + } + }, + { + "tag": "morphism", + "name": "", + "id": "01938868-f03f-7aed-bd6a-e78476e29a99", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "01938828-bdf1-7d20-9dd7-7ad4c8dd78bf" + }, + "dom": { + "tag": "Basic", + "content": "01938868-dfd1-74b2-8b89-32d599ef94e6" + }, + "cod": { + "tag": "Basic", + "content": "01938867-f97d-7189-97e5-e5cd19500108" + } + }, + { + "tag": "morphism", + "name": "", + "id": "01938869-a775-7129-9daa-40b8709365b9", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "01938829-0414-7b82-a0f9-f714d17e6ff4" + }, + "dom": { + "tag": "Basic", + "content": "01938867-f97d-7189-97e5-e5cd19500108" + }, + "cod": { + "tag": "Basic", + "content": "01938869-c0ae-7b1f-9179-35d64102bc6e" + } + }, + { + "tag": "morphism", + "name": "", + "id": "0193886c-3f95-7922-919b-7cd9a582d312", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "0193885a-f4e6-7dc9-9ca7-9b318d93b2c0" + }, + "dom": { + "tag": "Basic", + "content": "0193886b-bba3-7b84-9311-4412f4c0718f" + }, + "cod": { + "tag": "Basic", + "content": "0193886c-55b4-71a9-816c-3c30896be5d6" + } + }, + { + "tag": "morphism", + "name": "", + "id": "0193886a-1fbb-7049-ad9e-5d78097f2045", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "01938829-62c4-7897-b6d6-0a3563ee643d" + }, + "dom": { + "tag": "Basic", + "content": "01938869-c0ae-7b1f-9179-35d64102bc6e" + }, + "cod": { + "tag": "Basic", + "content": "01938867-e117-7be0-b32a-1e4cd0af5cd9" + } + }, + { + "tag": "morphism", + "name": "", + "id": "0193886e-eb24-767e-ab8e-1b346b0c1797", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "0193886b-48c1-788c-9ec4-767eaa555fb6" + }, + "dom": { + "tag": "Basic", + "content": "01938fb0-449c-7614-8086-2384cf6bc784" + }, + "cod": { + "tag": "Basic", + "content": "01938dad-a5b0-7ce3-a0bf-615e4111ce02" + } + }, + { + "tag": "morphism", + "name": "", + "id": "01938868-8384-7b3c-a56b-d5e4e91d3ee6", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "01938828-8662-7ac6-8e68-58c0cd035384" + }, + "dom": { + "tag": "Basic", + "content": "01938fb0-449c-7614-8086-2384cf6bc784" + }, + "cod": { + "tag": "Basic", + "content": "01938868-dfd1-74b2-8b89-32d599ef94e6" + } + }, + { + "tag": "morphism", + "name": "", + "id": "0193886b-a079-7129-bf3e-6af6be5907b4", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "01938829-d95c-74be-ba4d-17a19a569a64" + }, + "dom": { + "tag": "Basic", + "content": "01938867-e117-7be0-b32a-1e4cd0af5cd9" + }, + "cod": { + "tag": "Basic", + "content": "0193886b-bba3-7b84-9311-4412f4c0718f" + } + }, + { + "tag": "morphism", + "name": "", + "id": "01938dad-834e-708a-9439-d76797e399ce", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "01938dac-e35f-77ed-9081-6d926d3e366f" + }, + "dom": { + "tag": "Basic", + "content": "0193886f-3fe1-740e-9a0e-2f3d8add2fa1" + }, + "cod": { + "tag": "Basic", + "content": "01938dad-a5b0-7ce3-a0bf-615e4111ce02" + } + }, + { + "tag": "morphism", + "name": "", + "id": "0193886c-9d8a-7c1a-b32a-a7f03283df8d", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "01938829-62c4-7897-b6d6-0a3563ee643d" + }, + "dom": { + "tag": "Basic", + "content": "0193886c-55b4-71a9-816c-3c30896be5d6" + }, + "cod": { + "tag": "Basic", + "content": "0193886c-f9bd-7872-8d57-1e3c88e7c31b" + } + }, + { + "tag": "morphism", + "name": "", + "id": "0193886f-0477-7b28-9f71-ed4ea8014758", + "morType": { + "tag": "Basic", + "content": "Nonscalar" + }, + "over": { + "tag": "Basic", + "content": "0193882a-9357-7283-97e7-ad1171ae957e" + }, + "dom": { + "tag": "Basic", + "content": "0193886c-f9bd-7872-8d57-1e3c88e7c31b" + }, + "cod": { + "tag": "Basic", + "content": "0193886f-3fe1-740e-9a0e-2f3d8add2fa1" + } + } + ], + "model": [ + { + "id": "01938827-972c-7fd6-a4f7-bbf1e0ee52fc", + "name": "0-Form", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "tag": "object" + }, + { + "id": "01938827-ac3d-711a-a508-35087405d2eb", + "name": "Dual 1-Form", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "tag": "object" + }, + { + "id": "01938827-b4c3-7996-850f-2fecb9dc5343", + "name": "Dual 2-form", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "tag": "object" + }, + { + "cod": { + "content": "01938827-972c-7fd6-a4f7-bbf1e0ee52fc", + "tag": "Basic" + }, + "dom": { + "content": "01938827-b4c3-7996-850f-2fecb9dc5343", + "tag": "Basic" + }, + "id": "01938828-8662-7ac6-8e68-58c0cd035384", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "⋆₀⁻¹", + "tag": "morphism" + }, + { + "cod": { + "content": "01938827-972c-7fd6-a4f7-bbf1e0ee52fc", + "tag": "Basic" + }, + "dom": { + "content": "01938827-972c-7fd6-a4f7-bbf1e0ee52fc", + "tag": "Basic" + }, + "id": "01938828-bdf1-7d20-9dd7-7ad4c8dd78bf", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "Δ⁻¹", + "tag": "morphism" + }, + { + "cod": { + "content": "01938829-293b-726a-8e86-f7462a8f6f07", + "tag": "Basic" + }, + "dom": { + "content": "01938827-972c-7fd6-a4f7-bbf1e0ee52fc", + "tag": "Basic" + }, + "id": "01938829-0414-7b82-a0f9-f714d17e6ff4", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "d01", + "tag": "morphism" + }, + { + "id": "01938829-293b-726a-8e86-f7462a8f6f07", + "name": "Primal 1-Form", + "obType": { + "content": "Object", + "tag": "Basic" + }, + "tag": "object" + }, + { + "cod": { + "content": "01938827-ac3d-711a-a508-35087405d2eb", + "tag": "Basic" + }, + "dom": { + "content": "01938829-293b-726a-8e86-f7462a8f6f07", + "tag": "Basic" + }, + "id": "01938829-62c4-7897-b6d6-0a3563ee643d", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "⋆₁", + "tag": "morphism" + }, + { + "cod": { + "content": "01938827-ac3d-711a-a508-35087405d2eb", + "tag": "Basic" + }, + "dom": { + "content": "01938827-ac3d-711a-a508-35087405d2eb", + "tag": "Basic" + }, + "id": "01938829-d95c-74be-ba4d-17a19a569a64", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "lamb", + "tag": "morphism" + }, + { + "cod": { + "content": "01938827-b4c3-7996-850f-2fecb9dc5343", + "tag": "Basic" + }, + "dom": { + "content": "01938827-ac3d-711a-a508-35087405d2eb", + "tag": "Basic" + }, + "id": "0193882a-9357-7283-97e7-ad1171ae957e", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "d̃₁", + "tag": "morphism" + }, + { + "cod": { + "content": "01938829-293b-726a-8e86-f7462a8f6f07", + "tag": "Basic" + }, + "dom": { + "content": "01938827-ac3d-711a-a508-35087405d2eb", + "tag": "Basic" + }, + "id": "0193885a-f4e6-7dc9-9ca7-9b318d93b2c0", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "♭♯", + "tag": "morphism" + }, + { + "cod": { + "content": "01938827-b4c3-7996-850f-2fecb9dc5343", + "tag": "Basic" + }, + "dom": { + "content": "01938827-b4c3-7996-850f-2fecb9dc5343", + "tag": "Basic" + }, + "id": "0193886b-48c1-788c-9ec4-767eaa555fb6", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": " ∂ₜ", + "tag": "morphism" + }, + { + "cod": { + "content": "01938827-b4c3-7996-850f-2fecb9dc5343", + "tag": "Basic" + }, + "dom": { + "content": "01938827-b4c3-7996-850f-2fecb9dc5343", + "tag": "Basic" + }, + "id": "01938dac-e35f-77ed-9081-6d926d3e366f", + "morType": { + "content": "Nonscalar", + "tag": "Basic" + }, + "name": "-", + "tag": "morphism" + } + ], + "domain": "Sphere", + "mesh": "Icosphere6", + "initialConditions": { + "01938fb0-449c-7614-8086-2384cf6bc784": "TaylorVortex" + }, + "plotVariables": [ + "01938fb0-449c-7614-8086-2384cf6bc784" + ], + "scalars": {} +} \ No newline at end of file diff --git a/packages/catlog/src/simulate/ode/dynamic.rs b/packages/catlog/src/simulate/ode/dynamic.rs index 38f5d2a8..e71d112c 100644 --- a/packages/catlog/src/simulate/ode/dynamic.rs +++ b/packages/catlog/src/simulate/ode/dynamic.rs @@ -30,7 +30,7 @@ impl<'a, 'b> VectorFieldEnv<'a, 'b> { } } -impl<'a, 'b> Env for VectorFieldEnv<'a, 'b> { +impl Env for VectorFieldEnv<'_, '_> { type Var = Quantity; fn lookup(&self, t: &Self::Var) -> f32 { diff --git a/packages/catlog/src/stdlib/analyses/ode/stock_flow.rs b/packages/catlog/src/stdlib/analyses/ode/stock_flow.rs index 88ca4d9e..f26b4911 100644 --- a/packages/catlog/src/stdlib/analyses/ode/stock_flow.rs +++ b/packages/catlog/src/stdlib/analyses/ode/stock_flow.rs @@ -170,7 +170,7 @@ pub struct StockFlowSystem { struct DVectorEnv<'a>(&'a DVector); -impl<'a> Env for DVectorEnv<'a> { +impl Env for DVectorEnv<'_> { type Var = usize; fn lookup(&self, t: &Self::Var) -> f32 { diff --git a/packages/frontend/src/components/fixed_table_editor.tsx b/packages/frontend/src/components/fixed_table_editor.tsx index 54eff3a6..2545d186 100644 --- a/packages/frontend/src/components/fixed_table_editor.tsx +++ b/packages/frontend/src/components/fixed_table_editor.tsx @@ -3,15 +3,27 @@ import { For, Match, Show, Switch, createEffect, createSignal } from "solid-js"; import "./fixed_table_editor.css"; -/** Schema for a column in a table editor. */ -export type ColumnSchema = { +type ContentType = "string" | "boolean" | "enum"; + +type BaseColumnSchema = { + /** Type of content displayed in the column. */ + contentType: ContentType; + /** Name of column. */ name?: string; /** Is the column a header? */ header?: boolean; +}; + +/** Schema for a text column in a table editor. + +Each cell editor is a text input field. + */ +export type TextColumnSchema = BaseColumnSchema & { + contentType: "string"; - /** Content of the column at the given row. */ + /** Text content of the column at the given row. */ content: (row: Row) => string; /** Is the text valid as content for the column at the given row? @@ -20,14 +32,57 @@ export type ColumnSchema = { */ validate?: (row: Row, text: string) => boolean; - /** Sets the content for the columns at the given row, if possible. + /** Sets the content for the column at the given row, if possible. Returns whether setting the content was successful. If not specified, the column is not editable. */ - setContent?: (row: Row, text: string) => boolean; + setContent?: (row: Row, content: string) => boolean; }; +/** Schema for a boolean column in a table editor. + +Each cell editor is a checkbox. + */ +export type BooleanColumnSchema = BaseColumnSchema & { + contentType: "boolean"; + + /** Content of the column at the given row, if defined. */ + content: (row: Row) => boolean | null; + + /** Sets the content for the column at the given row. + + If not specified, the column is not editable. + */ + setContent?: (row: Row, content: boolean) => void; +}; + +/** Schema for an enum column in a table editor. + +Each cell editor is a select box. The enum variants are assumed to be strings. + */ +export type EnumColumnSchema = BaseColumnSchema & { + contentType: "enum"; + + /** List of variants comprising the enum. */ + variants: (row: Row) => string[]; + + /** Content of the column at the given row, if defined. */ + content: (row: Row) => string | null; + + /** Sets the content for the column at the given row. + + If not specified, the column is not editable. + */ + setContent?: (row: Row, content: string | null) => void; +}; + +/** Schema for a column in a table editor. */ +export type ColumnSchema = + | TextColumnSchema + | BooleanColumnSchema + | EnumColumnSchema; + /** Create schema for a column with numerical (floating point) data. */ export const createNumericalColumn = (args: { name?: string; @@ -36,10 +91,11 @@ export const createNumericalColumn = (args: { default?: number; validate?: (row: Row, data: number) => boolean; setData?: (row: Row, data: number) => void; -}): ColumnSchema => ({ +}): TextColumnSchema => ({ + contentType: "string", name: args.name, header: args.header, - content: (row) => { + content(row) { let value = args.data(row); if (value === undefined) { value = args.default ?? 0; @@ -47,7 +103,7 @@ export const createNumericalColumn = (args: { } return value.toString(); }, - validate: (row, text) => { + validate(row, text) { const parsed = Number(text); return !Number.isNaN(parsed) && (args.validate?.(row, parsed) ?? true); }, @@ -112,17 +168,24 @@ function Cell(props: { row: Row; schema: ColumnSchema; }) { - const { row, schema } = destructure(props); return ( - - - + + + {(schema) => } + + + {(schema) => } + + + {(schema) => } + + ); } -function CellEditor(props: { +function TextCellEditor(props: { row: Row; - schema: ColumnSchema; + schema: TextColumnSchema; }) { const { row, schema } = destructure(props); @@ -138,17 +201,64 @@ function CellEditor(props: { const [isValid, setIsValid] = createSignal(true); createEffect(() => setIsValid(schema().validate?.(row(), text()) ?? true)); + return ( + + setText(evt.target.value)} + onChange={(evt) => applyText(evt.target.value)} + /> + + ); +} + +function BooleanCellEditor(props: { + row: Row; + schema: BooleanColumnSchema; +}) { + const { row, schema } = destructure(props); + return ( setText(evt.target.value)} - onChange={(evt) => applyText(evt.target.value)} + type="checkbox" + checked={schema().content(row()) ?? false} + disabled={schema().setContent === undefined} + onInput={(evt) => schema().setContent?.(row(), evt.currentTarget.checked)} /> ); } + +function EnumCellEditor(props: { + row: Row; + schema: EnumColumnSchema; +}) { + const { row, schema } = destructure(props); + + return ( + + ); +} diff --git a/packages/frontend/src/stdlib/analyses/decapodes.css b/packages/frontend/src/stdlib/analyses/decapodes.css new file mode 100644 index 00000000..8423fd79 --- /dev/null +++ b/packages/frontend/src/stdlib/analyses/decapodes.css @@ -0,0 +1,6 @@ +.decapodes-domain { + display: flex; + flex-direction: row; + gap: 1ex; + margin-bottom: 1ex; +} diff --git a/packages/frontend/src/stdlib/analyses/decapodes.tsx b/packages/frontend/src/stdlib/analyses/decapodes.tsx index b2ce598d..4d0dbda6 100644 --- a/packages/frontend/src/stdlib/analyses/decapodes.tsx +++ b/packages/frontend/src/stdlib/analyses/decapodes.tsx @@ -1,4 +1,5 @@ -import { Match, Switch, createMemo, createResource, onCleanup } from "solid-js"; +import type { IReplyErrorContent } from "@jupyterlab/services/lib/kernel/messages"; +import { For, Match, Show, Switch, createMemo, createResource, onCleanup } from "solid-js"; import { isMatching } from "ts-pattern"; import type { DiagramAnalysisProps } from "../../analysis"; @@ -11,19 +12,30 @@ import { Warning, createNumericalColumn, } from "../../components"; -import { type DiagramJudgment, type LiveDiagramDocument, fromCatlogDiagram } from "../../diagram"; +import { + type DiagramJudgment, + type DiagramObjectDecl, + type LiveDiagramDocument, + fromCatlogDiagram, +} from "../../diagram"; import type { ModelJudgment, MorphismDecl } from "../../model"; import type { DiagramAnalysisMeta } from "../../theory"; +import { uniqueIndexArray } from "../../util/indexing"; import { PDEPlot2D, type PDEPlotData2D } from "../../visualization"; import Loader from "lucide-solid/icons/loader"; import RotateCcw from "lucide-solid/icons/rotate-ccw"; import baseStyles from "./base_styles.module.css"; +import "./decapodes.css"; import "./simulation.css"; /** Configuration for a Decapodes analysis of a diagram. */ export type DecapodesContent = JupyterSettings & { + domain: string | null; + mesh: string | null; + initialConditions: Record; + plotVariables: Record; scalars: Record; }; @@ -48,6 +60,10 @@ export function configureDecapodes(options: { description, component: (props) => , initialContent: () => ({ + domain: null, + mesh: null, + initialConditions: {}, + plotVariables: {}, scalars: {}, }), }; @@ -56,6 +72,7 @@ export function configureDecapodes(options: { /** Analyze a DEC diagram by performing a simulation using Decapodes.jl. */ export function Decapodes(props: DiagramAnalysisProps) { + // Step 1: Start the Julia kernel. const [kernel, { refetch: restartKernel }] = createResource(async () => { const jupyter = await import("@jupyterlab/services"); @@ -67,34 +84,57 @@ export function Decapodes(props: DiagramAnalysisProps) { const kernelManager = new jupyter.KernelManager({ serverSettings }); const kernel = await kernelManager.startNew({ name: "julia-1.11" }); + return kernel; + }); + + onCleanup(() => kernel()?.shutdown()); + + // Step 2: Run initialization code in the kernel. + const startedKernel = () => (kernel.error ? undefined : kernel()); + + const [options] = createResource(startedKernel, async (kernel) => { + // Request that the kernel run code to initialize the service. const future = kernel.requestExecute({ code: initCode }); - const reply = await future.done; + // Look for simulation options as output from the kernel. + let options: SimulationOptions | undefined; + future.onIOPub = (msg) => { + if (msg.header.msg_type === "execute_result") { + const content = msg.content as JsonDataContent; + options = content["data"]?.["application/json"]; + } + }; + + const reply = await future.done; if (reply.content.status === "error") { await kernel.shutdown(); - throw new Error(reply.content.evalue); + throw new Error(formatError(reply.content)); } - - return kernel; + if (!options) { + throw new Error("Allowed options not received after initialization"); + } + return { + domains: uniqueIndexArray(options.domains, (domain) => domain.name), + }; }); - onCleanup(() => kernel()?.shutdown()); - - const maybeKernel = () => (kernel.error ? undefined : kernel()); + // Step 3: Run the simulation in the kernel! + const initedKernel = () => + kernel.error || options.error || options.loading ? undefined : kernel(); - const [result, { refetch: rerunSimulation }] = createResource(maybeKernel, async (kernel) => { + const [result, { refetch: rerunSimulation }] = createResource(initedKernel, async (kernel) => { // Construct the data to send to kernel. const simulationData = makeSimulationData(props.liveDiagram, props.content); if (!simulationData) { return undefined; } - + console.log(JSON.parse(JSON.stringify(simulationData))); // Request that the kernel run a simulation with the given data. const future = kernel.requestExecute({ code: makeSimulationCode(simulationData), }); - // Handle output from the kernel. + // Look for simulation results as output from the kernel. let result: PDEPlotData2D | undefined; future.onIOPub = (msg) => { if ( @@ -108,7 +148,7 @@ export function Decapodes(props: DiagramAnalysisProps) { const reply = await future.done; if (reply.content.status === "error") { - throw new Error(reply.content.evalue); + throw new Error(formatError(reply.content)); } if (!result) { throw new Error("Result not received from the simulator"); @@ -116,6 +156,10 @@ export function Decapodes(props: DiagramAnalysisProps) { return result; }); + const obDecls = createMemo(() => + props.liveDiagram.formalJudgments().filter((jgmt) => jgmt.tag === "object"), + ); + const scalarDecls = createMemo(() => { const liveModel = props.liveDiagram.liveModel; return liveModel.formalJudgments().filter((jgmt) => @@ -134,6 +178,7 @@ export function Decapodes(props: DiagramAnalysisProps) { const scalarSchema: ColumnSchema[] = [ { + contentType: "string", header: true, name: "Scalar constant", content: (mor) => mor.name, @@ -148,17 +193,54 @@ export function Decapodes(props: DiagramAnalysisProps) { }), ]; - const header = () => ( + const variableSchema: ColumnSchema[] = [ + { + contentType: "string", + header: true, + name: "Variable", + content: (ob) => ob.name, + }, + { + contentType: "enum", + name: "Initial/boundary", + variants() { + if (!props.content.domain) { + return []; + } + return options()?.domains.get(props.content.domain)?.initialConditions ?? []; + }, + content: (ob) => props.content.initialConditions[ob.id] ?? null, + setContent: (ob, value) => + props.changeContent((content) => { + if (value === null) { + delete content.initialConditions[ob.id]; + } else { + content.initialConditions[ob.id] = value; + } + }), + }, + { + contentType: "boolean", + name: "Plot", + content: (ob) => props.content.plotVariables[ob.id] ?? false, + setContent: (ob, value) => + props.changeContent((content) => { + content.plotVariables[ob.id] = value; + }), + }, + ]; + + const Header = () => (
Simulation - + - + ) {
); + const DomainConfig = (domains: Map) => ( +
+ Domain: + + + {(name) => ( + <> + Mesh: + + + )} + +
+ ); + return (
- + + {(options) => DomainConfig(options().domains)}
+
- {"Loading the AlgebraicJulia service..."} + + {"Loading the AlgebraicJulia service..."} + {(error) => ( - - {error().message} + + 
{error().message}
+ + )} + + + {(error) => ( + + 
{error().message}
 )} {"Running the simulation..."} - {(error) => {error().message}} + {(error) => ( + + 
{error().message}
+ + )} @@ -206,6 +346,10 @@ export function Decapodes(props: DiagramAnalysisProps) { ); } +const formatError = (content: IReplyErrorContent): string => + // Trackback list already includes `content.evalue`. + content.traceback.join("\n"); + /** JSON data returned from a Jupyter kernel. */ type JsonDataContent = { data?: { @@ -213,7 +357,25 @@ type JsonDataContent = { }; }; -/** Data send to the Julia kernel defining a simulation. */ +/** Options supported by Decapodes, defined by the Julia service. */ +type SimulationOptions = { + /** Geometric domains supported by Decapodes. */ + domains: Domain[]; +}; + +/** A geometric domain and its allow discretizations. */ +type Domain = { + /** Name of the domain. */ + name: string; + + /** Supported meshes that discretize the domain. */ + meshes: string[]; + + /** Initial/boundary conditions supported for the domain. */ + initialConditions: string[]; +}; + +/** Data sent to the Julia kernel defining a simulation. */ type SimulationData = { /** Judgments defining the diagram, including inferred ones. */ diagram: Array; @@ -221,7 +383,19 @@ type SimulationData = { /** Judgments defining the model. */ model: Array; - /** Mapping from IDs of scalar operations to numerical values. */ + /** The geometric domain to use for the simulation. */ + domain: string; + + /** The mesh to use for the simulation. */ + mesh: string; + + /** Mapping from variable UUIDs to enum variants for initital conditions. */ + initialConditions: Record; + + /** Variables to plot, a list of UUIDs. */ + plotVariables: Array; + + /** Mapping from UIIDs of scalar operations to numerical values. */ scalars: Record; }; @@ -231,19 +405,22 @@ import IJulia IJulia.register_jsonmime(MIME"application/json"()) using AlgebraicJuliaService + +JsonValue(supported_decapodes_geometries()) `; /** Julia code run to perform a simulation. */ -const makeSimulationCode = (data: SimulationData) => ` -system = System(raw"""${JSON.stringify(data)}"""); -simulator = evalsim(system.pode); +const makeSimulationCode = (data: SimulationData) => + ` + system = PodeSystem(raw"""${JSON.stringify(data)}"""); + simulator = evalsim(system.pode); -f = simulator(system.dualmesh, system.generate, DiagonalHodge()); + f = simulator(system.geometry.dualmesh, system.generate, DiagonalHodge()); -soln = run_sim(f, system.init, 100.0, ComponentArray(k=0.5,)); + soln = run_sim(f, system.init, 50.0, ComponentArray(k=0.5,)); -JsonValue(SimResult(soln, system)) -`; + JsonValue(SimResult(soln, system)) + `; /** Create data to send to the Julia kernel. */ const makeSimulationData = ( @@ -254,11 +431,21 @@ const makeSimulationData = ( if (validatedDiagram?.result.tag !== "Ok") { return undefined; } + + const { domain, mesh, initialConditions, plotVariables, scalars } = content; + if (domain === null || mesh === null || !Object.values(plotVariables).some((x) => x)) { + return undefined; + } + return { diagram: fromCatlogDiagram(validatedDiagram.diagram, (id) => liveDiagram.objectIndex().map.get(id), ), model: liveDiagram.liveModel.formalJudgments(), - scalars: content.scalars, + domain, + mesh, + initialConditions, + plotVariables: Object.keys(plotVariables).filter((v) => plotVariables[v]), + scalars, }; }; diff --git a/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx b/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx index d8ae180a..782e9f49 100644 --- a/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx +++ b/packages/frontend/src/stdlib/analyses/lotka_volterra.tsx @@ -70,6 +70,7 @@ export function LotkaVolterra( const obSchema: ColumnSchema[] = [ { + contentType: "string", header: true, content: (ob) => ob.name, }, @@ -94,6 +95,7 @@ export function LotkaVolterra( const morSchema: ColumnSchema[] = [ { + contentType: "string", header: true, content: (mor) => mor.name, }, diff --git a/packages/frontend/src/visualization/pde_plot.tsx b/packages/frontend/src/visualization/pde_plot.tsx index 71de6293..7ee24ebf 100644 --- a/packages/frontend/src/visualization/pde_plot.tsx +++ b/packages/frontend/src/visualization/pde_plot.tsx @@ -1,6 +1,7 @@ import { makeTimer } from "@solid-primitives/timer"; import type { EChartsOption } from "echarts"; import { createMemo, createSignal, lazy } from "solid-js"; +import invariant from "tiny-invariant"; const ECharts = lazy(() => import("./echarts")); @@ -16,7 +17,7 @@ export type PDEPlotData2D = { y: number[]; /** Values of the state variable over time. */ - state: StateVarAtTime[]; + state: Record; }; /** The data of a state variable at a given time. */ @@ -30,11 +31,25 @@ export function PDEPlot2D(props: { const min = (x: number, y: number) => Math.min(x, y); const max = (x: number, y: number) => Math.max(x, y); + const firstState = (): StateVarAtTime[] => { + // FIXME: Shouldn't just take the first one! + const keys = Object.keys(props.data.state); + invariant(keys.length === 1); + const key = keys[0]; + const state = key && props.data.state[key]; + invariant(state); + return state; + }; + const minValue = createMemo(() => - props.data.state.map((data) => data.map((triple) => triple[2]).reduce(min)).reduce(min), + firstState() + .map((data) => data.map((triple) => triple[2]).reduce(min)) + .reduce(min), ); const maxValue = createMemo(() => - props.data.state.map((data) => data.map((triple) => triple[2]).reduce(max)).reduce(max), + firstState() + .map((data) => data.map((triple) => triple[2]).reduce(max)) + .reduce(max), ); const [timeIndex, setTimeIndex] = createSignal(0); @@ -86,7 +101,7 @@ export function PDEPlot2D(props: { { name: "Value", type: "heatmap", - data: props.data.state[idx], + data: firstState()[idx], emphasis: { itemStyle: { borderColor: "black",