@ti.kernel
def grid_op():
    for (i, j, k) in grid_m_in:
        inv_m = (1 / (grid_m_in[(i, j, k)] + 1e-10))
        v_out = (inv_m * grid_v_in[(i, j, k)])
        v_out[1] -= (dt * gravity)
        if ((i < bound) and (v_out[0] < 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((i > (n_grid - bound)) and (v_out[0] > 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((k < bound) and (v_out[2] < 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((k > (n_grid - bound)) and (v_out[2] > 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((j < bound) and (v_out[1] < 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
            normal = ti.Vector([0.0, 1.0, 0.0])
            lsq = ti.sqr(normal).sum()
            if (lsq > 0.5):
                if ti.static((coeff < 0)):
                    v_out[0] = 0
                    v_out[1] = 0
                    v_out[2] = 0
                else:
                    lin = (ti.transposed(v_out) @ normal)(0)
                    if (lin < 0):
                        vit = (v_out - (lin * normal))
                        lit = (vit.norm() + 1e-10)
                        if ((lit + (coeff * lin)) <= 0):
                            v_out[0] = 0
                            v_out[1] = 0
                            v_out[2] = 0
                        else:
                            v_out = ((1 + ((coeff * lin) / lit)) * vit)
        if ((j > (n_grid - bound)) and (v_out[1] > 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        grid_v_out[(i, j, k)] = v_out
@ti.kernel
def grid_op():
    for (i, j, k) in grid_m_in:
        inv_m = (1 / (grid_m_in[(i, j, k)] + 1e-10))
        v_out = (inv_m * grid_v_in[(i, j, k)])
        v_out[1] -= (dt * gravity)
        if ((i < bound) and (v_out[0] < 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((i > (n_grid - bound)) and (v_out[0] > 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((k < bound) and (v_out[2] < 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((k > (n_grid - bound)) and (v_out[2] > 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        if ((j < bound) and (v_out[1] < 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
            normal = ti.Vector([0.0, 1.0, 0.0])
            lsq = ti.sqr(normal).sum()
            if (lsq > 0.5):
                if ti.static((coeff < 0)):
                    v_out[0] = 0
                    v_out[1] = 0
                    v_out[2] = 0
                else:
                    lin = (ti.transposed(v_out) @ normal)(0)
                    if (lin < 0):
                        vit = (v_out - (lin * normal))
                        lit = (vit.norm() + 1e-10)
                        if ((lit + (coeff * lin)) <= 0):
                            v_out[0] = 0
                            v_out[1] = 0
                            v_out[2] = 0
                        else:
                            v_out = ((1 + ((coeff * lin) / lit)) * vit)
        if ((j > (n_grid - bound)) and (v_out[1] > 0)):
            v_out[0] = 0
            v_out[1] = 0
            v_out[2] = 0
        grid_v_out[(i, j, k)] = v_out