One of the nice things about #SPH (as in #SmoothedParticleHydrodynamics —the acronym has some other meanings) is that its standard formulation is “weakly compressible” (pressure is derived from an equation of state, but for the most part the compressibility can be ignored), which makes the method itself “embarrassingly parallel” with trivial implementations: there's a natural map from particle to work-item that fits GPUs really well (or does it? https://www.youtube.com/watch?v=zxV846HIKP)
1/n
The problem with WCSPH and the explicit integration it's normally used with is that explicit integration has stringent requirements for the time step. When simulating inviscid (or low-viscosity) fluids, like water, the more stringent condition typically comes from the sound speed: water, for example, has a sound speed of ~1.5km/s, so at an inter-particle spacing of ~10cm, you'd have timesteps in the order of 10^-5.
2/n
In WCSPH this is generally worked around by “cheating”: instead of using the actual speed of sound of the fluid, one uses a fictitious speed of sound which is “just enough” to keep the fluid weakly compressible, which depends on what the largest speed you expect from your problem is, including the “free fall” speed: pick at least 10 time that as your speed of sound and your results will be good enough.
3/n
Remember now yesterday's thread about dimensionality? The speed of sound limitation on the time-step is *linear*: double the resolution, you get half the time step.
If you're doing VISCOUS flow, though, you get a time step condition that is QUADRATIC: double the resolution, you get a time step that is FOUR TIMES smaller! It's trivial to get to time steps in the order of 10^-8 like this.
Guess what we need for our lava simulations?
4/n
Hint: water has a kinematic viscosity of around 10^-6 m^2/s, lava *at its fluidest* is in the order of 10^-2 m^2/s, and can get to hundreds of time that. So yeah, we need that. But while you can get around the sound speed limitation with an artificially lower sound speed if you don't actually care about modeling acoustics, you cannot model diffusive problem (such as those with viscous fluids) with a different viscosity: so there's no cheating there!
5/n
This, BTW, is why we developed a semi-implicit formulation that handles the viscous part (only) implicitly:
https://www.sciencedirect.com/science/article/abs/pii/S002199911830593X
https://www.sciencedirect.com/science/article/abs/pii/S0021999122004752
(the latter may be of interest to you if you use #BiCGSTAB or #ConjugateGradient methods to solve linear systems that have nothing to do with SPH, BTW).
And now comes the new difficulty: THERMAL!
6/n
Thermal conduction and radiation are also diffusive processes. So the time step constraints also go with the resolution SQUARED.
In our case we don't worry too much about conduction actually, because lava has “nice” thermal diffusivity coefficients: the viscous time-step dominates.
Radiation OTOH? Oh boy. That's a process whose diffusive coefficient goes with the (absolute) temperature to the FOURTH POWER.
7/n
Radiative losses behave following the Stefan-Boltzmann law
https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_law
and while the leading coefficient is ~5.6*10^-8 W/m^2/K^4, at 400K we're still talking about ~1400 W/m^2, i.e. 1400 W/m^2/s. And that's not really a problem when working with bodies whose size can be measured in meters, but if you need fine resolutions (say, around 1mm) because the thickness of the body is measured in centimeters … well, there you start having problems.
8/n
So … it looks like we'll have to start looking implicit treatment for radiative losses too. I'm not sure I'm looking forward to that.
9/9
Giuseppe Bilotta