Using NDSolve Iteratively to Solve a Boundary Value Problem

Question

I am trying to solve a simple nonlinear boundary value problem with FEM using the in built routine NDSolveValue. What I am doing is painfully slow, and I am sure you know how to do better. The equations I want to solve are way more complicated, so this one just serves as an illustration of what I want to do.

The equation in question takes the following form

$\partial^2_x\psi+\partial^2_{y}\psi=e^{\psi}$

with the integration domain being $(x,y)\in[0,1]\times[0,1]$. The boundary conditions are pretty simple:

$\partial_x \psi =0$ at $x=1$

and

$\psi=0$ at x =0, y=0 and y =1.

First, I create a grid:

<< NDSolve`FEM`
Om = Rectangle[];
Bmesh = ToBoundaryMesh[Om,MaxCellMeasure -> {"Length" -> 0.01}];
mesh = ToElementMesh[Om,MaxCellMeasure -> {"Length" -> 0.01}];
pts = mesh["Coordinates"];

Then I solve the linear problem iteratively:

u0[x_, y_] = 0;
u0data = Apply[u0, pts, 2];
Errmax = 10;
Errmin = 10^-6;
dE = (Errmax + Errmin)/2;
Print["Error = ", dE // Dynamic]
While[Errmin < dE < Errmax,
du = NDSolveValue[{
\!\(\*SuperscriptBox[\(du\), 
TagBox[
RowBox[{"(", 
RowBox[{"0", ",", "2"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] + 
\!\(\*SuperscriptBox[\(du\), 
TagBox[
RowBox[{"(", 
RowBox[{"2", ",", "0"}], ")"}],
Derivative],
MultilineFunction->None]\)[x, y] - E^u0[x, y] du[x, y] == -(\!\(
\*SubscriptBox[\(\[PartialD]\), \(x, x\)]\(u0[x, y]\)\) + \!\(
\*SubscriptBox[\(\[PartialD]\), \(y, y\)]\(u0[x, y]\)\) - 
      Exp[u0[x, y]]) + NeumannValue[0, x == 1], 
 DirichletCondition[
  du[x, y] == 0, ({x, y} \[Element] Bmesh) && x != 1]}, 
du, {x, y} \[Element] mesh, 
Method -> {"PDEDiscretization" -> {"FiniteElement"}}];
dudata = Apply[du, pts, 2];
u0data = u0data + dudata;
u0 = ElementMeshInterpolation[{mesh}, u0data];
dE = Norm[dudata, Infinity];
Clear[du, dudata];
];
ufinal = u0;
Clear[u0, Errmin, Errmax, dE]

The code is just to slow (I mean, I have implemented this in Fortran spectral collocation methods and using a simple Newton's method and runs in seconds). Is there a way to speed this up in Mathematica?

Answer 1

If I understand your question right, you try to find an iterative solution of your pde:

NestList[NDSolveValue[{Laplacian[\[Psi][x, y], {x, y}] ==Exp[#] (1 - # + \[Psi][x, y]) + NeumannValue[0, x == 1], \[Psi][0, y] == 0, \[Psi][x, 0] == 0, \[Psi][x, 1] ==0}, \[Psi][x, y], {x, 0, 1}, {y, 0, 1}] &, 0, 5]

the result (3 iterations) is obtained in 0.3s ( error 10^-9).

Answer 2

In version 12.0 this is straight forward:

sol = NDSolveValue[{Laplacian[u[x, y], {x, y}] - E^u[x, y] == 0, 
     DirichletCondition[u[x, y] == 0, x != 1]}, 
   u, Element[{x, y}, Rectangle[]]];
Plot3D[sol[x, y], Element[{x, y}, Rectangle[]]]