Hydrostatic initial condition does not provide steady flow

[rodionstepanov]

Hi! The flow in thermal convection case with thermal exp. koef. = 0 should remain zero. However test simulation revels some perturbations. Is it a mistake in setup or normal for LBM?

[ProjectPhysX]

Hi @rodionstepanov,

it's a bit confusing with the many LBM constructors, but here the very last parameter is beta, which is set to 1 and gets correctly passed along internally:

//      Nx , Ny  , Nz , Dx, Dy, Dz, nu   , fx  , fy  , fz     ,sigma,alpha, beta
LBM lbm(32u, 196u, 60u, 1u, 1u, 1u, 0.02f, 0.0f, 0.0f, -0.001f, 0.0f, 1.0f, 1.0f);

Did I understand your question correctly?

Kind regards, Moritz

[rodionstepanov]

Hi @ProjectPhysX ! Sorry for misleading. My point is that case beta=0 means zero buoyancy force. So initial state provided by V=0 and rho=rho_hydrostatic must remain unchanged. But it does not as you can find from the setup below. Why there is no equilibrium?

void main_setup() { // thermal convection; required extensions in defines.hpp: FP16S, VOLUME_FORCE, TEMPERATURE, INTERACTIVE_GRAPHICS
    // ################################################################## define simulation box size, viscosity and volume force ###################################################################
    LBM lbm(32u, 196u, 60u, 1u, 1u, 1u, 0.01f, 0.0f, 0.0f, -0.001f, 0.0f, 0.01f, 0.0f);
    // ###################################################################################### define geometry ######################################################################################
    const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); parallel_for(lbm.get_N(), [&](ulong n) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
        if(y==1) {
            lbm.T[n] = 1.8f;
            lbm.flags[n] = TYPE_T;
        } else if(y==Ny-2) {
            lbm.T[n] = 0.3f;
            lbm.flags[n] = TYPE_T;
        }
        lbm.rho[n] = units.rho_hydrostatic(0.001f, (float)z, (float)Nz-2.0f); // initialize density with hydrostatic pressure
        if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_S; // all non periodic
    }); // ####################################################################### run simulation, export images and data ##########################################################################
    lbm.graphics.visualization_modes = VIS_FLAG_LATTICE|VIS_STREAMLINES;
    lbm.run(0u);
    lbm.u.read_from_device(); println(lbm.u.z[lbm.index(Nx / 2u, Ny / 2u, Nz / 2u)]);
    lbm.run(10u); 
    lbm.u.read_from_device(); println(lbm.u.z[lbm.index(Nx/2u, Ny/2u, Nz/2u)]); wait(); // test for binary identity
}

[rodionstepanov]

@ProjectPhysX , I have a look at free decaying 3D flow in triple periodic cube. It is even simpler case to test conservation laws. I expect to get constant total mass, total momentum and monotonic decrease in total kinetic energy (kinetic helicity too but it is more complex). I didn't find it. I believe code is conservative but what is wrong with my setup below.

void main_setup() { // test free decay; 
    // ################################################################## define simulation box size, viscosity and volume force ###################################################################
    LBM lbm(128u, 128u, 128u, 1u, 1u, 1u, 0.01f);
    // ###################################################################################### define geometry ######################################################################################
    const uint threads = (uint)thread::hardware_concurrency();
    vector<uint> seed(threads);
    for (uint t = 0u; t < threads; t++) seed[t] = 42u + t;
    const uint Nx = lbm.get_Nx(), Ny = lbm.get_Ny(), Nz = lbm.get_Nz(); parallel_for(lbm.get_N(), threads, [&](ulong n, uint t) { uint x = 0u, y = 0u, z = 0u; lbm.coordinates(n, x, y, z);
    lbm.u.x[n] = random_symmetric(seed[t], 0.1f); // initialize velocity with random noise
    lbm.u.y[n] = random_symmetric(seed[t], 0.1f);
    lbm.u.z[n] = random_symmetric(seed[t], 0.1f);
    lbm.rho[n] = 1.0f; 
        }); 
    lbm.graphics.visualization_modes = VIS_FLAG_LATTICE | VIS_FIELD;
    lbm.run(0u); // initialize simulation
    while (true) { // main simulation loop
        lbm.u.read_from_device();
        lbm.rho.read_from_device();
        double en = 0, mx = 0, my = 0, mz = 0,mass=0;
        parallel_for(lbm.get_N(), [&](ulong n) {
            en += lbm.rho[n] * (lbm.u.x[n] * lbm.u.x[n] + lbm.u.y[n] * lbm.u.y[n] + lbm.u.z[n] * lbm.u.z[n]);
            mx += lbm.rho[n] * lbm.u.x[n];
            my += lbm.rho[n] * lbm.u.y[n];
            mz += lbm.rho[n] * lbm.u.z[n];
            mass += lbm.rho[n];
            });
        string s = to_string(lbm.get_t()) + " " + to_string(mass) + " " + to_string(mx) + " " + to_string(my) + " " + to_string(mz) + " " + to_string(en) + "\n";
        print_info(s);
        lbm.run(1u);
    }

}

