Merge pull request #283 from ToposInstitute/cm/navier-stokes

Improved Decapodes integration: domains, meshes, initial conditions, plot variables

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"

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

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

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))