Skip to content

Instantly share code, notes, and snippets.

@michaltakac

michaltakac/cfd.md

Forked from Wumpf/cfd.md
Created November 19, 2020 22:02
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save michaltakac/ea0c34e641ebd18e0ddfa47970910ca0 to your computer and use it in GitHub Desktop.
Save michaltakac/ea0c34e641ebd18e0ddfa47970910ca0 to your computer and use it in GitHub Desktop.
Download ZIP
Notes on CFD
Raw
cfd.md

Computational Fluid Dynamics (in graphics mostly)

Disclaimer: I have no professional education in Compulational Fluid Dynamics, these are various notes on what I found/learned when trying to learn a bit on the topic. Some of the notes lack sources as I wrote them up later often just remembering some StackExchange answer. Some things I might have misinterpreted.

General / Overview

Overarching resources

General Concepts & Keywords

  • Density ρ: Mass per volume, m / V
  • Pressure p: Force per area, F / A
  • Viscosity: The more a fluid has the "thicker" it is
    • Dynamic viscosity μ: Relation of viscous stresses to rate of change
    • Kinematic viscosity ν: μ / ρ
  • Reynolds Number Re: A measure for how turbulent a flow is (higher is more turbulent)
  • Courant–Friedrichs–Lewy condition (CFL): Condition for convergence of a numerical solution for a partial differential equation. Wiki
  • Poisson Equation: Differential equation on laplacian, ∇²p=f

Overview of types of simulations

Solving Navier Stokes: Either Lagrangian or Eulerian point of view (wiki). Eulerian: Grid. Lagrangian: Particles

  • Lagragian:
    • Smoothed Particle Hydrodynamics (SPH): Particles moving around
  • Eulerian:
    • Finite volume method (FVM): Eularian
    • Finite element method (FEM): Eularian
    • FLIP / PIC (particle in cell): Particles but solve with eulerian grid
  • Other:
    • Lattice Boltzmann methods (LBM) do not operate on Navier-Stokes directly but on grid of particle distributions governed by Boltzmann equations

Pro/Con (pure) Lagragian:

  • advection is fairly trivial
  • incompressibility harder to enforce. Can have issues with nonuniform particle distribution

Pro/Con (pure) Eulerian:

  • incompressibility condition relatively easy to enforce
  • advection harder. Tecniques like vorticity confinement used to counter dissipation

SPH (Lagragian)

Resources

Basic idea: Reconstruct continous signal with gaussians (== smoothed particles). This makes every quantity in every point a weighted sum of particles/kernel functions. Generally, first need to estimate the density in every particle. Then, compute total change of velocity caused by pressure gradient, viscosity and direct forces.

Kernel function used to be typically a cubic spline. But there are others (PySPH implements quite a bunch) Many simulations use different kernel functions for different parts of the velocity update for more stability.

Muller et al. 2003, Particle-Based Fluid Simulation for Interactive Applications seemingly introduced computer graphics community to the concept of SPH. Kernel functions defined in there (poly6, spiky, viscosity) are widespread.

Wikipedia write about accuracy for fluid mechanics:

One drawback over grid-based techniques is the need for large numbers of particles to produce simulations of equivalent resolution. [...] However, accuracy can be significantly higher with sophisticated grid-based techniques, especially those coupled with particle methods (such as particle level sets), since it is easier to enforce the incompressibility condition in these systems. SPH for fluid simulation is being used increasingly in real-time animation and games where accuracy is not as critical as interactivity.

Pressure solver

Incompressible fluid keeps density constant <=> velocity field divergence free.
(compressible fluid allows for a certain amount of density change, so problem still similar)

Methods can generally be distinguished whether they ...

Mentioning a few selected approaches selected here:

Weakly Compressible SPH (WCPSH)

Solving with the standard approach, weakly compressible SPH (WCSPH), requires to allow some compressiblity, usually done by reducing speed of sound (-> lower pressure than there should actually be since this governs how pressure is derived from density), otherwise required timesteps get too small. Guidance on how to choose speed of sound and timestep slide 35 here

Modern paper typically cite Becker & Teschner 2007, Weakly compressible SPH for free surface flows which presents treatments on a full framework (including surface tension). See equation (9) and following on timestamp / speed of sound / compressibility relations.

Position Based Fluids (PBF)

Maklin & Muller 2013, Position Based Fluids
SPH technique with constraints as in a position based physics system (i.e. more or less all game physics - [citation needed]). Uses constraints between particles to enforce incompressibility. Constraints are solved iteratively.

Divergence-Free Smoothed Particle Hydrodynamics (DFSPH)

Bender & Koschier 2015, Divergence-Free Smoothed Particle Hydrodynamicss
Bender & Koschier 2017, Divergence-Free SPH for Incompressible and Viscous Fluids

