5

I am trying to solve the following PDE: $$ u_{xx} + u_{yy} = \begin{cases} - \cos(x) \quad -\pi/2 \leq x \leq \pi/2, \\ 0 \quad \text{otherwise} \end{cases} $$

The domain is $\Omega = [-\pi,\pi] \times [0,2]$. The boundary conditions are periodic at $x = \pm \pi$ - $u(-\pi,y) = u(\pi,y)$ - and $u(x,0) = 0$ (Dirchlet) and $u_y(x,2) = 0$ (Neumann).

I am able to solve the PDE if I choose all of conditions to be Dirichlet. I cannot, however, solve the PDE in the case where I have a combination of Neumann, Dirichlet, Periodic boundary conditions. I have no idea how to implement the Neumann boundary condition appropriately. Also my implementation of periodic boundary conditions gives errors.

Ω = Rectangle[{-Pi, 0}, {Pi, 2}];

op = Laplacian[u[x, y], {x, y}] - Piecewise[{{-Cos[x], -Pi/2 < x < Pi/2}, {0, x < -Pi/2}, {0, x > Pi/2}}, {x, -Pi, Pi}];

Subscript[Γ, D] = {DirichletCondition[u[x, y] == 0, y == 0], DirichletCondition[u[x, y] == 0, y == 2 ]};

mapping = Last[FindGeometricTransform[{{-Pi, 0}, {-Pi, 2}}, {{Pi, 0}, {Pi, 2}}]];

Subscript[Γ, P] = PeriodicBoundaryCondition[u[x, y], x == -Pi, mapping];

ufun = NDSolveValue[{op == 0, Subscript[Γ, D], Subscript[Γ, P]}, u, {x, y} ∈ Ω];

ContourPlot[ufun[x, y], {x, y} ∈ Ω, ColorFunction -> "Temperature", AspectRatio -> Automatic];

user21
  • 39,710
  • 8
  • 110
  • 167
user82261
  • 151
  • 2

2 Answers2

5

4 issues here, 2 are simple mistakes, 2 are not.

  1. Though it happens not to cause problem in this case, your understanding for syntax of Piecewise is wrong, please check its document carefully. The correct one should be Piecewise[{{-Cos[x], -Pi/2 < x < Pi/2}}].

  2. You've inverted the direction in FindGeometricTransform.

  3. Currently PeriodicBoundaryCondition and DirichletCondition cannot overlap even at points, personally I think this design is a bit unnecessary and inconvenient.

  4. There's an issue related to periodic b.c. of FiniteElement method, which is discussed in this post. To obtain the desired solution, we need to set 2 PeriodicBoundaryConditions and a ghost layer.

To sum up:

eps = 0.1;
Ω = Rectangle[{-Pi, 0}, {Pi + eps, 2}];
offset = 1.5;
op = Laplacian[u[x, y], {x, y}] - 
   Piecewise[{{-Cos[x], -Pi/2 + offset < x < Pi/2 + offset}}];

Subscript[Γ, D] = DirichletCondition[u[x, y] == 0, y == 0];

mapping = Last[FindGeometricTransform[{{Pi, 0}, {Pi, 2}}, {{-Pi, 0}, {-Pi, 2}}]]; mappinginverse = Last[FindGeometricTransform[{{-Pi, 0}, {-Pi, 2}}, {{Pi, 0}, {Pi, 2}}]]; Subscript[Γ, P] = {PeriodicBoundaryCondition[u[x, y], x == -Pi && 0 < y < 2, mapping], PeriodicBoundaryCondition[u[x, y], x == Pi + eps && 0 < y < 2, mappinginverse]};

ufun = NDSolveValue[{op == 0, Subscript[Γ, D], Subscript[Γ, P]}, u, {x, y} ∈ Ω]; ufunwrong = NDSolveValue[{op == 0, Subscript[Γ, D], Last@Subscript[Γ, P]}, u, {x, y} ∈ Ω];

