Unable to solve nonlinear PDE with NDSolve

Question

Lately, I've been trying to solve the following PDE: \begin{equation} -v_0 |\nabla F| + {\bf f}\cdot \nabla F +D\nabla^2F = -1 \end{equation} inside a 2D region between two disks both centered in the origin with radii $r=1$ and $l=5$ respectively. The boundary conditions are $F=0$ on the inner circle and $\hat{n}\cdot\nabla F=0$ on the external one. Here is the code I have written:

v0 = 1.; D1 = 0.01; f = 0.7;
r2 := 1; l2 := 5; cell2 := 0.001;
\[CapitalOmega]2 = 
  RegionDifference[Disk[{0, 0}, l2], Disk[{0, 0}, r2]];
pde2 = D1 Laplacian[FF[x, y], {x, y}] + f D[FF[x, y], x] - 
   v0 Sqrt[D[FF[x, y], x]^2 + D[FF[x, y], y]^2] ;
dcond2 = DirichletCondition[FF[x, y] == 0, x^2 + y^2 == r2^2];

Fsol = NDSolveValue[{pde2 == -1 + NeumannValue[0., x^2 + y^2 == l2^2],
     dcond2}, {FF[x, y]}, {x, y} \[Element] \[CapitalOmega]2, 
   Method -> {"PDEDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> {"MaxCellMeasure" -> cell2}}}];

However, as you can see I am getting some error messages but I don't know why nor how to deal with them. I actually think this is coming from the sqrt in the $|\nabla F|$ term, but I can't get rid of that. Hope it is clear enough and thanks in advance for your help!

Answer 1

This is not a complete answer but a start; specifically the value for the convection term f is quite large, but we will get to that later.

The definitions:

v0 = 1; D1 = 0.01; f = 0.7;
r2 = 1; l2 = 5; cell2 = 0.002;
dcond2 = DirichletCondition[FF[x, y] == 0, x^2 + y^2 == r2^2];

The mesh:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[Annulus[{0, 0}, {r2, l2}], 
  MaxCellMeasure -> cell2]

We start with the linear problem:

vars = {FF[x, y], {x, y}};
f = 0.007;

pde1 = DiffusionPDETerm[vars, D1] + 
  ConvectionPDETerm[vars, -{f, 0}] - 1
Fsol1 = NDSolveValue[{pde1 == 0, dcond2}, 
  vars[[1]], {x, y} \[Element] mesh]
Plot3D[Fsol1, {x, y} \[Element] mesh, PlotRange -> All]

Next, we set up the nonlinear problem:

pde2 = pde1 + 
  SourcePDETerm[vars, 
   v0 Sqrt[Total[Grad[FF[x, y], {x, y}]^2] + $MachineEpsilon]]

Solve it. Note, that we use the solution of the linear problem as an initial seed for the nonlinear solver.

Fsol2 = NDSolveValue[{pde2 == 0, dcond2}, 
   vars[[1]], {x, y} \[Element] mesh, 
   InitialSeeding -> {FF[x, y] == Fsol1}];

Visualize:

Plot3D[Fsol2, {x, y} \[Element] mesh, PlotRange -> All]

This is now your starting point. Slowly increase f and use the solution from the last nonlinear problem as an initial seed.

Answer 2

First we map region to the anulus of unit outer radius and make high quality mesh. Then we use standard algorithm of the false transient (in analogue to the diffusion equation) as follows

v0 = 1; D1 = 0.01; f = 0.7;
r2 = 1/5; l2 = 5; 
reg = RegionDifference[Disk[{0, 0}, 1], Disk[{0, 0}, r2]];
Needs["NDSolve`FEM`"]

mesh0 = ToElementMesh[reg, MaxCellMeasure -> .0001] 



ff[0][x_, y_] := -2.9411764705882355 x; ff[-1][x_, y_] := -2.9411764705882355 x;
Do[ff[i] = 
   NDSolveValue[{ 
     D1/l2^2 Laplacian[FF[x, y], {x, y}] + f/l2 D[FF[x, y], x] - 
       v0/l2/2 (Sqrt[
           D[ff[i - 1][x, y], x]^2 + D[ff[i - 1][x, y], y]^2] + 
          Sqrt[D[ff[i - 2][x, y], x]^2 + D[ff[i - 2][x, y], y]^2]) + 
       1 - 20 (FF[x, y] - ff[i - 1][x, y]) == 0, 
     DirichletCondition[FF[x, y] == 0, x^2 + y^2 == r2^2]}, FF, 
    Element[{x, y}, mesh0], 
    Method -> {"PDEDiscretization" -> {"FiniteElement"}}];, {i, 1, 6}]

As initial condition we use exact solution $F=c x$ where c is solution of equation

NSolve[f/l2 c - v0/l2 Abs[c] + 1 == 0, c] // Quiet

 {{c -> -2.94118}}

The iterations converge very fast and finally we have

Table[DensityPlot[ff[i][x, y], Element[{x, y}, mesh0], 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  PlotRange -> All], {i, 1, 5}]