Function collide_stream

examples/pde/boltzmann_cfd.cpp:253–363

* @brief Updates the particle distribution functions for the new simulation * frame * */

 *
 */
void collide_stream(Simulation& sim) {
    auto& ux             = sim.ux;
    auto& uy             = sim.uy;
    auto& rho            = sim.rho;
    auto& sigma          = sim.sigma;
    auto& f              = sim.f;
    auto& feq            = sim.feq;
    auto& set_boundaries = sim.set_boundaries;

    auto& ex   = sim.ex;
    auto& ey   = sim.ey;
    auto& wt   = sim.wt;
    auto& ex_  = sim.ex_;
    auto& ey_  = sim.ey_;
    auto& ex_T = sim.ex_T;
    auto& ey_T = sim.ey_T;
    auto& wt_T = sim.wt_T;

    auto density  = sim.density;
    auto velocity = sim.velocity;
    auto reynolds = sim.reynolds;
    auto xcount   = sim.grid_width;
    auto ycount   = sim.grid_height;

    const float viscosity =
        velocity * std::sqrt(static_cast<float>(xcount * ycount)) / reynolds;
    const float tau  = 0.5f + 3.0f * viscosity;
    const float csky = 0.16f;

    auto edotu = ex_ * ux + ey_ * uy;
    auto udotu = ux * ux + uy * uy;

    // Compute the new distribution function
    feq =
        rho * wt * (edotu * edotu * 4.5f - udotu * 1.5f + edotu * 3.0f + 1.0f);

    auto taut =
        af::sqrt(sigma * (csky * csky * 18.0f * 0.25f) + (tau * tau * 0.25f)) -
        (tau * 0.5f);

    // Compute the shifted distribution functions
    auto fplus = f - (f - feq) / (taut + tau);

    // Compute new particle distribution according to the corresponding D2N9
    // weights
    for (int i = 0; i < 9; ++i) {
        int xshift = static_cast<int>(ex_vals[i]);
        int yshift = static_cast<int>(ey_vals[i]);

        fplus(af::span, af::span, i) =
            af::shift(fplus(af::span, af::span, i), xshift, yshift);
    }

    // Keep the boundary conditions at the borders the same
    af::replace(fplus, af::tile(!set_boundaries, af::dim4(1, 1, 9)), f);

    // Update the particle distribution
    f = fplus;