paper reference discrete lie derivative?

[spinynormal]

Hi eelco,
it seams the advection technique wich is used here is differently from Stable, Circulation-Preserving, Simplicial Fluids?

[EelcoHoogendoorn]

Yeah; while it definitely draws its inspiration from there, there are two differences:

• one small difference: incompressibility can be solved by either pressure projection or divergence constraint; should give the same result up to numerical error; but the practicalities of it become quite different when boundary conditions are involved. What I would like to do here instead is implement a full time-dependent stokes equation instead, since it would allow for the same unified treatment of boudary conditions as in the steady-state stokes example

• the bigger difference: the suggested dual barycentric interpolator (warren 2007) is difficult to implement efficiently in python. I therefore came up with my own interpolation scheme for interpolating over dual polyhedra. The implementation of this starts at complex.simplicial.sample_dual_0; and much of the crux lies in complex.weighted_average_operators, and complex.simplicial.pick_fundamental. The resulting method is efficient and elegantly generalizes to nd-power-complexes and cubical complexes. Which is all great on paper; the interpolator has all the desired theoretical properties*, can be elegantly implemented, and it works exactly as I want it to subjectively (see examples/tracing/pick_sphere_2.py for illustration); but the actual flow simulation so far only works correctly for cubical complexes. For triangular complexes, distortions appear in the flowfield. This might be due to some under-appreciated properties of the warren2007 interpolator employed elcot2007; or a bug in my advection code. I havnt had the time recently to get to the bottom of it. Trying to sort this out is why the Advector class is such a mess of commented out code.

If you do know of an elegant/efficient method to implement the warren2007-type interpolator in python/numpy, let me know! It would be great to have a baseline to compare to, so I can see where the appoaches diverge in practice.

*that is, aside from the standard barycentric properties:

• piecewise linear over fundamental domains (term im using in my code for the intersection of an overlapping primal and dual cell)
• interpolation over any lower-dimensional sub-complex of the interpolation domain, is only a function of the geometry participating in that sub-complex; all other coefficients go to zero; that is, interpolation reduces to triangular barycentric on the face of a tet, and so on.

[spinynormal]

in the last chapter its mentioned to project the interpolated vectorfield on curvature normals for continuous vectorfield, https://github.com/bertholet/DEC_2013/blob/master/DEC%20in%20a%20Week.pdf

dunno regarding py/numpy stuff, i am trying to avoid the computation of the flowlines, is the dual_flux_to_dual_velocity your owen discretization in the advection.py, i cant relate it,

[EelcoHoogendoorn]

dual_flux_to_dual_velocity maps from a dual-1-form-tangential-flux description of the flowfield, to an nd-vector-per-dual-vertex description of the flowfield. This step is identical to elcot2007 and the chapter you link to (nice resource, btw). "The authors of [ETK+07] propose the following interpolation scheme to get a continuous vector field. First, velocity vectors are computed at the positions of dual vertices; this can be done using the sharp operator from 6.4.1.".

Not sure they really explain it there, but the idea is that for an incompressible field, there exists a unique vector at each dual vertex of an n-simplex which exactly projects onto the tangent flux along each of its n +1 incident dual edges. The multiplication with pinv is supposed to find that unique vector. If you uncomment line 88-89 in advection.py, it is checked at runtime that these relations between the vector field and the tangent flux 1-form are indeed satisfied; the assert passed for all circumstances that I tried.

Good point about the reprojection; I looked into that but setting complex_type = 'simplex' in euler_flow.py gives the same problems, and that is a flat triangle mesh scenario.

Not sure what you mean by avoiding the computation of the flowlines. It can be avoided in the sense of doing a mathematically equivalent pressure projection step instead (as per VorticityAdvector.pressure_projection); but you are still solving a large sparse system, and the divergence constraint method tends to be numerically better conditioned in my experience. Still, I think it is more elegant since you do not have to deal with the harmonic component seperately. But if you want to simulate incompressible flows there isnt really any way around solving a coupled system of equations somehow, since molecules being strongly coupled is kinda the essence of incompressibility.

[EelcoHoogendoorn]

i cant relate it,

With respect to that; I have made some effort to keep this code well documented and clean, but it isnt to the same standard everywhere yet. So your feedback with regard to what is and isnt clear is most welcome.

Note that I have set myself quite ambitious goals with this project in terms of generality; for example, the euler flow code is intended to work exactly the same, regardless of if you work with a triangle mesh, a tet mesh, or a cubical mesh. This naturally leads to some divergence in coding style and naming relative to what is common in most DEC literature; like the product with the pseudo-inverse of the matrix of dual edge vectors to transform from dual flux to velocity vector, for instance. It isnt trivial to see why that should be the same as the way this is typically done in the DEC literature, and this is not the only place where some more documentation certainly wouldnt hurt.

[spinynormal]

