Notes Level Set Distortion Please read Enright et al., Animation - PDF document

Notes Level Set Distortion � Please read Enright et al., “Animation and � Assuming even no numerical diffusion rendering of complex water surfaces”, problems in level set advection (e.g. well- SIGGRAPH’02 resolved on grid), level sets still have problems � Initially equal to signed distance � After non-rigid motion, stop being signed distance • E.g. points near interface will end up deep underwater, and vice versa cs533d-winter-2005 1 cs533d-winter-2005 2 Fixing Distortion Reinitialization � Idea: we have a distorted � , | �� | � 1 � Remember it’s only zero isocontour we care about - free to change values away � Want to return to | �� |=1 without disturbing the location of the interface from interface � If we’re not too far from | �� |=1, makes sense to � Can reinitialize to signed distance use an iterative method (“redistance”) • We can even think of each iteration as a pseudo-time • Without moving interface, change values to be step the signed distance to the interface • Information should flow outward from interface • Advection in direction sign( � )n and with rate of change sign( � ): � � � t + sign ( � ) � � � � � � � = sign ( � ) � � � � cs533d-winter-2005 3 cs533d-winter-2005 4 Reinitialization cont’d Discretization � Simplifying this we get: � When we discretize (e.g. with semi- Lagrangian) we’ll end up interpolating with ( ) = 0 � t + sign ( � ) � � � 1 values on either side of interface � This is another Hamilton-Jacobi � Limit the possibility for weird stuff to equation… happen, like � changing sign • If we want | �� |=1 to very high order accuracy, � So instead of sign( � ), use S( � 0 ) can use high-order HJ methods • Can never flip sign • Sign function smeared out to be smooth: � 0 S ( � 0 ) = 2 � x 2 + � � 0 2 ( ) � 0 cs533d-winter-2005 5 cs533d-winter-2005 6

Aside: initialization Problems � This works well if we’re already close to signed � Reinitialization will unfortunately slightly distance move the interface (less than a grid cell) � What if we start from scratch at t=0? � Errors look like, as usual, extra diffusion or • For very simple geometry, may construct � smoothing analytically • In addition to the errors we’re making in • More generally, need to numerically approximate advection… � One solution - if we can at least get inside/outside on the grid, can run reinitialization equation from there (1st order accurate) cs533d-winter-2005 7 cs533d-winter-2005 8 Fast methods Velocity extrapolation � Problem with reinitialization from scratch - to get full field, � We can exploit level set to extrapolate velocity need to take O(n) steps, each costs O(n 3 ) field outside water � Can speed up with local level set method • Not a big deal for pressure solve - can directly handle • Only care about signed distance near interface, so only compute extrapolation there those O(n 2 ) values in O(1) steps • But a big deal for advection - with semi-Lagrangian • Gives optimal O(n 2 ) complexity (but the constant might be big!) method might be skipping over, say, 5 grid cells � If we really want full grid, but fast: • So might want velocity 5 grid cells outside of water • Fast Marching Method O(n 3 log n) (Tsitsiklis, Sethian) � Simply take the velocity at an exterior grid point • Fast Sweeping Method O(n 3 ) (Zhao) • Other more geometric ideas (e.g. Tsai, Mauch) to be interpolated velocity at closest point on � Nice property: more careful about not letting the interface interface move • Alternatively, propagate outward to solve � u � � � = 0 similar to redistancing methods cs533d-winter-2005 9 cs533d-winter-2005 10 Particle-Level Set Level set advancement � Last time - mentioned marker particles (MAC) � Put marker particles with values of � attached in method great for rough surfaces a band near the surface • We’re also storing � on the grid, so we don’t need � But if we want surface tension (which is particles deep in the water strongest for rough flows!) or smooth water • For better results, also put particles with � >0 (“air” surfaces, we need a better technique particles) on the other side � Hybrid method: particle-level set � After doing a step on the grid and moving � , also • [Fedkiw and Foster], [Enright et al.] move particles with (extrapolated) velocity field • Level set gives great smooth surface - excellent for � Then correct the grid � with the particle � getting mean curvature • Particles correct for level set mass (non-)conservation � Then adjust the particle � from the grid � through excessive numerical diffusion cs533d-winter-2005 11 cs533d-winter-2005 12

