Move sign inside wedge pd, and solve for pressure in euler

[lukem12345]

Previously, we were multiplying by the sign of a dual 2-form after performing the primal-dual 1-1 wedge product at the outermost "test level." It makes sense to curry the sign-accommodation code inside these wedge product operators themselves.

This PR contains two further improvements to the tests: solving for pressure in the test for Euler's equations, and a standalone interior product dual-dual test set.

[GeorgeR227]

Since you're already making changes, this PR should also address issue #76. At the very least moving the only.() call outside the returned function.

[lukem12345]

@jpfairbanks Can you give this a review this afternoon or Monday? This branch is used for the current NS simulations under investigation

[jpfairbanks]

@lukem12345, we had a race condition, I noticed that you have an -a - b instead of -(a+b).

[lukem12345]

I will run some @btimes and report back if I find an improvement!

[lukem12345]

I did not find a difference in runtime or allocs.

julia> using CombinatorialSpaces, GeometryBasics, BenchmarkTools;                                                                                                   
                                                                                                                                                                    
julia> Point3D = Point3{Float64};                                                                                                                                   
                                                                                                                                                                    
julia> s = loadmesh(Icosphere(5));                                                                                                                                  
                                                                                                                                                                    
julia> sd = EmbeddedDeltaDualComplex2D{Bool, Float64, Point3D}(s);                                                                                                  
                                                                                                                                                                    
julia> subdivide_duals!(sd, Barycenter());                                                                                                                          
                                                                                                                                                                    
julia> ldd = ℒ_dd(Tuple{1,1}, sd);                                                                                                                                  
                                                                                                                                                                    
julia> function myℒ_dd(::Type{Tuple{1,1}}, s::SimplicialSets.HasDeltaSet)                                                                                           
         # ℒ := -diuv - iduv                                                                                                                                        
         d0 = dec_dual_derivative(0, s)                                                                                                                             
         d1 = dec_dual_derivative(1, s)                                                                                                                             
         i1 = interior_product_dd(Tuple{1,1}, s)                                                                                                                    
         i2 = interior_product_dd(Tuple{1,2}, s)                                                                                                                    
                                                                                                                                                                    
         (f,g) ->                                                                                                                                                   
           -(d0 * i1(f,g) +                                                                                                                                         
             i2(f,d1 * g))                                                                                                                                          
       end;                                                                                                                                                         
                                                                                                                                                                    
julia> mdd = myℒ_dd(Tuple{1,1}, sd);                                                                                                                                
                                                                                                                                                                    
julia> u = ones(ne(sd));                                                                                                                                            
                                                                                                                                                                    
julia> all(mdd(u,u) .== ldd(u,u))                                                                                                                                   
true                                                                                                                                                                
                                                                                                                                                                    
julia> @btime ldd(u,u);                                                                                                                                             
  820.285 μs (52 allocations: 808.38 KiB)                                                                                                                           
                                                                                                                                                                    
julia> @btime mdd(u,u);                                                                                                                                             
  822.990 μs (52 allocations: 808.38 KiB)

The types given by @code_warntype is similar for both, just reordered:

f::Vector{Float64}                                                                                                                                                                                 
  g::Vector{Float64}                                                                                                                                                                                 
Body::Any                                                                                                                                                                                            
1 ─ %1  = Core.getfield(#self#, :d0)::SparseArrays.SparseMatrixCSC{Int8, Int32}                                                                                                                      
│   %2  = Core.getfield(#self#, :i1)::CombinatorialSpaces.FastDEC.var"#32#33"{LinearAlgebra.Diago                                                                                                    
ol, Float64, Point{3, Float64}}, SparseArrays.SparseMatrixCSC{Float64, Int64}, CombinatorialSpace                                                                                                    
astDEC.var"#30#31"{SparseArrays.UMFPACK.UmfpackLU{Float64, Int32}}}                                                                                                                                  
│   %3  = (%2)(f, g)::Any                                                                                                                                                                            
│   %4  = (%1 * %3)::Any                                                                                                                                                                             
│   %5  = Core.getfield(#self#, :i2)::CombinatorialSpaces.FastDEC.var"#34#35"{SparseArrays.Sparse                                                                                                    
C{Float64, Int32}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}}                                                                                                                                 
│   %6  = Core.getfield(#self#, :d1)::SparseArrays.SparseMatrixCSC{Int8, Int32}                                                                                                                      
│   %7  = (%6 * g)::Vector{Float64}                                                                                                                                                                  
│   %8  = (%5)(f, %7)::Vector{Float64}                                                                                                                                                               
│   %9  = (%4 + %8)::Any                                                                                                                                                                             
│   %10 = -%9::Any                                                                                                                                                                                   
└──       return %10

[lukem12345]

The i1 function returns Any...

[lukem12345]

Tracing back through the function calls, the sign function, which returns a vector of positive and negative ones corresponding to the sign of a simplex, is returning ::Any. This is inside of the primal-dual 1,1 wedge product:

Body::Any                                                                                                                                                                                            
1 ─ %1 = CombinatorialSpaces.FastDEC.:*::Core.Const(*)                                                                                                                                               
│   %2 = Core.getfield(#self#, :sd)::EmbeddedDeltaDualComplex2D{Bool, Float64, Point{3, Float64}}                                                                                                    
│   %3 = CombinatorialSpaces.FastDEC.sign(2, %2)::Any                                                                                                                                                
│   %4 = Core.getfield(#self#, :Λ_cached)::CombinatorialSpaces.FastDEC.var"#7#8"{Tuple{Matrix{Int3                                                                                                   
2}, Matrix{Float64}, UnitRange{Int64}}}                                                                                                                                                              
│   %5 = Core.getfield(#self#, :♭♯_m)::SparseArrays.SparseMatrixCSC{Float64, Int64}                                                                                                                  
│   %6 = (%5 * g)::Vector{Float64}                                                                                                                                                                   
│   %7 = (%4)(f, %6)::Vector{Float64}                                                                                                                                                                
│   %8 = Base.broadcasted(%1, %3, %7)::Any                                                                                                                                                           
│   %9 = Base.materialize(%8)::Any                                                                                                                                                                   
└──      return %9

[lukem12345]

Tracing back through the function calls, the sign function, which returns a vector of positive and negative ones corresponding to the sign of a simplex, is returning ::Any.

This is because the type of the orientation is encoded as a type parameter to EmbeddedDeltaDualComplex2D.

One way to work around this is to hard-code Bool a la

CombinatorialSpaces.jl/src/DiscreteExteriorCalculus.jl

Line 1101 in ab1ea66

lazy_edge_orient .= one(Bool)