ProjectPhysX / FluidX3D

The fastest and most memory efficient lattice Boltzmann CFD software, running on all GPUs via OpenCL. Free for non-commercial use.
https://youtube.com/@ProjectPhysX
Other
3.48k stars 281 forks source link

Vibrating Lines #200

Closed Meerkov closed 1 week ago

Meerkov commented 1 week ago

image

Using one of the default setups, I noticed that the initial state of random noise results in visible grid patterning.

Using recommended setup for Rayleigh Benard Convection demo: TEMPERATURE VOLUME_FORCE D3Q19 FP16S

Is this expected behavior? It looks wrong

The default coloration might be hard to see, so here it is when executed in my custom environment. The vibrating and self-orienting pattern of the random noise should become apparent. image

Edit: In my sample image of my setup I'm using velocity 0.15 instead of 0.015 as in the demo setup. However it still reproduces, just with lower visibility. Here's example data from a single pixel: Format is "Frame Number, FPS" then the Velocity/c

478 FPS: 253
0.009542987
479 FPS: 251
0.008611688
480 FPS: 250
0.00957845
481 FPS: 250
0.008645571
482 FPS: 251
0.009601751
483 FPS: 251
0.0086692115
484 FPS: 252

I picked a pixel at random, and it's changing by about 10% up and down every frame, in a vibrating fashion.

Meerkov commented 1 week ago

I removed all the extensions, to try to reduce the problem. You may need to raise the contrast of your monitor to clearly see the patterns, but you can also export the data as I did above to see the raw values are indeed fluctuating, somewhere between 10% and 100% every other pixel. image Here, I ran 200 frames on six separate instances with torus wrapping (appended into a 6x1 image). The second image had a large droplet of density added.

I have done something special here, j[ 5] = n/*xp+yp+z0*/; . This obviously dampens the fluid, but quickly the vibrations are removed.

Please note that there is a macro scale diagonal directional preference, but zero micro scale vibrating lines.

image Here is an image using the original D2Q9 setup with the j[5] neighbor as normal. Notice here there is both large scale waves, and microscale velocity waves. The density is actually not impacted (!!). I think this might mean that there is some issue with how the velocity fields are exported, but not how the calculations are performed. otherwise I would expect to see these ripples in the density field. But no such artifacting exists. image

This is not limited to D2Q9. Here it is again, on a D3Q27 (with 1 z layer). image And with 10 z layers image Going back to where we started, at D2Q9, and negating the j[7] direction by setting that neighbor to n, we see the direction of our macroscale ripple change, but again we return to removing all microscale rippling artifacts. I lowered the large density blob so it's more clear that all artifacting is removed. image

However, it is not sufficient to remove any direction. Removing the j[3] direction results in the microscale artifacts remaining. image

You can try using an odd dimensionality (81 pixels instead of 80) but this won't work. You can try changing the nu variable to 0.1 which reduces the effect but does not negate it.

However, if you negate any not straight neighbor, the vibrational lines disappear nearly instantaneously.

Given that they show up in both my UI, and FluidX3D's UI, it seems like it's actually a problem with the exporting of the uxn, uyn, and uzn variables

ProjectPhysX commented 1 week ago

Hi @Meerkov,

the Rayleigh-Benard setup intentionally initializes the velocity with a bit of random noise:

        lbm.u.x[n] = random_symmetric(seed[t], 0.015f); // initialize velocity with random noise
        lbm.u.y[n] = random_symmetric(seed[t], 0.015f);
        lbm.u.z[n] = random_symmetric(seed[t], 0.015f);

This is to trigger the Rayleigh-Benard instability. It is not a bug.

Kind regards, Moritz

Meerkov commented 1 week ago

You didn't read what I wrote at all, did you.

Meerkov commented 1 week ago

Plotting X velocity while cross section looking at the X Y plane of a 64x64x64 grid: image Plotting the X velocity while looking through the Y Z plane image

Plotting X velocity through the X Z plane image

It should be clear that the X velocity is aligned to the X direction in a pattern. It's not random.

ProjectPhysX commented 1 week ago

Oh sorry, I didn't read correctly.

I can reproduce it, there is some horizontal/vertical stripe patterns visible in the velocity slice when setting GRAPHICS_U_MAX very low. grafik

This is expected behavior with LBM. A cartesian grid by its nature is not anisotropic and there is always some leftover grid artifacts. The SRT collision operator also does some over-relaxation that can result in cells oscillating in state, both in even/odd time step cadence and in spatial neighbors. That causes the stripe patterns.

At such low velocities this is usually not a problem. The same phenomenon becomes macroscopically visible and problematic when setting lbm_u>>0.1.

Kind regards, Moritz

Meerkov commented 1 week ago

Thanks for re-reading.

This is expected behavior with LBM. A cartesian grid by its nature is not anisotropic and there is always some leftover grid artifacts.

