lbm.calculate_force_on_boundaries only works on a stationary mesh

[trparry]

Based on my testing, lbm.calculate_force_on_boundaries only works on a stationary mesh. With the default D3Q19 velocity set, I moved a mesh with mesh->rotate or mesh->translate, then I revoxelize and flag with type S and run for 50 LBM timesteps. Then I call lbm.calculate_force_on_boundaries. I do this repeatedly in a while loop as shown. The simulation diverges after 3 cycles of the while loop where the forces are NaN, and the flow field disappears (I think it might be NaN).

I can remove lbm.calculate_force_on_boundaries and the simulation runs without issues.

When I rotate less per loop (ex: radians(0.01f)) I'm able to do 5-6 while loops but then I still diverge. I think its probably the mid-grid bounce back BC where a large movement of the mesh causes issues with lbm.calculate_force_on_boundaries. However, this BC is obviously working without calculating the forces, so why would calculating the forces break it?

Any idea what's going on? It's strange that if I calculate the forces after the while loop its fine, but when this command is inside the while loop problems occur.

I also tried FP16C, and lbm.calculate_force_on_boundaries works fine for each loop, but the fluid domain looks "unphysical" where the entire domain is turbulence.

You can try this on any mesh:

void main_setup() {
    // make a list of all variables you have in SI units (m, kg, s)
    const float si_x = 0.25f; // 0.25m cylinder length
    const float si_u = 5.0f; //wind speed
    const float si_rho = 1.225f; // 1.225kg/m^3 air density
    const float si_nu = 1.48E-5f; // 1.48E-5f m^2/s kinematic shear viscosity of air

    // set velocity in LBM units, should be between 0.01-0.2
    const float lbm_u = 0.1f;
    const float lbm_rho = 1.0f; // density in LBM units always has to be 1
    const float lbm_x = 256u;

    // set unit conversion factors between SI units and LBM units
    units.set_m_kg_s(lbm_x, lbm_u, lbm_rho, si_x, si_u, si_rho);

    // compute kinematic shear viscosity in LBM units
    const float lbm_nu = units.nu(si_nu);

    // set grid resolution based on lbm_x
    const uint Nx = to_uint(1.25*lbm_x);
    const uint Ny = to_uint(2*lbm_x);
    const uint Nz = to_uint(lbm_x/2);

    LBM lbm(Nx, Ny, Nz, lbm_nu); // create LBM object

    // load geometry from stl file, mark all grid points that belong to your geometry with (TYPE_S)
    const float size = 1.0f * (float)lbm_x;
    float3 center = float3(lbm.center().x, lbm.center().y, lbm.center().z);
    const float3x3 rotation = float3x3(float3(0, 0, 1), radians(90.0f));
    Mesh* mesh = read_stl("<path to your stl file>", lbm.size(), center, rotation, size);
    voxelize_mesh_hull(lbm, mesh, TYPE_S);

    // set box boundary conditions
    for (uint n = 0u, x = 0u, y = 0u, z = 0u; n < lbm.get_N(); n++, lbm.coordinates(n, x, y, z)) {
        if (!(lbm.flags[n] & TYPE_S)) lbm.u.y[n] = lbm_u; // initial velocity
        if (x == 0u || x == Nx - 1u || y == 0u || y == Ny - 1u || z == 0u || z == Nz - 1u) lbm.flags[n] = TYPE_E;
        if (y == 0u) lbm.u.y[n] = lbm_u; //velocity inlet
    }

lbm.run(0u);
    while (lbm.get_t() < 50000u) {
        for (uint n = 0u; n < lbm.get_N(); n++) lbm.flags[n] &= ~TYPE_S; // clear flags
        const float3x3 rotation = float3x3(float3(0.0f, 0.0f, 1.0f), radians(0.9f)); // create rotation matrix to rotate mesh 
        mesh->rotate(rotation); // rotate mesh
        //voxelize_mesh_hull(lbm, mesh, TYPE_S); // alternative but issue still occurs with this  
        lbm.voxelize_mesh(mesh, TYPE_S);
        lbm.flags.write_to_device(); // lbm.flags on host is finished, write to device now
        lbm.run(50u); 

                // calculate forces on boundaries on GPU, then copy force field to CPU memory
            lbm.calculate_force_on_boundaries(); //remove this and the simulation runs fine.
                // sum force over all boundary nodes marked with TYPE_S
        lbm.F.read_from_device();

                //transition to si units
        const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S);
        const float si_force_x = units.si_F(lbm_force.x); // force in Newton
        const float si_force_y = units.si_F(lbm_force.y);
        const float si_force_z = units.si_F(lbm_force.z);
                // print force components
        print_info("z force = " + to_string(si_force_z) + " N");
        print_info("y force = " + to_string(si_force_y) + " N");
        print_info("x force = " + to_string(si_force_x) + " N");
           }
}

[ProjectPhysX]

I think I know what's going on here. There is this lbm.F data field that stores the 3D force per volume for each grid point. This is used for 2 things:

1. If lbm.F is set unequal to 0 on a fluid node, it is applied as a local force per volume in every time step via Gui forcing. This allows to set inhomogeneous forces depending on the position, i.e. like the gravitational field of a planet.
2. lbm.calculate_force_on_boundaries() computes and sets lbm.F (on the GPU side) is for all solid boundary node. These don't do the Gui forcing as they are not computed. Put 1. and 2. together and you get the crash: lbm.F is set on a boundary node by lbm.calculate_force_on_boundaries(), then the mesh is moved and re-voxelized and the boundary node turns into a fluid node, and suddenly the (large) boundary force is read as a local force per volume on the fluid and the simulation blows up.

Solution: Reset lbm.F to 0 before re-voxelization, so call lbm.F.reset(); lbm.F.write:to:device();. A faster way would be to reset lbm.F directly on the GPU side without slow PCIe memory transfer. I might add a short kernel for this later on.

Alternative solution (more elegant, but please test if this works I'm not sure): Comment out src/defines.hpp line 73. This let's you enable FORCE_FIELD but keep Guo forcing on the fluid disabled.

[ProjectPhysX]

I just tested the alternative solution (it works), and updated the repository. FORCE_FIELD and VOLUME_FORCE can now be used independently.

[trparry]

@ProjectPhysX yep it works, extremely helpful thank you!