Does not only focus on constant density but also on divergence free velocity field. CFL condition for maximum timestep applies.

Viscosity

Physical viscosity formulations typically take dynamic viscosity μ parameter. Muller et al. 2003, Particle-Based Fluid Simulation for Interactive Applications formulates physical viscosity as: fi = μ sum_over_particles_j(m * velocity_difference / density_j laplacian_of_kernel_func)

Artificial approaches are used to dampen unphysical oscilations on shocks (or generally high Reynolds numbers). Mentioning a few selected approaches selected here:

Artificial Viscosity

Term refers in general to any non-physical viscosity formulation, usually refers though to a special term added in the computation for pressure related accelleration. Typical formulations set artificial viscosity to zero for particles that do not move towards each other. See also D.J. Price PHD thesis, 3.5 Shocks

XSPH

By some classified as a kind of artificial viscosity.

Added on velocity integration. Instead of changing the position via r += v * dt, neighbors are taken into account using the SPH Kernel: r_a += (v + eps * ForAllParticlesB(m_b * (v_ba / ((ρ_a + ρ_b) * 0.5)) W_ab) )*dt (see 1992 review by Monaghan)

Can replace physical viscosity treatment completely:

The goal of artificial viscosity is just to damp out non-physical oscillations that plague the undamped method: a full treatment of physical viscosity, with a nonphysical coefficient tuned to stabilize the discretization, is overkill." Ghost SPH for Animating Water, Schechter et al. (this paper also uses the XSPH viscosity on the velocity, which makes it easier to stick to the usual framework - SPlisHSPlasH impl takes this further and puts it onto the accelleration by dividing by dt)

Neighborhood search

Need to access particles in neighborhood a lot (bound by smoothing distance). All approaches are based on some kind of (sparse) grid binning (at the very least conceptually!). A few general ideas:

  • spatial consistency
    • all algorithms profit from spatial coherency. employ space filling curve!
    • sorting doesn't need to be done every frame
      • can be costly since there might be many attributes
  • spatial hashing
  • sorting
    • sort by bin and remember offsets per bin
    • nice for GPUs since everything is in place
    • tends to rewrite all the data every frame though
    • GPU Tech Conference slides on solving this with counting sort on gpu
  • cell linked lists (CLL approach)
    • every particle has a linked list to its "true" neighbors
      • during query no need to look at cell data structure, but for builtup cells are still important
    • if particles are spatially (i.e. cell index) sorted then this linked list can be heavily compressed!
    • General approach:
      • sort all particles by grid cell index
        • keep a list of the first particle of all grid cells
      • for each particle create linked list, sped up by grid cell search
      • compress linked list and store it

Foam

Unified Spray, Foam and Bubbles for Particle-Based Fluids, Ihmsen et al 2012 (also talks a bit about rendering)

Hybrid methods (mostly Eulerian)

Approaches that combine eularian grid and lagrangian particles: Use particles to resolve transport, use grid to compute forces and collision response. (i.e. "fix" the neighborhood problem) Generally, transfer between representation creates error!

Important Methods:

  • PIC: Particle in Cell.
    • The "original"
    • keep quantities like velocity completely in the grid
    • leads to very strong damping effect, angular momentum not preserved
  • FLIP: FLuid Implicit Particle. Brackbill et. al. 1986/88. Widely cited is Zhu & Bridson 2005, Animating Sand as a Fluid
    • Very popular in visual effects - primary solution for fluids in Houdini, very popular Blender plugin exists
    • Less dissipative than PIC but occasionally less stable
    • transfer increments of velocity/displacement from grid to particles instead of interpolating from grid
    • suffers from instabilities, preserves some angular momentum, sometimes mixed with PIC for more stability (which brings in drawbacks as well)
  • APIC: Affine Particle In Cell. Jiang et al., 2015, video
    • less noisy than FLIP and retains low dissipation
    • still uses MAC grid like PIC/APIC but transfer happens from/to not only a velocity but also an locally affine 3x3 velocity matrix

