2
$\begingroup$

Question: I want to solve a PDE using NDEigensystem. The region is a rectangle Rectangle[{0,-70}, {20, 70}]. How can I make a mesh that satisfies the following two conditions?

(1) The mesh is finer near $y \approx 0$.

(2) The line $y=0$ should be the boundary of the mesh.

Motivation: The PDE that I want to solve can be written in the following command: First, I set the mesh as below:

mesh = ToElementMesh[Rectangle[{0, -70}, {20, 70}],  
  MeshRefinementFunction -> 
   Function[{vertices, area}, 
    Block[ {x, y}, {x, y} = Mean[vertices]; 
     If[-10 < y < 10, area > 0.1, area > 10]]] ]

and then write down the PDE:

NDEigensystem[{{3*psi3[x, y] - 10*psi2[x, y]*Sign[y] - Derivative[0, 1][psi2][x, y] - 
    I*Derivative[1, 0][psi2][x, y], 3*psi4[x, y] - 10*psi1[x, y]*Sign[y] + Derivative[0, 1][psi1][x, y] - 
    I*Derivative[1, 0][psi1][x, y], 3*psi1[x, y] - 10*psi4[x, y]*Sign[y] - Derivative[0, 1][psi4][x, y] - 
    I*Derivative[1, 0][psi4][x, y], 3*psi2[x, y] - 10*psi3[x, y]*Sign[y] + Derivative[0, 1][psi3][x, y] - 
    I*Derivative[1, 0][psi3][x, y]}, PeriodicBoundaryCondition[psi1[x, y], x == 0, 
   TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi2[x, y], x == 0, 
   TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi3[x, y], x == 0, 
   TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], PeriodicBoundaryCondition[psi4[x, y], x == 0, 
   TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], DirichletCondition[psi1[x, y] == 0, 
   0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi2[x, y] == 0, 
   0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi3[x, y] == 0, 
   0 < x < 20 && (y == 70 || y == -70)], DirichletCondition[psi4[x, y] == 0, 
   0 < x < 20 && (y == 70 || y == -70)]}, {psi1, psi2, psi3, psi4}, Element[{x, y}, mesh], 10]

The single important fact about the PDE is that the differential operator contains $\mathrm{sgn}(y)$, hence it is not continuous at $y=0$. This is why I want the mesh boundary to include $y=0$. The mesh I created is visualized as follows:

enter image description here

As wanted, the mesh is finer. However, $y=0$ is not the boundary of the mesh.

The result of one solution obtained by solving the PDE is plotted as below:

enter image description here

In the above plot, the horizontal direction is the $x$-direction. Note that the solution is localized near $y\approx 0$. This is why I want to make the mesh finer near $y\approx 0$. I don't believe the result is the good representation of the solution, since the graph looks jagged, which I don't expect since the differential operator has a translation symmetry in $x$-direction.

$\endgroup$
1
  • $\begingroup$ Line y=0 is not shown in the mesh. What boundary condition supposed to be on the line y=0? $\endgroup$ Mar 20 at 21:46
3
$\begingroup$

You can achieve your objectives with an anisotropic quad mesh that you construct manually. There are no diffusive terms in your operators, so NDEigensystem will complain about the system being convection dominated. Therefore, you may want to take a closer look at your system of equations and see if there is an opportunity to add diffusive terms.

Helper functions

The following helper functions will help construct an anisotropic Quad mesh (note that not all functions will be used).

(*Import required FEM package*)
Needs["NDSolve`FEM`"];
(*Define Some Helper Functions For Structured Meshes*)
pointsToMesh[data_] := 
  MeshRegion[Transpose[{data}], 
   Line@Table[{i, i + 1}, {i, Length[data] - 1}]];
unitMeshGrowth[n_, r_] := 
 Table[(r^(j/(-1 + n)) - 1.)/(r - 1.), {j, 0, n - 1}]
meshGrowth[x0_, xf_, n_, r_] := (xf - x0) unitMeshGrowth[n, r] + x0
firstElmHeight[x0_, xf_, n_, r_] := 
 Abs@First@Differences@meshGrowth[x0, xf, n, r]
lastElmHeight[x0_, xf_, n_, r_] := 
 Abs@Last@Differences@meshGrowth[x0, xf, n, r]
findGrowthRate[x0_, xf_, n_, fElm_] :=(*Quiet@*)
 Abs@FindRoot[
    firstElmHeight[x0, xf, n, r] - fElm, {r, 0.00000001, 
     100000000/fElm}, Method -> "Brent"][[1, 2]]
meshGrowthByElm[x0_, xf_, n_, fElm_] := 
 N@Sort@Chop@meshGrowth[x0, xf, n, findGrowthRate[x0, xf, n, fElm]]
