5
$\begingroup$

I need to solve non-linear Poisson equation

Laplacian[u[x, y], {x, y}] == u[x, y]^2

Over a non-rectangular domain

The problem in short: non-linear Poisson equation over rectangular domain runs OK, and linear Poisson equation over non-rectangular domain runs OK, but not the non-linear over non-rectangular.

The domain is

boundaries = {-y, .25^2 - (x)^2 - y^2, -x, y - 1, x - 1};

\[CapitalOmega]in = 
ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}];

Show[RegionPlot[\[CapitalOmega]in], 
ContourPlot[
Evaluate[Thread[boundaries == 0]], {x, 0., 1}, {y, 0, 1.}, 
ContourStyle -> {Purple, Green, Red, Blue, Purple}], 
PlotRange -> {{0.0, 1}, {0., 1.}}, AspectRatio -> Automatic] 

with simple boundary conditions

  Conditions = {DirichletCondition[u[t, x, y] == 1, 
  boundaries[[1]] == 0.],
  DirichletCondition[u[t, x, y] == 1, boundaries[[2]] == 0],
  DirichletCondition[u[t, x, y] == 1, boundaries[[3]] == 0.],
  DirichletCondition[u[t, x, y] == 1, boundaries[[4]] == 0.],
  DirichletCondition[u[t, x, y] == 1, boundaries[[5]] == 0.],
  u[0, x, y] == 1};

I try to run a relaxation scheme

 Eq = Laplacian[u[t, x, y], {x, y}] - u[t, x, y]^2

sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
Conditions}, u, {t, 0, 1}, {x, y} \[Element] \[CapitalOmega]in, 
Method -> {"MethodOfLines", Method -> "Automatic", 
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

The problem is "Nonlinear coefficients are not supported in this version of NDSolve".

The linear Poisson equation runs OK

Eq = Laplacian[u[t, x, y], {x, y}] - u[t, x, y]

sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
 Conditions}, u, {t, 0, 1}, {x, y} \[Element] \[CapitalOmega]in, 
Method -> {"MethodOfLines", Method -> "Automatic", 
"DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

ContourPlot[sol[1, x, y], {x, y} \[Element] \[CapitalOmega]in, 
 ColorFunction -> "TemperatureMap", Contours -> 50, 
 AspectRatio -> Automatic]

Also, non-linear over a rectangular domain runs OK:

boundaries = {-y, -x, y - 1, x - 1};
\[CapitalOmega]in = 
ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}];

 Show[RegionPlot[\[CapitalOmega]in], 
 ContourPlot[
 Evaluate[Thread[boundaries == 0]], {x, 0., 1}, {y, 0, 1.}, 
 ContourStyle -> {Purple, Green, Red, Blue, Purple}], 
 PlotRange -> {{0.0, 1}, {0., 1.}}, AspectRatio -> Automatic] 

Conditions = {DirichletCondition[u[t, x, y] == 1, 
 boundaries[[1]] == 0.],
 DirichletCondition[u[t, x, y] == 1, boundaries[[2]] == 0],
 DirichletCondition[u[t, x, y] == 1, boundaries[[3]] == 0.],
 DirichletCondition[u[t, x, y] == 1, boundaries[[4]] == 0.],
 u[0, x, y] == 1};

 sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
  Conditions}, u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, 
 Method -> {"MethodOfLines", Method -> "Automatic", 
 "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

 ContourPlot[sol[1, x, y], {x, y} \[Element] \[CapitalOmega]in, 
 ColorFunction -> "TemperatureMap", Contours -> 50, 
 AspectRatio -> Automatic]
$\endgroup$
  • $\begingroup$ Your last block of code doesn't work. I guess it's because you have changed the definiton of Conditions . It works fine with Conditions = {u[t, 0, y] == 1,u[t, 1, y] == 1,u[t, x, 0] == 1,u[t, x, 1] == 1,u[0, x, y] == 1};. Please correct your question. $\endgroup$ – andre314 Mar 22 '17 at 21:43
  • $\begingroup$ @andre - yes the boundaries and Conditions are a bit different in the last block from the previous - there 4 boundaries and 4 Conditions (not 5 like before). I think your and mine are identical. If you run the last block as a whole it works OK. $\endgroup$ – Maxim Lyutikov Mar 22 '17 at 22:06
  • $\begingroup$ I have posted a question related to your problem (very likely) $\endgroup$ – andre314 Mar 22 '17 at 23:06
  • $\begingroup$ Mac, Mathematica 10.2.0.0 $\endgroup$ – Maxim Lyutikov Mar 22 '17 at 23:29
  • $\begingroup$ It's Strange. I haved tried my code below on the cloud (Wolfram Development Paltform, Mathematica 11.1 Unix). It doesn't work too. $\endgroup$ – andre314 Mar 22 '17 at 23:33
6
$\begingroup$

If I understand the question correctly, you are looking for the stationary solution of the non-linear PDE over a region. Unfortunately, NDSolve can not handle this out of the box (V11.1) but you can use the low-level FEM functions to get a solution. To do so we roughly follow the idea presented here. The idea is to create an interpolating function in every non-linear step and feed that into the linearized coefficients.

Setup: We use a linear Laplace equation, so get the linear discretization:

Needs["NDSolve`FEM`"]
\[CapitalOmega]in = 
  ImplicitRegion[
   And @@ (# <= 0 & /@ {-y, .25^2 - (x)^2 - y^2, -x, y - 1, 
       x - 1}), {x, y}];