MAC Grid (also staggered Grid): Marker and cell grid (named after how it's used originally). Instead of storing velocities at the grid cells, x/y/z components are offset which makes finite differences at the cells much easier and avoids skipping over mid point. Avoids pressure oscillations ("checkerboard pressure") caused by skipping mid cell in central differences. Makes boundary equations a bit harder. Another way of avoiding oscillations in a collocated(/uniform/cell-centered) grid is Rhie-Chow interpolation.

Excellent explanation for how to solve with MAC grids (but also very helpful for understanding grid based approaches in general): Bridson & Müller-Fischer, Fluid Simulation Siggraph 2007 Course Notes

Houdini docs:

The Splashy Kernel corresponds to the FLIP/PIC method of velocity transfer, while the Swirly Kernel is implemented with the newer APIC method.

Some old docs talk a bit about pro/con of SPH versus FLIP: Bottom line FLIP is more stable such less substepping required and easier to work with albeit bad for slow moving (-> APIC alleviates that).

Suffers typically from volume loss, recent approach aiming to fix this: Kugelstadt et al., Implicit Density Projection for VolumeConserving Liquids

Basic algorithm

  • transfer velocity to grid
  • apply forces to grid
  • solve pressure equation & enforce boundary conditions
  • transfer velocities back to particles
    • for PIC: Replace
    • for FLIP: Update particle velocities with difference on grid since transfer to grid
  • advect particles

Pressure solver

Unlike with SPH, solving equation as Ax=b (with x as pressure and b as divergence) problem is typical. A is a quite simple coefficient matrix that describes relationship between neighboring pressures. Note that the dimension of A for a 3D fluid problem is (width * height * depth)², but it is very sparse, symmetric and positive definite.

Most popular solver for this is Preconditioned Conjugate Gradient Method (PCG), but Jacobi iterations work well enough for getting started on a toy example. Implementing it requires only some vector math, understanding it is a bit harder. Jonathan Richard Shewchuk 1994, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain Edition 1 1/4 explains it step by step and in detail. Conjugate Gradient Method, 2007, a handout for a lecture of theIowa State University, gives a shorter overview.

Chentanez & Müller 2011, A Multigrid Fluid Pressure SolverHandling Separating Solid Boundary Conditions describe a more complex multi-grid solver which also avoids the typical stickiness problems of the "standard" boundary condition.

GPU implementations

Transfer Particles->Grid is surprisingly tricky! Can do either gather or scatter.
Gather particles per grid cell: Need to find relevant particles somehow Scatter particles to grid: Read/write contention - there's atomic float addition is not widely available, and even when it is it's slow.

Some ideas for ideas for conflict resolution in Ming et al 2018, GPU Optimization of Material Point Methods

Doing staggered grids in particular can get very expensive since now particles scatter over an even wider area.

Some pressure solvers are not well suited for GPU impl: For example Preconditioned Conjugate Gradient Method with Incomplete Cholesky preconditioner is very effective and popular (see Bridson, Fluid Simulation for Computer Graphics), but different pre-conditioners are much better suited for a GPU impl see e.g. Austin Eng's WebGL FLIP implementation.

MPM

Generalization to wider range of physical phenomena: Material Point Method (MPM). Excellent overview and redirect to various resources, by Niall Tl.

State of the art (and easier and faster than previous!) is Moving Least Squares MPM (MLS-MPM) by SIGGRAPH 2018, Yahnming Hu et al..

In a way (MLS-)MPM is closer again to SPH as more computation is happening on per-particle basis than rather on the grid as is the case with "pure" APIC fluids. Most notably there is (typically) no pressure solver or similar applied to the grid (rather state equations like in WCSPH). This makes the implementation easier & faster but also requires lower time steps. Note that due to the different gradient requirements the grid in (MLS-)MPM is not staggered!
From the MLS-MPM paper:

Interestingly, our derivation shows that MPM, although seemingly quite different from purely Lagrangian meshless methods, can be treated as a modified element-free Galerkin (EFG) method [Belytschko et al.1994], where the background Eulerian grid merely acts as a helper structure for accelerating MLS interpolation from particle neighbor regions.

Finite Difference/Volume/Element (Eulerian)

Finite Element

Popular in structural analysis but not so much in fluid dynamics since it requires higher cell resolution wiki

Lattice Boltzmann methods

In pure form simple to implement grid propagation scheme. Looks at fluid with Boltzmann equations instead of Navier Stokes -> probability distributions of particles. Simulates steps of particles moving to neighbor gridcells and colliding with statistical functions. Needs very very small time steps and high grid resolution for more macroscopic phenomena. Scales very well with GPUs since can operate on homogenous grid cells individiually.

Rendering

Screen Space from particles

Screen space splatting of fluid particles, described well in these GDC 2010 slides, some improvements & paper here: Screen space fluid rendering with curvature flow And even more improvements: Adding a foam layer that can also be underwater! A Layered Particle-Based Fluid Model for Real-Time Rendering of Water, Bagar et al 2010 adding handling of a foam layer.
Improved screen space filtering A Narrow-Range Filter for SCreen-Space Fluid Rendering, Truong et al. 2018

Screen Space Foam Reandering, Akinci et al. 2013

Real-Time Screen-Space Liquid Rendering with Complex Refractions, Imai et al. 2016 (Sounds at first like they do two layers each with two refactions, but that's only technically true since one of the layers consists entirely of splashes)

todo:

level sets? surface extraction?

Important/big/famous/interesting simulation software packages

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment