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]
Improve this question
$\endgroup$
12
  • $\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$ Commented Mar 22, 2017 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$ Commented Mar 22, 2017 at 22:06
  • $\begingroup$ I have posted a question related to your problem (very likely) $\endgroup$ Commented Mar 22, 2017 at 23:06
  • $\begingroup$ Mac, Mathematica 10.2.0.0 $\endgroup$ Commented Mar 22, 2017 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$ Commented Mar 22, 2017 at 23:33

3 Answers 3

Reset to default
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$
22
  • $\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$ Commented Mar 23, 2017 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$ Commented Mar 23, 2017 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$ Commented Mar 23, 2017 at 17:35
  • $\begingroup$ You can ignore the Neumann 0. The code should work as is else I'll have a look tomorrow. $\endgroup$ Commented Mar 23, 2017 at 18:06
  • $\begingroup$ @MaximLyutikov, could you figure this out? $\endgroup$ Commented Mar 24, 2017 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$
1
  • $\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$ Commented Mar 22, 2017 at 23:21

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.