FEM breaks because of Abs(f'(x))

Question

I am trying to solve a 1d differential equation in Mathematica using the finite element method. I have successfully implemented the problem in a commercial software called COMSOL, but I like to use Mathematica much more due to the flexibility.

I quickly realised that a Abs[f'[x]] function causes the DiscretizePDE to break. My question: Does it have to do with the fact that the derivative of Abs[] is not defined? I tried to define my own function myAbs[x]=If[x<0, -x, x] but that does not seem to help. Has anybody an idea?

If that helps I can also provide a small working example that shows the error.

The differential equation is a Poisson Boltzmann equation which is already dimensionless and has the form:

$$\nabla (\epsilon(r) \nabla \phi(r)) + \kappa ( n_a(r)-n_c(r)) = 0$$

where $\kappa = \frac{e_0^2}{k_B T \epsilon_0 a0}\approx 1,3*10^4$ is some constant from non-dimensionalization.

The space dependent dielectric constant is

$$\epsilon(r)=\epsilon_{op}+\kappa*p^2*n_s(r)*L(p*\vert\nabla \phi\vert)$$

$$L(A)= \frac{\coth(A)}{A}-\frac{1}{A^2}$$

where $p\approx 2$ and $\epsilon_{op}=6$.

The particle densities are given by Boltzmann factors: $$n_{a/c}(r)= e^{\pm \phi(r)- W(r)}$$

and

$$n_s(r)=e^{-W(r)\cdot}\frac{\sinh(p*\vert\nabla\phi\vert)}{p*\vert\nabla\phi\vert}$$

where the term W(r) is simply a term which is very large when $x<1$ which prevents the particles to come very close to the surface.

The problem is that the function L(r) and the solvent density ns(r) depends on the absolute value of the electric field Abs[\phi'[x]] which I replaced by the myAbs function which throws the error:" The FEMStiffnessElements operator failed". Since I know that in this problem $\phi'[x]$ is always negative, I just replaced it by -$\phi'[x]$ which gives me a beautiful result. I don't understand why the if expression breaks the code since the case $\phi'[x]>0$ never happens.

The parameters are

a0 = 5.29*10^(-11);(*[m] Bohr radius*)
nref = a0^(-3);
NA = 6.02*10^23;(*[1/mol] Avogadro number*)

rsol = 2.56*1*10^(-10);  (*m,diameter of water,taken as the lattice \
size*)
rc = 1.38*1*10^(-10) + rsol; (*m,radius of cations,K+*)
ra = 3*1*10^(-10);(* m,radius of anions,1.81 for Cl-,2.2 for I-*)
cbulk = 100;(* bulk concentration,0.1 M*)
chiv = 0.05; (*density fraction of vacancy in solution bulk*)

\[Gamma]a = (2 ra/rsol)^3; (*relative size of anions*)
\[Gamma]c = (2 rc/rsol)^3; (*relative size of anions*)
\[Gamma]s = 1;
nsbulk = 
  55.6*10^3*
   NA/nref; (*dimensionless number density of free solvent (water)*)
nabulk = 
  cbulk*NA/nref; (*dimensionless bulk concentration of anions,0.1 M;*)
ncbulk = 
  cbulk*NA/nref; (*dimensionless bulk concentration of cations,0.1 M;*)
nmax = (\[Gamma]a*nabulk + \[Gamma]c*ncbulk + \[Gamma]s*nsbulk)/(1 - 
     chiv); (*dimensionless lattice density in solution phase;*)
chia = nabulk/nmax ;(*density fraction of anions in solution bulk*)
chic = ncbulk/nmax ;(*density fraction of cations in solution bulk*)
chis = nsbulk/nmax ;(*density fraction of solvents in solution bulk*)

systemSize = 170
epsopS = 
  6.0; (*relative optical dielectric permittivity in solution phase*)
psd = 1.86673; (*dimensionless dipole moment*)
kappa = 13296;

One helper function is the step function which also defines the shortRangeRepulsion from the metal:

StepFunction[x_, d_] := 
  Piecewise[{{0, x <= -d}, {1, 
     x >= d}, {0.5 + 0.9375*x/d - 0.625*(x/d)^3 + 0.1875*(x/d)^5, 
     True}}];
Plot[StepFunction[x, 1], {x, -2, 2}]


shortRangeRepulsionParameter = 10;(*short range forces,nd w.r.t.kB*T*)
\[Sigma]1 = 4;(*distance in short range force relations*)
\[Sigma]2 = (\[Sigma]1/
    5);  (*nd distance in short range force relations*)


shortRangeRepulsion[distanceMetal_] := 
 shortRangeRepulsionParameter*(\[Sigma]1/\[Sigma]2)^6*
   StepFunction[\[Sigma]1 - distanceMetal, 10^-3] + 
  shortRangeRepulsionParameter*(\[Sigma]1/(\[Sigma]2 + 
     Abs[distanceMetal - \[Sigma]1]))^6 StepFunction[
    distanceMetal - \[Sigma]1, 10^-3]
Plot[shortRangeRepulsion[distanceFromMetal], {distanceFromMetal, -10, 
  10}, FrameLabel -> {"distance From Metal [A]", "Potential"}, 
 Frame -> True, PlotRange -> Full]

The problem needs an initial guess:

b = 0.4;
\[Phi]Ini = (kappa*0.1)/(
   epsopS*8*
    b^2)*((1 + 2*b^2*x^2)*(1 - Erf[b*x]) - (2*b*x)/Sqrt[Pi]*
      Exp[-b^2*x^2]);
Plot[\[Phi]Ini, {x, 0, 100}, PlotRange -> All]

And the DGL is simply:

Ecr = 1*10^(-3);
Elb = 1*10^(-6);

Efield = Min[If[\[Phi]'[x] < 0, -\[Phi]'[x], \[Phi][x]] + Elb, 50] (*Working but wrong*);


\[Theta]a = 
  Min[Exp[Min[+\[Phi][x], 30] - shortRangeRepulsion[x]], 10^20];
\[Theta]c = 
  Min[Exp[Min[-\[Phi][x], 30] - shortRangeRepulsion[x]], 10^20];
\[Theta]s = 
  Exp[-shortRangeRepulsion[
      x]] (Min[Sinh[psd*Efield]/(psd*Efield), 10^20] StepFunction[
       psd*Efield - Ecr, Elb] + (1 + (psd*Efield)^2/6) StepFunction[
       Ecr - psd*Efield, Elb]);
\[CapitalOmega] = \[Gamma]a*chia*\[Theta]a + \[Gamma]c*
    chic*\[Theta]c + \[Gamma]s*chis*\[Theta]s + chiv;


na = nmax*chia*\[Theta]a/\[CapitalOmega]; nc = 
 nmax*chic*\[Theta]c/\[CapitalOmega]; ns = 
 nmax*chis*\[Theta]s/\[CapitalOmega];

Plot[#, {x, 0, 100}, 
   PlotRange -> {0, 1.2}] & /@ ({na/nabulk, nc/ncbulk, 
    ns/nsbulk} /. {\[Phi][x] -> \[Phi]Ini, \[Phi]'[x] -> 
     D[\[Phi]Ini, x]})

A = ComplexExpand[Norm[\[Phi]'[x]]]*psd + Elb; (*not working but what I want*)
A = If[\[Phi]'[x] < 0, -\[Phi]'[x], \[Phi][x]]*psd + Elb;(*working*)
A = -\[Phi]'[x]*psd + Elb(*working*);


langevinFunction = (Coth[A]/A - 1/A^2)*StepFunction[A - Ecr, 100] + 
   1/3*StepFunction[Ecr - A, 100];


epsoEff = epsopS + kappa*psd^2*ns*langevinFunction;
testPlot[
   x_] = (epsoEff /. {\[Phi][x] -> \[Phi]Ini, \[Phi]'[x] -> 
      D[\[Phi]Ini, x]});
Plot[%, {x, 0, 100}, PlotRange -> All]


potentialDiffusionTerm = DiffusionPDETerm[{\[Phi][x], {x}}, {epsoEff}];
potentialSourceTerm = +kappa*(na - nc);

ode\[Phi] = potentialDiffusionTerm + potentialSourceTerm;

The mesh is :

coordintates = {{0.}, {systemSize}};
bmesh = ToBoundaryMesh["Coordinates" -> coordintates, 
  "BoundaryElements" -> {PointElement[{{1}, {2}}]}]

mesh = ToElementMesh[bmesh, MaxCellMeasure -> .2]

And the solving is done via:

Timing[sol = 
  NDSolveValue[{ode\[Phi] == 0, 
    DirichletCondition[\[Phi][x] == 188, x == 0],
    DirichletCondition[\[Phi][x] == 0, x >= systemSize]}, {\[Phi]}, 
   x \[Element] mesh, InitialSeeding -> {\[Phi][x] == \[Phi]Ini}, 
   Method -> FiniteElements]]
Plot[sol[[1]][x], {x, 0, systemSize}, PlotRange -> All]

This code works for the two last A's in the code.

Answer 1

In Mathematica the function RealAbs[] is an absolute value function where it assumes the argument to be real. With this function one can write the absolute value of the electric field as

If[RealAbs[ϕ'[x]] < Elb, Elb, RealAbs[ϕ'[x]]]

where the cutoff Elb is introduced to remove the singularity in the function L(A) at A=0.

The point is that the derivative of RealAbs can be computed

In[1]:= D[RealAbs[x], x]

Out[1]= x/RealAbs[x]

While the derivative of Abs

In[2]:= D[Abs[x], x]

Out[2]= Derivative[1][Abs][x]

is not computed.