5
$\begingroup$

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];
$\endgroup$

2 Answers 2

5
$\begingroup$

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

$\endgroup$
4
  • $\begingroup$ 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. $\endgroup$
    – user82261
    Mar 15, 2023 at 6:46
  • $\begingroup$ @u Which version are you in? Anyway, try the editted code tested in 13.2.1. $\endgroup$
    – xzczd
    Mar 15, 2023 at 8:49
  • $\begingroup$ 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? $\endgroup$
    – user82261
    Mar 15, 2023 at 13:04
  • $\begingroup$ @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. $\endgroup$
    – xzczd
    Mar 17, 2023 at 4:33
3
$\begingroup$

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

$\endgroup$
9
  • $\begingroup$ 1. 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. $\endgroup$
    – xzczd
    Mar 15, 2023 at 10:13
  • $\begingroup$ Thanks, I modified my answer -Cos[x]. $\endgroup$ Mar 15, 2023 at 10:16
  • $\begingroup$ @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? $\endgroup$ Mar 15, 2023 at 10:41
  • $\begingroup$ 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. $\endgroup$
    – xzczd
    Mar 15, 2023 at 10:52
  • $\begingroup$ @xzczd I thought so too, but it doesn't work on v12.2 for this example! $\endgroup$ Mar 15, 2023 at 11:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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