You're closer, but I think there is still an issue.

This applies at all scales of random noise. lbm = Lattice(D=3, Q=27, res=64, Nz=64, sides=6, nu=0.02) I'll show you. Each image will be 15 frames after start, showing only the X velocity through the X Y frame. This isn't just any pattern, but an alternating pattern every other frame, which is directly related to the Esoteric Pull algorithm in some way. lbm.velocity[:]=cp.random.uniform(-0.15,0.15, (lbm.N*3), dtype=cp.float32) image lbm.velocity[:]=cp.random.uniform(-0.55,0.55, (lbm.N*3), dtype=cp.float32) image lbm.velocity[:]=cp.random.uniform(-0.05,0.05, (lbm.N*3), dtype=cp.float32) image

Meerkov commented 1 week ago

image Density calculation, 15 frames, lbm.velocity[:]=cp.random.uniform(-0.05,0.05, (lbm.N*3), dtype=cp.float32) Same as the last image in the previous, but notice there is no grid pattern. The density field isn't impacted, only the velocity field.

ProjectPhysX commented 1 week ago

The alternating pattern between even/odd time steps comes from the SRT collision operator when kinematic viscosity lbm.nu is small, that is a known phenomenon. Try increasing kinematic viscosity lbm.nu in initialization. Patterns should be vastly reduced.

The streaming scheme does not affect the arithmetic at all. You can replace the streaming algorithm with any other and it is binary identical. It just changes the memory locations where the DDFs are stored. I provide some other streaming algorithms here in the appendix.

Meerkov commented 1 week ago

100 frames with nu=0.02 image with nu=0.22 image

The artifacts are reduced, sorta, but the fluid also behaves differently. TRT didn't seem to change much. TRT with nu=0.02 image

I keep looking at that density field, and how it doesn't have these artifacts. I'm not sure how it could be that if the velocity truly has obvious alternating patterns, that those patterns wouldn't create the same ripples in the density distribution image

Meerkov commented 1 week ago

I took for a long time reading the code, and I can't determine what causes it.

I decided to use a simple work around, and average the uxn, uyn, uzn variables over 2 frames for output purposes only (not impacting the internal calculations). This removes the artifact and provides the expected smooth gradient. image

ProjectPhysX commented 1 week ago

Averaging 2 consecutive time steps is indeed a way to smooth over the even/odd oscillations of LBM. Adjusting viscosity lbm_nu alone of course changes the physical behavior. But unit conversion usually leaves one free parameter, so you can adjust lbm_nu and velocity lbm_u together without changing the physics; the Reynolds number $Re=\frac{u}{\nu}\cdot L$ is kept unchanged then. This comes with limits though, as too small lbm_u will blow up runtime and too large lbm_u will be unstable.

I still owe you an explanation of where exactly the time oscilations come from. The SRT collision operator (implemented here) reads:

$$f_i(\vec{x}, t+\Delta t) = (1-\frac{\Delta t}{\tau})\ f_i^{temp}(\vec{x}, t)+\frac{\Delta t}{\tau}\ f_i^{eq}(\vec{x}, t)$$

For very small viscosuty $\nu\rightarrow 0$, the the relaxation time $\tau=\frac{\nu}{c^2}+\frac{\Delta t}{2}\rightarrow\frac{1}{2}$. Inserting that gives:

$$f_i(\vec{x}, t+\Delta t) = - f_i^{temp}(\vec{x}, t)+2\ f_i^{eq}(\vec{x}, t)$$

If we assume small equilibrium DDFs $f_i^{eq}(\vec{x}, t)\approx 0$ here for simplicity, then the DDFs

$$f_i(\vec{x}, t+\Delta t) = - f_i^{temp}(\vec{x}, t)$$

oscillate between +/- sign in even/odd time steps.

Kind regards, Moritz

Meerkov commented 1 week ago

I appreciate the effort to explain it. Switching to TRT did not change it (this is not SRT specific).

I tried changing just nu and u but this didn't work.

Can you provide more specific guidance on how to manipulate the default settings to remove the effect? Does the density also need to increase?

Regardless, I think I will stay with my simple hack for now, but I will continue to experiment with these settings based on your feedback.

Thanks for taking a look

ProjectPhysX commented 1 week ago

The oscillation behavior is universal also with other collision operators, SRT is just the easiest to explain.

Parametrization tuning is an art in itself - you can only experiment here: vary u between 0.001f and 0.1f, and adjust nu accordingly to keep Re constant, and see what works best. There is no one-choice-fits-all, and sometimes there is no decent choice at all.

The density should not be changed, it always has to be close to 1.0f in LBM, otherwise LBM starts to behave strange and unstable. The choice of rho_avg=1 goes back to floating-point arithmetic being most accurate arounf 1.0f.

Good luck!