meshGrowthByElm0[len_, n_, fElm_] := meshGrowthByElm[0, len, n, fElm]
flipSegment[l_] := (#1 - #2) & @@ {First[#], #} &@Reverse[l];
leftSegmentGrowth[len_, n_, fElm_] := meshGrowthByElm0[len, n, fElm]
rightSegmentGrowth[len_, n_, fElm_] := 
 Module[{seg}, seg = leftSegmentGrowth[len, n, fElm];
  flipSegment[seg]]
reflectRight[pts_] := 
 With[{rt = ReflectionTransform[{1}, {Last@pts}]}, 
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
reflectLeft[pts_] := 
 With[{rt = ReflectionTransform[{-1}, {First@pts}]}, 
  Union[pts, Flatten[rt /@ Partition[pts, 1]]]]
extendMesh[mesh_, newmesh_] := Union[mesh, Max@mesh + newmesh]

Construct anisotropic Quad mesh

The following workflow will create a symmetric graded mesh about $y=0$ and a uniform mesh in the x-direction.

ht = 140;
len = 20;
ny = 100;
nx = 100;
Print["Mesh along the X-direction"]
rh = pointsToMesh@Subdivide[0, len, nx]
Print["Mesh along the Y-direction"]
rv = pointsToMesh@reflectLeft@leftSegmentGrowth[ht/2, ny/2, ht/2/ny/10]
(*Create 2D mesh with RegionProduct*)
rp = RegionProduct[rh, rv];
(* Build mesh based on region product *)
crd = MeshCoordinates[rp];
inc = Delete[0] /@ MeshCells[rp, 2];
mesh = ToElementMesh["Coordinates" -> crd, 
   "MeshElements" -> {QuadElement[inc]}];
Print["Mesh zoomed in around y=0"]
Show[mesh["Wireframe"], PlotRange -> {{0, 20}, {-5, 5}}]

Mesh

This mesh is refined and has a line about $y=0$.

Solve and plot solution

(*Solve and plot*)
{vals, funs} = 
  NDEigensystem[{{3*psi3[x, y] - 10*psi2[x, y]*Sign[y] - 
      Derivative[0, 1][psi2][x, y] - I*Derivative[1, 0][psi2][x, y], 
     3*psi4[x, y] - 10*psi1[x, y]*Sign[y] + 
      Derivative[0, 1][psi1][x, y] - I*Derivative[1, 0][psi1][x, y], 
     3*psi1[x, y] - 10*psi4[x, y]*Sign[y] - 
      Derivative[0, 1][psi4][x, y] - I*Derivative[1, 0][psi4][x, y], 
     3*psi2[x, y] - 10*psi3[x, y]*Sign[y] + 
      Derivative[0, 1][psi3][x, y] - I*Derivative[1, 0][psi3][x, y]}, 
    PeriodicBoundaryCondition[psi1[x, y], x == 0, 
     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], 
    PeriodicBoundaryCondition[psi2[x, y], x == 0, 
     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], 
    PeriodicBoundaryCondition[psi3[x, y], x == 0, 
     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], 
    PeriodicBoundaryCondition[psi4[x, y], x == 0, 
     TransformationFunction[{{1, 0, 20}, {0, 1, 0}, {0, 0, 1}}]], 
    DirichletCondition[psi1[x, y] == 0, 
     0 < x < 20 && (y == 70 || y == -70)], 
    DirichletCondition[psi2[x, y] == 0, 
     0 < x < 20 && (y == 70 || y == -70)], 
    DirichletCondition[psi3[x, y] == 0, 
     0 < x < 20 && (y == 70 || y == -70)], 
    DirichletCondition[psi4[x, y] == 0, 
     0 < x < 20 && (y == 70 || y == -70)]}, {psi1, psi2, psi3, psi4}, 
   Element[{x, y}, mesh], 10];
Table[Plot3D[
  funs[[i, 1]][x, y] Conjugate[funs[[i, 1]][x, y]], {x, y} \[Element] 
   mesh, PlotRange -> All, PlotLabel -> Chop@vals[[i]], Axes -> True, 
  AxesLabel -> Automatic], {i, Length[vals]}]

Solution

The solution no longer is jagged in the X-direction. However, there are many oscillations in the Y-direction. You may need to consider diffusive terms to dampen these oscillations.

$\endgroup$
2
  • $\begingroup$ Thank you so much for your detailed answer. Your answer does answer my question, but it is strange that the solution is not localized near $y\approx 0$. There is a good physical reason that the solution should be localized near $y\approx 0$. I may suspect that this is due to the "convection dominated PDE". How is the result reliable? $\endgroup$
    – eigenvalue
    Mar 21 at 8:32
  • $\begingroup$ @eigenvalue The result may not be reliable. There is a good chance, but not necessarily, that a convection-dominated problem is numerically unstable. I found a question regarding the semi-classical Pauli operator here 148252 that seems to have some similarity to yours, but it has a diffusive term. If there is an opportunity to add this diffusive term to your operator, it could stabilize the results. Alternatively, you could investigate adding artificial diffusion. $\endgroup$
    – Tim Laska
    Mar 21 at 17:59

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.