the hint for the sharp operator is really helpfully !
i meant the pathtrace for the dualvertices in the semi langrangian procedure, i read somewhere it is numerically unstable for triangle meshes, so i was pretty happy to see a different way for comparison purposes later on, i had to shot some questions =) but i really do appreciate your answers and thanks for clarifying

[EelcoHoogendoorn]

Do you have a reference for that; because that is indeed exactly my experience; that is isnt stable for triangle meshes; but works fine for quads. I do not recall anyone mentioning that in print though; and itd be a pretty big omission from elcot2007 if their implementation did share the same quirk. It could save me a lot of debugging time if this was actually expected behavior rather than a consequence of my interpolation or somesuch...

[spinynormal]

using the vorticity formulation of incompressible fluid flow. Their method
preserves circulation, while ours preserves vorticity. We improve
on their work by avoiding the computation of flow
lines of the velocity vector field, which tends to be both challenging
to implement and numerically unstable for triangle meshes
http://www.cs.technion.ac.il/~mirela/publications/ff-surf.pdf,

in the thesis from elcott is some talk about the appearance of discontinuous in the advection step _3.2.3 , 3.2.4 but i assume your aware of ?
http://www.multires.caltech.edu/pubs/SharifMSThesisDiscreteFluids.pdf,

[EelcoHoogendoorn]

ah! yeah ive read elcots thesis many times, but it was along time ago. Indeed he does mention stability. Before I started implementing this myself I didnt realize the extent to which he seems to be relying on mesh quality to get acceptable results; if you look at his triangle meshes, they are extraordinarily uniform. Thanks for the refresher! Seems like I can stop hunting for bugs in that aspect of my code.

That functional fluids paper looks extremely interesting; hope to take a closer look at it tonight. Now i get what you mean by avoiding the computation of the streamfunction.

[EelcoHoogendoorn]

Interesting, yes; but id rate their documentation worse than my own :P.

From their own reference it seem pycomplex already has all the low level operators to construct their Dv. That is, eq 11 from http://www.cs.technion.ac.il/~mirela/publications/fvf.pdf

There are many properties of their scheme I do not yet fully appreciate though. It seems to me that this can be viewed as a DEC-type construction of a classic antisymmetric finite difference local advection operator, having many of the same drawbacks that come with it, compared to a semi-lagrangian method. Yet they mention insensitivity to timestep; but how would they overcome the CFL condition, using only a local operator without any path-tracing?

[EelcoHoogendoorn]

the flux-to-vector code was infact quite badly broken; ive fixed that on current master and added a test for it that visualizes the current behavior; which i am happy with.

Still, the flow on triangle meshes isnt exactly as id like it to be. You can see it better in an animation, but here is a still from running the euler_flow.py example; and you can see lines in the vorticity field imprinted from the underlying icosahedral structure; which is a bug, not a feature. Maybe its something I can fix but I am starting to think it is something intrinsic to elcot2007 on structured meshes.

If you run this on an unstructured mesh it should serve as a quite solid comparison of any attempted implementation of that functional fluids paper. For a sphere with 10k vertices and BFECC enabled a single update step takes less than a second, and its also second order in time and space, with good conservation properties. So quite apples-to-apples.

[spinynormal]

i run the simulation it looks really good so far! thanks for debugging=)

can i might make a future request for adding interior boundaries on the grid, i started with the "warren interpolator" so there isnt really a rush but i would love to compare it!

i am not really aware of his stuff, but desbrun also uses a functional map similiar to azencots in: https://authors.library.caltech.edu/62269/2/LMHTD15.pdf

[EelcoHoogendoorn]

Happy to help; not sure what you mean by internal boundaries precisely; the euler_flow code currently has some limitations in the way it deals with boundary conditions, in the sense that velocity=0 is assumed at the boundary everywhere. But it does work fine with boundaries, internal or otherwise. I do plan on adding more flexible boundary handling eventually; my approach actually lends itself to that very well; and I never managed to figure out how this is supposed to work with the Warren method.

I just pushed some performance optimizations; the code is pretty close to optimal now performance wise I think. Every update step only requires applying a bunch of (sparse) linear operators, and a single KDTree lookup per advected vertex. The only main optimization left would be to add some caching to complex.sample_dual_0; most advected vertices map to the same triangle from timestep to timestep so the KDTree lookup only needs to happen for a fraction of the vertices. Havnt profiled it yet but with that optimization id be surprised if it was possible to do semi-lagrangian fluids on unstructured meshes on the cpu much faster.

Thats another very interesting reference! Its been more than 10 years since keeping up with this field was a fulltime occupation of mine, so im a little out of the loop, but its nice to see that cool things are still happening. If I had the time, id implement all these kind of papers. As well written as many of the papers in this community are, actually providing something resembling runnable code hasnt been a thing since Stam1999 I feel. If you are at all into python/numpy dont hestitate to make any pullrequests on this repo; id be happy to review.