Issue with solving coupled partial differential equations with NDSolve for certain boundary conditions

Question

I am trying to solve the coupled PDEs for functions u, v and w. The main problem I am facing is that I am unable to understand how to input boundary condition of this type in NDSolve: D[u[x, y], y] + D[v[x, y], x] == (-u[x, y]/constant) /. y -> ymax

I read about NeumannValue in the Mathematica Documentation, but still could not realize about implementing it in my context.

This is my code:

(*Equations: *)
eqn1 := D[u[x, y], x] + D[v[x, y], y] == 0
eqn2 := \[Alpha]*D[w[x, y], x] + \[Beta] Laplacian[u[x, y], {x, y}] ==
   u[x, y]
eqn3 := \[Alpha]*D[w[x, y], y] + \[Beta] Laplacian[v[x, y], {x, y}] ==
   v[x, y] 


(* CONSTANTS: *)
xmin := -50
xmax := 50
ymin := -5
ymax := 5

(* \[CapitalOmega] is the area we are working on*)
\[CapitalOmega] = Rectangle[{xmin, ymin}, {xmax, ymax}];


(*Boundary Conditions: *)
bcnd1 := v[x, ymax] == 0 
bcnd2 := v[x, 
   ymin] == \[Gamma] DiracDelta[x + x0] - \[Gamma] DiracDelta[
     x - x0] 
bcnd3 := D[u[x, y], y] + D[v[x, y], x] == -u[x, y]/\[Delta] /. 
  y -> ymax 
bcnd4 := D[u[x, y], y] + D[v[x, y], x] == u[x, y]/\[Delta] /. 
  y -> ymin 

(*other boundary condition to determine w*)
bcnd5 := w[xmax, y] == w0  
bcnd6 := w[-xmax, y] == w0




solN = NDSolve[{eqn1, eqn2, eqn3, bcnd1, bcnd2, bcnd3, bcnd4, bcnd5, 
     bcnd6} /. {\[Alpha] -> 1, \[Beta] -> 1, \[Gamma] -> 1, \[Delta] ->
       1, x0 -> 4, w0 -> 0.1}, {w, u, 
    v}, {x, y} \[Element] \[CapitalOmega]];


(*Output: *)
(*
NDSolve::fembdnl: The dependent variable in (u^(0,1))[x,5]+(v^(1,0))[x,5]==-u in the boundary condition DirichletCondition[(u^(0,1))[x,5]+(v^(1,0))[x,5]==-u,y==5.] needs to be linear.

NDSolve::fembdnl: The dependent variable in (u^(0,1))[x,5]+(v^(1,0))[x,5]==-u in the boundary condition DirichletCondition[(u^(0,1))[x,5]+(v^(1,0))[x,5]==-u,y==5.] needs to be linear.
*)

In the end, I have to obtain the vector {u[x,y],v[x,y]} and the function w[x,y] and plot it in the region [CapitalOmega] I have spent much time to try and do this, but could not. I will be very grateful, if somebody helps me to do so.

Answer 1

To reproduce Figure 3 from the paper Non-local transport and the hydrodynamic shear viscosity in graphene we use FDM code from here.

 (*CONSTANTS:*)rule = {\[Alpha] -> 1, \[Beta] -> 1, \[Gamma] -> 
   1, \[Delta] -> 1, x0 -> 1, w0 -> .0, sigma -> 2, 
  p0 -> 0.}; xmin = -5;
xmax = 5;
ymin = -.5;
ymax = .5;

XYgrid[dom_List, pts_List] := 
  MapThread[
   N@Range[Sequence @@ #1, Abs[Subtract @@ #1]/#2] &, {dom, 
    pts - 1}];
