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}}}];

enter image description here 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!


2 Answers 2


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:

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]

enter image description here

Next, we set up the nonlinear problem:

pde2 = pde1 + 
   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}];


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.

enter image description here

  • $\begingroup$ Thanks a lot for your answer! The code now works for small values of f. However, I have been trying to slowly increase f (using the last solution as the initial seed following your suggestion) but the result does not change at all. I got up to f=0.3 and the solution shouldn't look exactly the same as for f=0.007. Do you know what is going wrong here? Thank you once again! $\endgroup$
    – lorenzop
    Jun 10, 2021 at 14:28
  • 1
    $\begingroup$ One or more of several things could be going on. Try to compare the solution of Alex and this one with the same f and check that they are the same. It's very well possible that the nonlinear solver is not strong enough to solve the equation, the non-linearity is too strong. $\endgroup$
    – user21
    Jun 11, 2021 at 3:51

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]];

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

ff[0][x_, y_] := -2.9411764705882355 x; ff[-1][x_, y_] := -2.9411764705882355 x;
Do[ff[i] = 
     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}]

Figure 1

  • $\begingroup$ (+1), you may want to add a note on where you got the ff[0] from. $\endgroup$
    – user21
    Jun 10, 2021 at 7:16
  • 1
    $\begingroup$ @user21 Thank you for your remarks. The definition of initial condition ff[0] has been added. $\endgroup$ Jun 10, 2021 at 8:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.