With[{y = 1.5}, Plot[{ufun[Mod[x, 2 Pi, -Pi], y], ufunwrong[Mod[x, 2 Pi, -Pi], y]} // Evaluate, {x, -Pi, 3 Pi}, PlotStyle -> {Automatic, Dashed}]]

enter image description here

Notice I've set offset = 1.5 to illustrate the necessity of double periodic b.c., for your original problem, set offset = 0.

"Wait, where's the Neumann b.c.?" Your Neumann b.c. corresponds to a zero NeumannValue. As mentioned in Details section of document of NeumannValue:

...not specifying a boundary condition at all is equivalent to specifying a Neumann 0 condition.

So you don't need to write code for your Neumann b.c..

xzczd
  • 65,995
  • 9
  • 163
  • 468
  • When I run your code, I get the error: NDSolveValue::femnodpbc: DirichletCondition can not be present on the target boundary of a PeriodicBoundaryConditon. Seems like the same issue persists. – user82261 Mar 15 '23 at 06:46
  • @u Which version are you in? Anyway, try the editted code tested in 13.2.1. – xzczd Mar 15 '23 at 08:49
  • I am running version 12, hmm. Let me see if I can find an update. If not, what seems to be going on? Why do I get the error? – user82261 Mar 15 '23 at 13:04
  • @user82261 Have you tested the updated code in v12? If it still doesn't work, I'd say it's a bug of earlier version. – xzczd Mar 17 '23 at 04:33
3

modified

Thanks to @xzczd for the helpful discussion

Here my simplified version without ghost layer and mapping:

The PeriodicBoundaryCondition's have to be defined on both sides x==-Pi and x==Pi. y==0 must be excluded, otherwise there is a conflict with DirichletCondition.

Using a triangle-mesh NDSolve is able to correctly handle the periodic bc's

\[CapitalOmega] = Rectangle[{-Pi, 0}, {Pi , 2}];
offset = 1.5;
pde = Laplacian[u[x, y], {x,y }] == -Cos[x] Boole[-Pi/2 + offset < x< Pi/2 + offset];
\[CapitalGamma]D = DirichletCondition[u[x, y] == 0, y == 0];

[CapitalGamma]P = { PeriodicBoundaryCondition[u[x, y],x == [Pi] && 0 < y <= 2, TranslationTransform[{-2 Pi, 0}]], PeriodicBoundaryCondition[u[x, y], x == -[Pi] && 0 < y <= 2, TranslationTransform[{2 Pi, 0}]]}

ufun = NDSolveValue[{pde, [CapitalGamma]D, [CapitalGamma]P},u, {x, y} [Element] [CapitalOmega], Method -> {"FiniteElement","MeshOptions" -> {"MeshElementType" ->"TriangleElement", "MeshOrder" -> 1}}]

Plot3D[ufun[x, y], Element[{x, y}, [CapitalOmega]],AxesLabel -> {x, y, u[x, y]}, MeshFunctions -> {#2 &}]

enter image description here

The periodic boundary conditions are well fullfilled now

Plot [Evaluate[Table[ufun[x, y], {y, Subdivide[0, 2, 10]}]], {x, -Pi,Pi}]

enter image description here

The solution QP asked for follows with offset==0

Ulrich Neumann
  • 53,729
  • 2
  • 23
  • 55
  • It should be -Cos[x]. 2. One of your PeriodicBoundaryCondition is actually ignored, try pde = Laplacian[u[x, y], {x, y}] == -Cos[x] Boole[-Pi/2 + offset < x < Pi/2 + offset]; and then check the derivative as shown my answer. The ghost layer cannot be eliminated.
  • – xzczd Mar 15 '23 at 10:13
  • Thanks, I modified my answer -Cos[x]. – Ulrich Neumann Mar 15 '23 at 10:16
  • @xzczd In the link "ghost layer" you mentioned there was a hint: "For triangle-mesh epsilon==0 is possible". But here in this problem not? – Ulrich Neumann Mar 15 '23 at 10:41
  • That's another way to go, but you need to set "MeshElementType" -> TriangleElement. (This isn't the default setting. ) BTW, you forgot to remove the eps in your code. – xzczd Mar 15 '23 at 10:52
  • @xzczd I thought so too, but it doesn't work on v12.2 for this example! – Ulrich Neumann Mar 15 '23 at 11:04
  • How do you set it? – xzczd Mar 15 '23 at 11:19
  • I tried ufun = NDSolveValue[{pde , \[CapitalGamma]D , \[CapitalGamma]P }, u, {x, y} \[Element] \[CapitalOmega], Method -> {"FiniteElement", "MeshOptions" -> {"MeshElementType" -> "TriangleElement", "MeshOrder" -> 1}} – Ulrich Neumann Mar 15 '23 at 11:44
  • .. and got error NDSolveValue::bcnop: No places were found on the boundary where x==\[Pi]&&0<y<=2 was True, so PeriodicBoundaryCondition[u,x==\[Pi]&&0<y<=2,TransformationFunction[{{1,0,-2 \[Pi]},{0,1,0},{0,0,1}}]] will effectively be ignored. – Ulrich Neumann Mar 15 '23 at 11:45
  • I can reproduce the issue in v13.2.1. Instead of 0 < y <= 2, using 0 < y < 2 seems to circumvent it, try it out in v12.2 :) . – xzczd Mar 15 '23 at 12:11