Boundary conditions being ignored or not obeyed by NDSolveValue

Question

I have the following coupled ordinary differential equations (drift-diffusion with a loss term) that I am having trouble with solving. One of the boundary condition I supply is being ignored. When I provide a different one that should also hold, I flag an error.

mod1D[zMax_, {l0_, lb_, \[Gamma]s_}, {g_, \[Sigma]_}] := Module[{},
  
  

  \[Tau]s = 1/\[Gamma]s;
  
  (*dn/dt=0*)
  pde1 = -D[D[u[z], {z, 1}] + lb/l0 v[z], {z, 1}];
  (*ds/dt=0*)
  pde2 = -D[D[v[z], {z, 1}] + lb/l0 u[z], {z, 1}] + v[z] -
    ( 2 g \[Tau]s)/(Sqrt[2 \[Pi]] l0^2 \[Sigma])
      Exp[-(z^2/(2 \[Sigma]^2))];
  
  (* Neumann bc for charge*)
  Bu3 = NeumannValue[0, z == zMax];
  Bu4 = NeumannValue[0, z == -zMax];
  
  
  (* Dirichlet bc for spin*)
  R1 = DirichletCondition[v[z] == 0, z == -zMax];
  R2 = DirichletCondition[v[z] == 0, z == zMax];
  
  (* Dirichlet bc for charge*)
  R3 = DirichletCondition[u[z] == -u[-z], z == zMax];
  R4 = DirichletCondition[u[z] == 0, z == 0];
  
  (*Solve the PDE over the region*)
  sol = NDSolveValue[{pde1 == Bu3 + Bu4, pde2 == 0, R1, R2, 
     R3}, {u, v}, {z, -zMax, zMax}, 
    Method -> {"FiniteElement", 
      "MeshOptions" -> {"MaxCellMeasure" -> 0.0001}}, 
    MaxStepSize -> 0.005];

  ]

mod1D[10, {1, 2, 1}, {1, 0.001}]
Plot[sol[[1]][z], {z, -10, 10}, PlotStyle -> {{Black, Thick}}, 
 PlotRange -> {{-10, 10}, {-2, 2.1}}, 
 FrameLabel -> {"y", "z", "u"}, 
 PlotLabel -> "Steady-state solution for u"]

gives

which is not an odd function despite the boundary condition R3. If I use the boundary condition R4 instead, I get an error NDSolveValue::bcnop: No places were found on the boundary where z==0 was True, so DirichletCondition[u==0,z==0] will effectively be ignored. However, I know, from analytic solutions, solving same problem in Python using solve_bvp, and solving this same problem in Mathematica with an added dimension, that u[0] = 0.

What I expect to see if the same shape as shown above but shifted down to produce an odd function of z. I get same issue with NDSolve.

Using 13.2.

Answer 1

To have a boundary condition within the region you need to add an include point like so:

Needs["NDSolve`FEM`"]

mesh = ToElementMesh[Line[{{-zMax}, {zMax}}], 
   "IncludePoints" -> {{0}}, "MaxCellMeasure" -> 0.0001];
(*Solve the PDE over the region*)
sol = NDSolveValue[{pde1 == Bu3 + Bu4, pde2 == 0, R1, R2, R4}, {u, 
    v}, {z} \[Element] mesh];

To model a periodic boundary condition use PeriodicBoundaryCondition:

mod1D[zMax_, {l0_, lb_, \[Gamma]s_}, {g_, \[Sigma]_}] := 
 Module[{}, \[Tau]s = 1/\[Gamma]s;
  (*dn/dt=0*)pde1 = -D[D[u[z], {z, 1}] + lb/l0  v[z], {z, 1}];
  (*ds/dt=0*)
  pde2 = -D[D[v[z], {z, 1}] + lb/l0  u[z], {z, 1}] + 
    v[z] - (2  g  \[Tau]s)/(Sqrt[
         2  \[Pi]]  l0^2  \[Sigma])  Exp[-(z^2/(2  \[Sigma]^2))];
  (*Neumann bc for charge*)Bu3 = NeumannValue[0, z == zMax];
  Bu4 = NeumannValue[0, z == -zMax];
  (*Dirichlet bc for spin*)
  R1 = DirichletCondition[v[z] == 0, z == -zMax];
  R2 = DirichletCondition[v[z] == 0, z == zMax];
  (*Dirichlet bc for charge*)
  R3 = DirichletCondition[u[z] == -u[-z], z == zMax];
  R4 = DirichletCondition[u[z] == 0, z == 0];
  pbc1 = 
   PeriodicBoundaryCondition[-u[z], z == -zMax, 
    Function[z, z + 2 zMax]];
  pbc2 = 
   PeriodicBoundaryCondition[-u[z], z == zMax, 
    Function[z, z - 2 zMax]];
  mesh = 
   ToElementMesh[Line[{{-zMax}, {zMax}}], "IncludePoints" -> {{0}}, 
    "MaxCellMeasure" -> 0.0001];
  (*Solve the PDE over the region*)
  sol = NDSolveValue[{pde1 == Bu3 + Bu4, pde2 == 0, R1, R2, R4, pbc1, 
     pbc2}, {u, v}, {z} \[Element] mesh];
  ]

Needs["NDSolve`FEM`"]

mod1D[10, {1, 2, 1}, {1, 0.001}]
Plot[sol[[1]][z], {z, -10, 10}, PlotStyle -> {{Black, Thick}}, 
 PlotRange -> {{-10, 10}, {-2, 2.1}}, FrameLabel -> {"y", "z", "u"}, 
 PlotLabel -> "Steady-state solution for u"]