I’d like to compute the scalar streamfunction of a given irrotational vector field. As mentioned in the documentation (p. 538) the streamfunction satifies a laplacian equation. I don’t know how to prescribe the right boundary condition in the general case so that the gradient is orthogonal to the vector field.

Moreover instead of using the laplacian equation, is there a way to solve that problem using directly the equation grad(u)=(-Fy,Fx) ?

I will answer to your first question, that is how to impose boundary conditions for a
case where the streamfunction is defined by a laplace/poisson equation. In the following I will deal with the streamfunction of a fluid domain that is laplace(\Psi) = rot(u) = -vorticity. In this sense this problem is more general that what you propose and it is exactly what you want in the case vort = 0.
Consider for simplicity and without loss of generality the two dimensional cartesian case.
In such a case the velocity field is described by:
D_y \Psi = u and D_x \Psi = -v where D_i is the partial derivative operator and \Psi the streamfunction.

In such a case supposing that we know a priori the vector field (u,v) (velocity vector field)
we can impose the boundary conditions by integration over the domain.
For exemple:

where vort is the vorticity field, label is the labels of every boundary except 2 which is the label of a solid wall where you set a constant value for the streamfunction. If the Dirichlet boundary condition is not set the problem has unique solution up to a constant. Then we need to fix that constant, which is easily done in the case where one of the components of the vector field is equally zero, i.e. a solid wall, an axe of symmetry…