Multigrid

[KevinDCarlson]

I've tried to separate out the multigrid-specific work from the last couple of months from the work applying to simplicial complexes and geometric maps more generally which is covered in #101 . You might diff this with multiscale rather than main for easiest review.

The main output and good entry point so far is the stokes.md notebook.

[lukem12345]

I can do the algebra by hand with L1(x)=b on some analytically-validated solution as a diagnostic

[KevinDCarlson]

TODO: See whether today's 1-form prolongation and restriction functions work when you do give multigrid invertible operators, even if they're not physically interesting.

[lukem12345]

Decapodes related todos:

• Emit LinearOperators.jl-compatible functions by re-interpreting TVars

[KevinDCarlson]

TODO: See whether today's 1-form prolongation and restriction functions work when you do give multigrid invertible operators, even if they're not physically interesting.

(They do not.)

[lukem12345]

@KevinDCarlson Is the trie backend necessary for the multigrid feature?

[KevinDCarlson]

@lukem12345 It's not strictly necessary but GeometricMaps go through it and bypassing that would make construction of what would be just implicitly GeometricMaps less safe.

[lukem12345]

Does the notion of GeometricMap as given in this PR require a trie/ simplex tree backend? I'm just looking at how to modularize the code in this PR (after merging).

[KevinDCarlson]

Yes, essentially; a geometric map is fully defined by its action on 0-simplices but that means it's only well defined on simplicial complexes.

[lukem12345]

Gotcha. The requirement is that it is capable of being faithfully represented as a trie/ simplex tree, which are guaranteed to only represent simplicial complexes. This goes for the differential operators in this package too (assuming they are going to be used for DEC purposes), and further that they are manifold-like. Those functions typically just assume the simplicial complex requirements are met. Aside, I added a helper function for one of the manifold-like conditions that seemed more likely to occur in practice.

Is there a test for the case in which a SimplicialSet that is not a SimplicialComplex has the SimplicialComplex constructor called on it? What would occur?

[KevinDCarlson]

I believe it'll error with something about trying to duplicate a simplex.

[lukem12345]

Awesome. It sounds like a method that does is_a_simplicial_complex should be super quick to define. Although surely there's another way to check this property in constant memory.

[KevinDCarlson]

@lukem12345 I'll factor the simplicial complexes stuff out of this PR then.

[jpfairbanks]

It is time to start getting this code reviewed and merged. Lets factor the SimplicialComplex / Trie / Geometric Map stuff into one PR and then do the multigrid stuff in a second PR. The first PR should now be stable since we successfully used those APIs for the multigrid effort.

[jpfairbanks]

I think I resolved the merge conflicts correctly. The \flat_mat code was the part I am least confident in.

[KevinDCarlson]

@jpfairbanks The simplicial complexes etc are supposed to be covered in #101 , if we want to bring that one back up for review.

[lukem12345]

I think I resolved the merge conflicts correctly. The \flat_mat code was the part I am least confident in.

The flat operators from main are the ones that need to be used, for whoever rebases next.

[jpfairbanks]

I kept both versions since they had different signatures. This branch currently has three versions.

function ♭(s::AbstractDeltaDualComplex2D, X::AbstractVector, ::PPFlat)
  map(edges(s)) do e
  # Assume linear-interpolation of the vector field across the edge,
 # determined solely by the values of the vector-field at the endpoints.
    vs = edge_vertices(s,e)
    l_vec = 1/2*(X[vs][1]+X[vs][2])
    e_vec = (point(s, tgt(s,e)) - point(s, src(s,e))) * sign(1,s,e)
    dot(l_vec, e_vec)
 end
end

function ♭_mat(s::AbstractDeltaDualComplex2D)
  ♭_mat(s, ∂(2,s))
end

♭_mat(s::AbstractDeltaDualComplex2D, f::DPPFlat) =
  ♭_mat(s, ∂(2,s), f)

Should we delete any of them?

[lukem12345]

New type signatures were added to the old DPPFlat since it was previously the only DiscreteFlat. The middle one that you quoted here should get axed.

[jpfairbanks]

@KevinDCarlson I just did that directly on this branch.

[KevinDCarlson]

Something is wrong with loading the torus mesh (which I've commented out) and maybe with BuildKite, but otherwise everything seems to be working here, would be happy for somebody to start reviewing.

[lukem12345]

Likely guess is that it is saved as a JSON ACSet, and this PR renamed some Homs.

[KevinDCarlson]

Likely guess is that it is saved as a JSON ACSet, and this PR renamed some Homs.

Ahhh, nice. The D vs dual thing, maybe. I've never used the JSON ACSet API so I'd be much obliged if somebody else wants to fix it, but I can probably figure it out eventually if need be.

[KevinDCarlson]

Ah, now I remember; the renaming makes it much easier and more resilient to write automatic migrations from a dual complex to the schema for its associated primal type. Anyway, it's reverted since we weren't really using extract_dual for anything yet.

[lukem12345]

Calling this either binary or medial subdivision would work