[ProjectPhysX]

Hi @rodionstepanov,

in the first convection case, there is two opposing forces:

• the downward volume force fz
• the pressure/density gradient from units.rho_hydrostatic(...)

In theory both perfectly cancel out. In practice, they only cancel if force and density gradient both go to zero, and for finite values there still is a bit of wobble in the beginning of the simulation.

The second case, there is a mistake in the summation. Be cautious with multithreaded parallel_for when doing horizontal sums: you need one accumulator for every thread. Otherwise, concurrent threads constantly overwrite each other's accumulator sums.

Use either a single-threaded summation

double accumulator = 0.0;
for(ulong n=0ull; n<lbm.get_N(); n++) {
    accumulator  += (double)lbm.rho[n];
}

or multiple accumulators:

const uint threads = thread::hardware_concurrency();
vector<double> accumulators(threads, 0.0);
parallel_for(lbm.get_N(), threads, [&](ulong n, uint t) {
    accumulators[t] += (double)lbm.rho[n];
});
double accumulator = 0.0;
for(uint t=0u; t<threads; t++) {
    accumulator += accumulators[t];
}

Expectation for the setup is that mass and momentum (components) remain constant (apart from floating-point round-off), and kinetic energy decays approximately exponential. This is the results with fixed summation:

Kind regards, Moritz

[rodionstepanov]

In theory both perfectly cancel out. In practice, they only cancel if force and density gradient both go to zero, and for finite values there still is a bit of wobble in the beginning of the simulation.

I see your point but how do you know it is "wobble" or mistake in the code? Do you know any publication on LMB where such wobbling demonstrated analytically?

The second case, there is a mistake in the summation.

Thanks a lot. I follow your correct version.

Expectation for the setup is that mass and momentum (components) remain constant (apart from floating-point round-off), and kinetic energy decays approximately exponential.

I see floating-point round-off effect in mass and momentum but for energy an error is much larger. Look at t=2 the energy (last column) increases significantly.

Info: 0 2.0971520000000004E6 -1.7117407679557801E1 -2.8646887183189396E1    1.2938460651785136E1 2.0971056823777810E4              
Info: 1 2.0971520000183581E6 -1.7117073257551247E1 -2.8646542910754231E1    1.2938803720443979E1 3.5405405270094938E3                            
Info: 2 2.0971519998005035E6 -1.7117072651320768E1 -2.8646542668130546E1    1.2938808677505649E1 4.4315150465497251E3                       

[ProjectPhysX]

The density wobbling with volume force is neither an analytical phenomenon nor a mistake in the code. It's a result of LBM working only in the weakly compressible regime. A few % variation in density are already quite a lot for LBM, and the wobbling is expected. I don't know if this is documented in literature somewhere, your best bet is the Krüger LBM book.

What you see in the energy is not floating-point round-off, it's a result of the LBM collision operator. For very large velocity gradient (like in this setup when neighboring cells have almost opposite velocity), LBM does some over-relaxation, leading to oscillations between even and odd time steps. This is well-known behavior.

[rodionstepanov]

@ProjectPhysX I see your point. Thank you for explanation. I'm just wondering how I can distinguish acceptable violations of conservation laws from things going wrong.

[ProjectPhysX]

There is no simple answer to this question. Working with a particular numerical model for a long time helps, and also comparing results from different implementations of the same thing. LBM has its weirdness sometimes, even with fully correct implementation. Outside a certain parameter range it may completely break down. Most literature doesn't enlighten you that 99% of their simulation attempts went horribly wrong and what the paper presents is the cherry-picked 1% success.

Of course I made sure that every part of my implementation is correct, with at least triple redundant tests for everything, but don't trust my code alone. If you see something that seems wrong, and where other LBM codes (such als OpenLB) produce different results, please let me know.

Also don't blindly trust numerical models, they are only models after all, with limitations and sometimes characteristic artifacts. Experimental observation of nature is always the judge.

Anyways, I'm closing this now, thank you for the great discussion. Let me know when you encounter further problems and/or questions. Cheers!