op = -Laplacian[u[x, y], {x, y}] == 0;
conditions = {DirichletCondition[u[x, y] == 1, True]};

ProcessPDEEquations is a build in function (V11. in NDSolveFEM context) that does the same as the PDEtoMatrix function in the linked code.

{dPDE, dBC, vd, sd, md} = 
  ProcessPDEEquations[{op, conditions}, 
   u, {x, y} \[Element] \[CapitalOmega]in, 
   Method -> {"FiniteElement", 
     "MeshOptions" -> {"ImproveBoundaryPosition" -> False}}];

mesh2 = md["ElementMesh"];
linearLoad = dPDE["LoadVector"];
linearStiffness = dPDE["StiffnessMatrix"];
diriPos = dBC["DirichletRows"];

Now, we can start to write the non-linear loop. This version is a little different to the version in the other post. Mainly, how it handles boundary conditions.

The InitializedPDECoefficients gets the linearized coefficients at uOld. Finding this linearization is the hardest part. For those I'd like to refer you to Wolfgang Bangerth's Video lectures.

ClearAll[rhs]
rhs[uIn_] := Module[{uOld}, uOld = uIn;
  Do[ClearAll[u0];
   u0 = ElementMeshInterpolation[{mesh2}, uOld];
   nlPdeCoeff = InitializePDECoefficients[vd, sd
     , "LoadCoefficients" -> {
       {-u0[x, y]^2}
       }
     , "LoadDerivativeCoefficients" -> {
       {{-Derivative[1, 0][u0][x, y], -Derivative[0, 1][u0][x, y]}}
       }
     , "ReactionCoefficients" -> {
       {2 u0[x, y]}
       }
     ];
   nlsys = DiscretizePDE[nlPdeCoeff, md, sd];
   nlLoad = nlsys["LoadVector"];
   nlStiffness = nlsys["StiffnessMatrix"];
   ns = nlStiffness + linearStiffness;
   nl = nlLoad + linearLoad;

   nl[[diriPos]] = {0.};
   ns[[diriPos, All]] = 0.;
   ns[[All, diriPos]] = 0.;
   (ns[[#, #]] = 1.) & /@ diriPos;

   dU = LinearSolve[ns, nl];
   Print[i, " Residual: ", Norm[nl, Infinity], "  Correction: ", 
    Norm[dU, Infinity]];
   uOld = uOld + dU;
   , {i, 8}];
  uOld]

Another word on boundary conditions. In this version of the non-linear loop the initial guess has to satisfy the boundary conditions. Because the uOld satisfies the boundary conditions we modify the load vector and (tangent) stiffness matrix such that uOld = uOld + du will still respect the boundary conditions. (Again this is explained in the video lectures)

We setup the initial guess and set the boundary conditions:

guess = 1.;
uOld = ConstantArray[{guess}, md["DegreesOfFreedom"]];
uOld[[diriPos]] = dBC["DirichletValues"];

We run the code:

uNew = rhs[uOld];

1 Residual: 0.0016115  Correction: 0.0637695
2 Residual: 5.34982*10^-6  Correction: 0.000162486
3 Residual: 3.45194*10^-11  Correction: 8.79701*10^-10
4 Residual: 4.94049*10^-15  Correction: 1.27027*10^-15
5 Residual: 2.65066*10^-15  Correction: 1.02935*10^-15
6 Residual: 2.84928*10^-15  Correction: 1.21563*10^-15
7 Residual: 3.17975*10^-15  Correction: 1.16816*10^-15
8 Residual: 2.75214*10^-15  Correction: 1.10233*10^-15

We have quadratic convergence. (That will change if you omit the "ImproveBoundaryPosition" -> False - then region will be approximated better but the interpolation of the uOld on curved leads to some damping(?))

Well, the plot:

eif = ElementMeshInterpolation[{md["ElementMesh"]}, uNew];
Plot3D[eif[x, y], {x, y} \[Element] \[CapitalOmega]in, 
 PlotRange -> All]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ I am adapting your code to my problem - how do I modify guess = 1.; to more complicated boundary conditions with more complicated Dirichlet conditions and some Neumann? $\endgroup$ – Maxim Lyutikov Mar 23 '17 at 15:58
  • $\begingroup$ This line uOld[[diriPos]] = dBC["DirichletValues"]; is setting the boundary conditions to what you have. guess=1 is to the the initial guess inside the domain. - With the NeumannValue I don't know. You should create an example to which you know the solution and experiment (and tell me ;-) $\endgroup$ – user21 Mar 23 '17 at 16:17
  • $\begingroup$ My actual boundary conditions look something like that: boundaries = {-y, 0.25^2 - (x)^2 - y^2, -x, y - 1, x - 1}; op = -Laplacian[u[x, y], {x, y}] + NeumannValue[0, boundaries[[1]] == 0.] + NeumannValue[0, boundaries[[5]] == 0.] == 0; conditions = { DirichletCondition[u[x, y] == 0, boundaries[[4]] == 0.], DirichletCondition[u[x, y] == 0, boundaries[[3]] == 0.], DirichletCondition[u[x, y] == 1/(1 + x), boundaries[[2]] == 0.]}; I would like generate 'guess' array to be consistent with that. How would I do that? $\endgroup$ – Maxim Lyutikov Mar 23 '17 at 17:35
  • $\begingroup$ You can ignore the Neumann 0. The code should work as is else I'll have a look tomorrow. $\endgroup$ – user21 Mar 23 '17 at 18:06
  • $\begingroup$ @MaximLyutikov, could you figure this out? $\endgroup$ – user21 Mar 24 '17 at 6:52
4
$\begingroup$

In version 12.0 this works out of the box:

boundaries = {-y, .25^2 - (x)^2 - y^2, -x, y - 1, x - 1};
\[CapitalOmega]in = 
  ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}];
Conditions = {DirichletCondition[u[t, x, y] == 1, 
    boundaries[[1]] == 0.], 
   DirichletCondition[u[t, x, y] == 1, boundaries[[2]] == 0], 
   DirichletCondition[u[t, x, y] == 1, boundaries[[3]] == 0.], 
   DirichletCondition[u[t, x, y] == 1, boundaries[[4]] == 0.], 
   DirichletCondition[u[t, x, y] == 1, boundaries[[5]] == 0.], 
   u[0, x, y] == 1};

Nonlinear equation:

Eq = Laplacian[u[t, x, y], {x, y}] - u[t, x, y]^2;

Solve:

sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
   Conditions}, u, {t, 0, 1}, {x, y} \[Element] \[CapitalOmega]in]