Level set correction Artificial Compressibility � Look for escaped particles � Let’s make a quick detour… • Any particle on the wrong side (sign differs) by more than the � So far we’ve seen projection methods for particle radius | � | enforcing divergence-free constraint � Rebuild � <0 and � >0 values from escaped particles • Means solving Poisson equation for pressure (taking min/max’s of local spheres) • Big, sparse linear system - it’s slow, it’s the bottleneck � Merge rebuilt � <0 and � >0 by taking minimum- • Particularly on parallel architectures - global magnitude values � Reinitialize new grid � communication • Needs a weird staggered grid, or more complicated � Correct again problems and fixes � Adjust particle � values within limits � An alternative: artificial compressibility (never flip sign) cs533d-winter-2005 13 cs533d-winter-2005 14 Revisiting incompressibility Equation of state � Turn hard constraint � • u=0 into soft constraint � Real fluids are not incompressible � We just make the idealization of incompressibility • Allow the fluid to compress a little, but add restoring • Water, air are very close unless material velocity comparable to force to get it back sound speed (transonic or faster) � Real compressible flow does this, but with all • Simplifies math a lot sorts of complications from thermodynamics • Means we can ignore sound waves in numerical methods - terrible time step limit � We’ll fake it, simplify compressible flow � But we could go the other way • We don’t care about compressibility effects and • Replace real compressible physics with fake ones that still have ideally won’t even see them at all sound speed much faster than material velocity � Artificial equation of state: p=c 2 � � [Chorin ‘67] cs533d-winter-2005 15 cs533d-winter-2005 16 New equations What is c? � Need to include density again (continuity � Can derive acoustic wave equation eq. = conservation of mass) � We want to make sure that the maximum ( ) = 0 material speed (u) is much less than c � t + � � � u • E.g. choose c at least 10 |u| max � t + u � � � = � � � � u � Note that time step limit (for explicit � And momentum equation methods) will have � t< � x/c ( ) � � � µ � u + � u T u t + u � � u + 1 � � p = g + 1 • Hope is that 10+ times the number of steps is � And the new equation of state worth it for no pressure solve, easier p = c 2 � programming, etc. cs533d-winter-2005 17 cs533d-winter-2005 18

The flies in the ointment Mesh Free Methods � To make it stable without a staggered grid, need artificial viscosity, or sophisticated conservation law methods • Just like shallow water � We may have to give up a lot of space and time resolution to make it work cs533d-winter-2005 19 cs533d-winter-2005 20 Particle fluids Mesh-free? � Particles are great for advection (hence marker � Mathematically, particle-only methods are particles in MAC, particle-level set, etc.) independent of meshes � So get rid of the mesh - figure out how to do � p � Practically, need an acceleration structure to etc. with just the particles speed up finding neighbouring particles (to figure out forces) � Basic qualitative behaviour of fluids: resist density changes � Most popular structure (for non-adaptive codes, • When particles get too close, add repulsion forces between them i.e. where h=constant for all particles) is… • When they get just a little too far, add attraction forces a mesh (background grid) • When far, no force at all � Damp particle interactions • Otherwise we see small-scale vibration (“heat”) • Also accounts for viscosity cs533d-winter-2005 21 cs533d-winter-2005 22 SPH Kernel � Smoothed Particle Hydrodynamics � Need to define particle’s influence in surrounding � SPH can be interpreted as a particular way of choosing space (how we’ll build the basis functions) forces, so that you converge to solving Navier-Stokes � Choose a kernel function W � [Lucy’77], [Gingold & Monaghan ‘77], [Monaghan…], • Smoothed approximation to � [Morris, Fox, Zhu ‘97], … • W(x)=W(|x|) - radially symmetric � Similar to FEM, we go to a finite dimensional space of • Integral is 1 functions • W=0 far enough away - when |x|>2.5h for example • Basis functions now based on particles instead of grid elements • Can take derivatives etc. by differentiating the real function from � Examples: the finite-dimensional space • Truncated Gaussian • Splines (cubic, quartic, quintic, …) cs533d-winter-2005 23 cs533d-winter-2005 24