BoundaryIndex[xgridlen_, ygridlen_] := 
  Module[{tmp, left, right, bot, top}, 
   tmp = Table[(n - 1) ygridlen + Range[1, ygridlen], {n, 1, 
      xgridlen}]; {left, right} = tmp[[{1, -1}]]; {bot, top} = 
    Transpose[{First[#], Last[#]} & /@ tmp]; {top, right[[2 ;; -2]], 
    bot, left[[2 ;; -2]]}];
FDMat[deriv_, xygrid_, difforder_] := 
 Map[NDSolve`FiniteDifferenceDerivative[#, xygrid, 
     "DifferenceOrder" -> difforder]["DifferentiationMatrix"] &, deriv]
{domain, pts, difforder} = {{{-5, 5}, {-.5, .5}}, {100, 10}, 4};
xygrid = XYgrid[domain, pts]; {nx, ny} = 
 Map[Length, xygrid]; {top, right, bot, left} = 
 BoundaryIndex[nx, ny]; {dx, dy, dx2, dy2, dxy} = 
 FDMat[{{1, 0}, {0, 1}, {2, 0}, {0, 2}, {1, 1}}, xygrid, 
  difforder]; boundaries = Join[top, right, bot, left]; sgrid = 
 Flatten[Outer[List, Sequence @@ xygrid], 1]; 
Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}];
{u, v, w, p} = MakeVariables[{uu, vv, ww, pp}, nx ny];

sx1 = Table[
  If[Abs[sgrid[[bot]][[i, 1]] + x0] <= 1/20, 1, 0] /. rule, {i, 
   Length[bot]}]; sx2 = 
 Table[If[Abs[sgrid[[bot]][[i, 1]] - x0] <= 1/20, 1, 0] /. rule, {i, 
   Length[bot]}];

eqn1 = (dx2 + dy2) . p - dx . u - dy . v; eqn4 = (dx2 + dy2) . w;
eqn2 = \[Alpha]*dx . w + \[Beta] (dx2 + dy2) . u - u;
eqn3 = \[Alpha]*dy . w + \[Beta] (dx2 + dy2) . v - v;
(*bcnd1:*)eqn3[[top]] = v[[top]];
(*bcnd2:=*)eqn3[[bot]] = v[[bot]] - \[Gamma] sx1 + \[Gamma] sx2;
(*bcnd3:*)eqn2[[top]] = (dy . u + dx . v)[[top]] + u[[top]]/\[Delta];
(*bcnd4:*)eqn2[[bot]] = (dy . u + dx . v)[[bot]] - u[[bot]]/\[Delta];
(*bcnd5:=*)eqn1[[left]] = p[[left]] - p0;
(*bcnd6:=*)eqn1[[right]] = p[[right]] - p0; 
eqn1[[top]] = (dy . p)[[top]]; eqn1[[bot]] = (dy . p)[[bot]];
eqn2[[left]] = (dx . u)[[left]]; eqn2[[right]] = (dx . u)[[right]]; 
eqn3[[left]] = (dx . v)[[left]]; eqn3[[right]] = (dx . v)[[right]]; 
eqn4[[left]] = (dx . w)[[left]] - w0;
(*bcnd6:=*)eqn4[[right]] = (dx . w)[[right]] - w0; 
eqn4[[top]] = (dy . w)[[top]]; 
eqn4[[bot]] = (dy . w)[[bot]] - 
  sigma (\[Gamma]  sx1 - \[Gamma]  sx2); eqs = 
 Join[eqn1, eqn2, eqn3, eqn4]; var = Join[u, v, w, p]; {vec, mat} = 
 CoefficientArrays[eqs /. rule, var];
sol = LinearSolve[mat, -vec]; ns = Length[sol]; phi = 
 Take[sol, -ns/4]; sol1 = 
 sol - Join[dx . phi, dy . phi, 0  phi, 0  phi];

{U, V, W, P} = 
  Map[Interpolation@Join[sgrid, Transpose@List@#, 2] &, 
   Partition[sol1, Length[sgrid]]];

Visualization

Show[DensityPlot[W[x, y], {x, -2, 2}, {y, -.5, .5}, 
  PlotLegends -> Automatic, ColorFunction -> "Rainbow", 
  PlotRange -> Automatic, AspectRatio -> Automatic, 
  ColorFunctionScaling -> True], 
 StreamPlot[
  Evaluate[-Grad[W[x, y], {x, y}]], {x, -2, 2}, {y, -.5, .5}, 
  AspectRatio -> Automatic, StreamColorFunction -> None, 
  StreamStyle -> LightGray, PlotLegends -> Automatic, 
  PlotRange -> Automatic]]

Show[DensityPlot[W[x, y], {x, -2, 2}, {y, -.5, .5}, 
  AspectRatio -> Automatic, ColorFunction -> "Rainbow", 
  PlotLegends -> Automatic, PlotRange -> Automatic], 
 StreamPlot[-Evaluate[{U[x, y], V[x, y]}], {x, -2, 2}, {y, -.5, .5}, 
  AspectRatio -> Automatic, StreamColorFunction -> None, 
  StreamStyle -> LightGray, PlotLegends -> Automatic, 
  PlotRange -> Automatic]]