def METHOD_NAME(config, method, pdfs, shear_dir_normal, stencil):
    """Calculates the difference between quilibrium pdf distributions
    with (rho, u) and (rho, u+v) and applies them to out-flowing
    populations in the boundary layer. Returns an AssignmentCollection
    with one Assignment per stencil direction.
    """
    dim = len(stencil[0])
    default_dtype = config.data_type.default_factory()
    # Placeholders indicating a population flows up or down.
    # Will be replaced later using the component of the stencil direction
    # along the shear_dir_normal.
    points_up = sp.Symbol('points_up')
    points_down = sp.Symbol('points_down')
    # Symbol for the coordinate index within the field,
    # used to identify boundary layers
    counters = [LoopOverCoordinate.get_loop_counter_symbol(
        i) for i in range(dim)]
    grid_size = TypedSymbol("grid_size", dtype=default_dtype)
    # +,-1 for upper/lower boundary layers, 0 otherwise.
    # Based on symbolic counters defined above. Only becomes
    # non-zero if the corresponding points_up/down flags
    # are engaged (which is only done for out-flowing populations)
    layer_prefactor = sp.Piecewise(
        (-1,
         sp.And(type_all_numbers(counters[1] <= 0, default_dtype),
                points_down)),
        (+1,
         sp.And(type_all_numbers(counters[1] >= grid_size - 1, default_dtype),
                points_up)),
        (0, True)
    )
    # Start with an equilibrium distribution for a given density and velocity
    delta_pdf_eqs = macroscopic_values_setter(
        method, sp.Symbol("dens"), [
            sp.Symbol("v_0"), sp.Symbol("v_1"), sp.Symbol("v_2")], pdfs)
    # Replace the assignments of (rho,u) by (rho, u+v) - (rho,u)
    ma = []
    for a, c in zip(delta_pdf_eqs.main_assignments, method.stencil):
        # Determine direction of the stencil component in the
        # shear_dir_normal
        if c[shear_dir_normal] == 1:
            up = True
            down = False
        elif c[shear_dir_normal] == -1:
            up = False
            down = True
        else:
            up = False
            down = False
        # Replace (rho,u) by (rho,u+v) in boundary layers
        rhs = sp.simplify(
            a.rhs -
            a.rhs.replace(
                sp.Symbol("u_0"),
                sp.Symbol("u_0") +
                layer_prefactor *
                sp.Symbol("v_s")))
        # Only engage if the population is outflowing. See layer_prefactor
        rhs = rhs.replace(points_up, up)
        rhs = rhs.replace(points_down, down)
        new_a = Assignment(a.lhs, rhs)
        ma.append(new_a)
        print(c, ma[-1])
    # Plug in modified assignments
    delta_pdf_eqs.main_assignments = ma
    return delta_pdf_eqs.main_assignments