Plots:

ContourPlot[sol[1, x, y], {x, y} \[Element] \[CapitalOmega]in, 
 ColorFunction -> "TemperatureMap", Contours -> 50, 
 AspectRatio -> Automatic]

enter image description here

Plot3D[sol[1, x, y], {x, y} \[Element] \[CapitalOmega]in, 
 PlotRange -> All]

enter image description here

| improve this answer | |
$\endgroup$
1
$\begingroup$

Not an answer

I don't undestand your last comment.
If I take your last block of code with the addition of the missing definition Eq = Laplacian[u[t, x, y], {x, y}] - u[t, x, y]^2, it doesn't work on my machine (Mma 11.0.0.0 , Windows 7) . The first error message is :

NDSolveValue::femnonlinear: Nonlinear coefficients are not supported in this version of NDSolve.

Here is the whole code :

Eq = Laplacian[u[t, x, y], {x, y}] - u[t, x, y]^2

boundaries = {-y, -x, y - 1, x - 1};
\[CapitalOmega]in = 
ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y}];

 Show[RegionPlot[\[CapitalOmega]in], 
 ContourPlot[
 Evaluate[Thread[boundaries == 0]], {x, 0., 1}, {y, 0, 1.}, 
 ContourStyle -> {Purple, Green, Red, Blue, Purple}], 
 PlotRange -> {{0.0, 1}, {0., 1.}}, AspectRatio -> Automatic] 

Conditions = {DirichletCondition[u[t, x, y] == 1, 
 boundaries[[1]] == 0.],
 DirichletCondition[u[t, x, y] == 1, boundaries[[2]] == 0],
 DirichletCondition[u[t, x, y] == 1, boundaries[[3]] == 0.],
 DirichletCondition[u[t, x, y] == 1, boundaries[[4]] == 0.],
 u[0, x, y] == 1};

 sol = NDSolveValue[{Eq == Derivative[1, 0, 0][u][t, x, y], 
  Conditions}, u, {t, 0, 1}, {x, 0, 1}, {y, 0, 1}, 
 Method -> {"MethodOfLines", Method -> "Automatic", 
 "DifferentiateBoundaryConditions" -> {True, "ScaleFactor" -> 1}}]

 ContourPlot[sol[1, x, y], {x, y} \[Element] \[CapitalOmega]in, 
 ColorFunction -> "TemperatureMap", Contours -> 50, 
 AspectRatio -> Automatic]
| improve this answer | |
$\endgroup$
  • $\begingroup$ on my machine (Mac, Mathematica 10.2.0.0) both my and your code work OK, but only for rectangular grid. Yes, the general problem is the non-linearity. $\endgroup$ – Maxim Lyutikov Mar 22 '17 at 23:21

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.