2
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ I recall a similar question where ElementMeshInterpolation was the bottleneck... $\endgroup$ – Henrik Schumacher Dec 13 '17 at 0:05
  • $\begingroup$ Is the integration domain still a square in your actual problem? $\endgroup$ – xzczd Dec 13 '17 at 5:46
  • $\begingroup$ No, my integration domain will be more complicated, that is why I am using finite elements. $\endgroup$ – user12588 Dec 14 '17 at 8:52
3
$\begingroup$

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).

$\endgroup$
  • $\begingroup$ Is there any further interest for an answer of the question ??? $\endgroup$ – Ulrich Neumann Dec 15 '17 at 17:27
  • $\begingroup$ Your answer looks great, the only catch is that I cannot control the error, and that is not good for my purposes. But thank you for your time! $\endgroup$ – user12588 Jan 4 '18 at 8:02
  • $\begingroup$ @user12588: Thank you for your feedback. error control would be the same as you used in your approach $\endgroup$ – Ulrich Neumann Jan 4 '18 at 8:20
5
$\begingroup$

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

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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