Enforcing volumetric velocity condition throughout the domain for the entire simulation isn't physically accurate

Hello, I just wanted to point out that some of the setups in setup.cpp enforce a condition that might not be physically accurate:

In your setups, this is often used:

const ulong N=lbm.get_N(); const uint Nx=lbm.get_Nx(), Ny=lbm.get_Ny(), Nz=lbm.get_Nz(); for(ulong n=0ull; n<N; n++) { uint x=0u, y=0u, z=0u; lbm.coordinates(n, x, y, z);
        // ########################################################################### define geometry #############################################################################################
        if(lbm.flags[n]!=TYPE_S) lbm.u.y[n] = u;
        if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_E; // all non periodic
    }   // #########################################################################################################################################################################################

When you say

if(lbm.flags[n]!=TYPE_S) lbm.u.y[n] = u;

I think this means (please correct me) that during the entire simulation (this is not an IC), the y component of the velocity in every cell that isn't flagged TYPE_S will equal u. This is enforced on the entire volumetric domain. This is not physically accurate because near a solid wall a boundary layer will develop where the y component will gradually approach 0 (no slip wall). A better way:

Enforce this only on the inlet, where all of the boundaries are also TYPE_E:

if (y == 0u) lbm.u.y[n] = u; //inlet
if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==0u||z==Nz-1u) lbm.flags[n] = TYPE_E; // all non periodic

This way, the entire domain is initialized at zero velocity and the y component is free to be whatever the solution wills.

The first way that you've used is probably faster because you initialize the entire domain at a nonzero value. What I've proposed here requires time for the domain to reach a steady value.

To alleviate this issue, it would be nice to be able to initialize the domain at a nonzero velocity, and then allow it to change value. This would mean the first simulation step would use if(lbm.flags[n]!=TYPE_S) lbm.u.y[n] = u;, and then this constraint would be released on any subsequent simulation steps. But, I don't know how to do this.

This is only an initial condition. The velocity y-component is free to change in the entire box, except in the TYPE_E cells at the box boundary, where the initial velocity (0, u, 0) is enforced for all time steps.

If you enforce u only at the inlet side, it is not set at the other TYPE_E walls, and they enforce 0 velocity (but are non-reflecting, so allow tangential flow/outflow). This isn't wrong, but I haven't tested that myself yet.

Note that the initial state of the velocity field remains in lbm.u GPU memory unless lbm.update_fields() is called. Until then, the new velocity field state is only stored in the LBM density distribution functions. To update the velocity field in CPU memory, you also have to call lbm.u.read_from_device().

ProjectPhysX / FluidX3D

Enforcing volumetric velocity condition throughout the domain for the entire simulation isn't physically accurate #60