3
$\begingroup$

The documentation for PeriodicBoundaryCondition (https://reference.wolfram.com/language/ref/PeriodicBoundaryCondition.html) has: enter image description here Where it says $u ( x_{target} ) = a + b\ u ( f ( x_{target} ) )$, I think it should instead say $u( f ( x_{target} ) ) = a + b\ u ( x_{target} )$. I could be wrong but I believe this is demonstrated by this example: https://wolfram.com/xid/0bswu24h9fy656tmxe-jnf5k3. I have copied the code here, only modifying it from $a=-1/20$ to $a=0$ because this will demonstrate what I am talking about:

Ω = Rectangle[{0, 0}, {2, 1}];
pde = -\!\(
\*SubsuperscriptBox[\(∇\), \({x, y}\), \(2\)]\(u[x, y]\)\) == 
   If[1.25 <= x <= 1.75 && 0.25 <= y <= 0.5, 1., 0.];
Subscript[Γ, D] = 
  DirichletCondition[u[x, y] == 0, (y == 0 || y == 1) && 0 < x <= 2];
a = 0; b = 2;
pbc = PeriodicBoundaryCondition[a + b*u[x, y], x == 0 && 0 <= y <= 1, 
   TranslationTransform[{2, 0}]];
ufun = NDSolveValue[{pde, pbc, Subscript[Γ, D]}, 
   u, {x, y} ∈ Ω];
ContourPlot[ufun[x, y], {x, y} ∈ Ω, 
 ColorFunction -> "TemperatureMap", AspectRatio -> Automatic]

The left-hand edge, from $(0,0)$ to $(0,1)$, is the target (i.e. where the predicate in PeriodicBoundaryCondition is true) and the right-hand edge, from $(2,0)$ to $(2,1)$, is the source because $x_{source} = f (x_{target})$. Now as $b=2$ (and $a=0$) by the current definition we would expect $u(x_{target}) = 2 u(x_{source})$, meaning the values at the left-hand edge should be twice as large as those on the right-hand edge. But they're not. Instead, they are half as large. This implies that the definition should instead be $u( f ( x_{target} ) ) = a + b u ( x_{target} )$. I have explored this for many hours with many examples and keep arriving at the same conclusion.

$\endgroup$
  • $\begingroup$ There are some subtle issues with the PeriodicBoundaryCondition. A triangle mesh and a symmetrized PBC, may be required. You can look at my answer 223465 as well as the other answers. $\endgroup$ – Tim Laska Aug 10 at 15:59
  • $\begingroup$ I think you analysis is correct - this is a typo in the documentaaion and I have updated the documentation. Sorry for the the trouble and thanks for reporting this. You always have the option to report things like this to support AT wolfram.com. I may not see all issues I am responsible for if posted here. $\endgroup$ – user21 Aug 14 at 9:29
  • $\begingroup$ @user21 Thanks for that. Could you post your comment as an answer so I may select it as the best answer? Also I went through the "contact us" form on wolfram.com on the 6th of July, case 4565672, but as I received no update on it for some time I felt the need to post it here. Is the contact us form what I should be using for issues such as this? Thanks again. $\endgroup$ – Anthony Aug 15 at 16:35
  • $\begingroup$ @Anthony, oh well.... not sure why that report did not reach me yet; in principal this is a correct approach to contact support. $\endgroup$ – user21 Aug 17 at 4:23
  • $\begingroup$ @user21 No worries; thank you. I notice the documentation is still the same - does it take a while for your changes to take effect? $\endgroup$ – Anthony Aug 18 at 20:08
1
$\begingroup$

I think you analysis is correct - this is a typo in the documentation and I have updated the documentation. Sorry for the the trouble and thanks for reporting this. You always have the option to report things like this to support AT wolfram.com. I may not see all issues I am responsible for if posted here.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you. I notice the documentation is still the same - does it take a while to update? $\endgroup$ – Anthony Aug 25 at 12:56
  • $\begingroup$ @Anthony, well yes, it will take a new release for it to show :-( $\endgroup$ – user21 Aug 25 at 13:04
3
$\begingroup$

As suggested in the comments and elaborated in my answer 223465, you can use a triangle mesh and symmetrize the PeriodicBoundaryCondition by making the following workflow:

Create Triangle Mesh

Here we use ToElementMesh to create a triangle mesh with refinement on the boundaries.

Needs["NDSolve`FEM`"]
Ω = Rectangle[{0, 0}, {2, 1}];
(* Create Triangle Mesh *)
mesh = ToElementMesh[Ω, 
   "MaxCellMeasure" -> {"Length" -> 0.05}, 
   "MaxBoundaryCellMeasure" -> 0.0025, 
   "MeshElementType" -> TriangleElement];

Create a Plotting Function

Here we will create a function that simulates a parametric function, generates a contour plot, and an error plot of the periodic condition of the two boundaries.

plotFn[a_, b_][pfun_] := 
 Module[{ufun, uRange, legendBar, options, cp, error, assoc},
  ufun = pfun[a, b];
  uRange = MinMax[ufun["ValuesOnGrid"]];
legendBar = 
   BarLegend[{"TemperatureMap", uRange}, 50, 
    LegendLabel -> Style["u", Opacity[0.6`]]];
options = {PlotRange -> uRange, 
    ColorFunction -> ColorData[{"TemperatureMap", uRange}], 
    ContourStyle -> Opacity[0.1`], ColorFunctionScaling -> False, 
    Contours -> 30, PlotPoints -> All, FrameLabel -> {"x", "y"}, 
    PlotLabel -> 
     Style[StringTemplate["u(x,y) Field for a=`` and b=`` "][a, b], 
      18], AspectRatio -> Automatic, ImageSize -> 500};
  cp = Legended[
    ContourPlot[ufun[x, y], {x, y} \[Element] ufun["ElementMesh"], 
     Evaluate[options]], legendBar];
  cp = Rasterize@cp;
  error = 
   Plot[{a + b*ufun[0, y] - ufun[2, y]}, {y, 0, 1}, PlotPoints -> 200,
     PlotRange -> 1.*^-15 {-1, 1}];
  assoc = <|"cp" -> cp, "error" -> error|>
  ]

Set up ParametricNDSolveValue

It would be nice to look at the effect of the $a$ and $b$ parameters. So, let's use ParametricNDSolveValue to generate a parametric function so we can quickly test the parameters.

pde = -Laplacian[u[x, y], {x, y}] == 
      If[1.25 <= x <= 1.75 && 0.25 <= y <= 0.5, 1., 0.];
ΓD = 
    DirichletCondition[u[x, y] == 0, (y == 0 || y == 1) && 0 < x < 2];
(* Symmetrized PBCs *)
pbcf = PeriodicBoundaryCondition[a + b*u[x, y], x == 0 && 0 <= y <= 1, 
      TranslationTransform[{2, 0}]];
pbcr = PeriodicBoundaryCondition[-a /b + 1/b*u[x, y], 
   x == 2 && 0 <= y <= 1, 
      TranslationTransform[{-2, 0}]];
pfun = ParametricNDSolveValue[{pde, pbcf, pbcr, ΓD}, 
     u, {x, y} ∈ mesh, {a, b}]

Test Several $a$ and $b$ values

sim01 = plotFn[0, 1][pfun]
sim02 = plotFn[0, 2][pfun]
sim03 = plotFn[1/10, 1][pfun]
sim04 = plotFn[1/10, 2][pfun]

Simulation Series

The error between the left and right side is quite low. Introducing the offset parameter $a$ causes some ringing at the corner points. The ringing is most likely due to an inconsistency with the DirichletCondition, DC, and the PeriodicBoundaryCondition, PBC. The DC specifies zero on the top and bottom boundaries, but the PBC specifies an offset between the left and right boundary. At the corner points, there is a discontinuity between the DC and PBC.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Hi Tim, thanks for the response. So according to your other answer creating a triangle mesh is good because it doesn't require "extra padding of the domain"? Also what do you mean by "symmetrize" the Periodic boundary conditions? $\endgroup$ – Anthony Aug 10 at 18:40
  • $\begingroup$ @Anthony In theory, the solver should be able to accommodate triangles or quads, but there is something subtle going on that causes quads to fail without padding. By "symmetrize", I mean that you need to apply the PBC twice switching the source and target. I tried to indicate that with "f" and "r" for forward and reverse. When one does both, there is a good correspondence with other FEM solvers. $\endgroup$ – Tim Laska Aug 10 at 18:54
  • 1
    $\begingroup$ Thanks Tim. For symmetry though should it not read ``` pbcr = PeriodicBoundaryCondition[-a + (1/b)*u[x, y], x == 2 && 0 <= y <= 1, TranslationTransform[{-2, 0}]]; ``` ? Also as far as I can tell it is still not applying $u ( x_{target} ) = a + b\ u ( x_{source} )$ as the definition suggests but $u( x_{source} ) = a + b\ u ( x_{target} )$ ? Finally, I'm assuming we must symmetrize and use a triangle mesh when using FindGeometricTransform as well? (as in this question: mathematica.stackexchange.com/questions/228305/…) $\endgroup$ – Anthony Aug 10 at 19:20
  • $\begingroup$ @Anthony For the {a,b}={0,1} case it does not matter. In the general case, you are probably right except that it should read -a/b + (1/b)u[x,y]. The previous answer only considered the {a,b}={0,1}, so I do not know how general the solution is. It would be good to have some verification cases to test. With regard to your other question, I have used NDEigensystem much. Hopefully, someone can chime in that has more experience. $\endgroup$ – Tim Laska Aug 11 at 0:54
  • $\begingroup$ thanks very much for the thorough edit to your answer. To me your sims still show the definition to be incorrect though. When $a=0$ and $b=2$ for example, pbcf has the target at $x=0$ and $0\le y \le 1$, and the source at $x=2$ and $0\le y \le 1$. We see on your graph that $u(x_{source})=2u(x_{target})$, which is opposite to what the definition of PeriodicBoundaryCondition would suggest (namely that $u(x_{target})=u(x_{source})$. Your pbcr also shows opposite to the definition, as it shows $u(x_{source})=(1/2)u(x_{target})$ with the source and target at $x=0$ and $x=2$ respectively. $\endgroup$ – Anthony Aug 11 at 14:35

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.