Solutions of NDSolve with FEM don't solve the PDE's

Question

I have 2 mixed-order PDE's that I need to solve in a finite computational space. The first equation is $$\frac{2 A \mu ^2 e^{-\frac{3 r^2}{\text{r0}^2}} \left(r^2 e^{\frac{2 r^2}{\text{r0}^2}}-\frac{4 A^2 \lambda }{\left(u(r,\theta ,\phi )+\frac{M}{2 r}+1\right)^{12}}\right)}{r^3}+\partial_r \text{Ef}(r,\theta ,\phi ) + 2\text{Ef}(r,\theta ,\phi )/r=0 $$ and the second is $$ \frac{576 a^2 M^2 r \sin ^2(\theta )}{(2 r u(r,\theta ,\phi )+M+2 r)^7}+\frac{6 \pi A^4 \lambda \mu ^2 e^{-\frac{4 r^2}{\text{r0}^2}}}{r^4 \left(u(r,\theta ,\phi )+\frac{M}{2 r}+1\right)^{19}}-\frac{128 \pi A^2 \mu ^2 r^5 e^{-\frac{2 r^2}{\text{r0}^2}}}{(2 r u(r,\theta ,\phi )+M+2 r)^7}-\frac{\pi \text{Ef}(r)^2}{\left(u(r,\theta ,\phi )+\frac{M}{2 r}+1\right)^7}+\nabla^2 u(r,\theta,\phi) = 0$$ The first PDE is a single derivative in the radial direction, and the second is a spherical laplacian operator is a complicated source, which depends on the scalar field.

The boundary conditions are two Neumann boundary conditions for Ef and u: $$ \partial_r Ef(r,\theta, \phi) = \frac{-2Ef(r, \theta, \phi)}{r} $$ $$ \partial_r u(r,\theta, \phi) = \frac{1 - u(r, \theta, \phi)}{r} $$

Everything else is a constant.

What I've tried:

I tried just using NDSOlve but I got the error The maximum derivative order of the nonlinear PDE coefficients for the Finite Element Method is larger than 1. It may help to rewrite the PDE in inactive form.

So I wrote it in inactive form as

ClearAll[r, \[Theta], \[Phi], \[Phi]f, Ef, A, \[Mu], \[Lambda], r0, E1]
\[Phi]f[r_] := A/r*Exp[-(r/r0)^2];
Subscript[\[CapitalPsi], bl] = 1 + M/(2*r);
Evec = {E1[r, \[Theta], \[Phi]], 0, 0};
EDiv = (2 E1[r, \[Theta], \[Phi]])/r + 
   Inactive[Derivative][1, 0, 0][E1][r, \[Theta], \[Phi]];
ULaplacian = (
   Csc[\[Theta]]^2 Inactive[Derivative][0, 0, 2][u][
     r, \[Theta], \[Phi]])/r^2 + (
   Cot[\[Theta]] Inactive[Derivative][0, 1, 0][u][
     r, \[Theta], \[Phi]])/r^2 + 
   Inactive[Derivative][0, 2, 0][u][r, \[Theta], \[Phi]]/r^2 + (
   2 Inactive[Derivative][1, 0, 0][u][r, \[Theta], \[Phi]])/r + 
   Inactive[Derivative][2, 0, 0][u][r, \[Theta], \[Phi]];
If[Simplify[
   Activate[ULaplacian] == 
    Laplacian[u[r, \[Theta], \[Phi]], {r, \[Theta], \[Phi]}, 
     "Spherical"]] != True, Print["Error. Check derivatives"]]
If[Simplify[
   Activate[EDiv] == Div[Evec, {r, \[Theta], \[Phi]}, "Spherical"]] !=
   True, Print["Error. Check derivatives"]]
EQ1 = (( 
     EDiv + 2*\[Mu]^2*\[Phi]f[r] - 
      8*\[Mu]^2*\[Lambda]*(Subscript[\[CapitalPsi], bl] + 
         u[r, \[Theta], \[Phi]])^-12*\[Phi]f[r]^3)) // Simplify;
EQ2 = (ULaplacian + 
    1/8 (Subscript[\[CapitalPsi], bl] + 
       u[r, \[Theta], \[Phi]])^-7*36*(a*M)^2/
     r^8 (r^2*Sin[\[Theta]]^2) - 
    2*\[Pi] (1/
        2 (Subscript[\[CapitalPsi], bl] + u[r, \[Theta], \[Phi]])^-7*
        Evec . Evec + \[Mu]^2/
        2*(Subscript[\[CapitalPsi], bl] + 
          u[r, \[Theta], \[Phi]])^-7*\[Phi]f[r]^2 - 
       3*\[Lambda]*\[Mu]^2 (Subscript[\[CapitalPsi], bl] + 
          u[r, \[Theta], \[Phi]])^-19*\[Phi]f[r]^4));
EIC = DirichletCondition[E1[r, \[Theta], \[Phi]] == 0, True];
uIC = NeumannValue[(1 - u[r, \[Theta], \[Phi]])/r, True];
EQs = Flatten[{EQ1, EQ2}];

Then solved it with

\[Epsilon] = 10^-6;
rmax = 30;
Block[{A = 1, r0 = 2, \[Mu] = 1, \[Lambda] = 0, M = 1, a = 0},
 sol = NDSolve[{Activate[EQ1] == EIC, Activate[EQ2] == uIC},
      {E1, u},
     {r, \[Epsilon], 
      rmax}, {\[Theta], \[Epsilon], \[Pi] - \[Epsilon]}, {\[Phi], \
\[Epsilon], 2 \[Pi] - \[Epsilon]},
     AccuracyGoal -> 20,
     PrecisionGoal -> 20,
     WorkingPrecision -> 35, 
     Method ->
      {"PDEDiscretization" ->
        {"FiniteElement", 
         
         "MeshOptions" -> {"MaxCellMeasure" -> 0.1, 
           "MeshOrder" -> 2},
         "IntegrationOrder" -> 5
         }
       }] // First;
 ]
E1sol = E1 /. sol;
Usol = u /. sol;

This finally produces solutions, however the solutions dont solve the original differential equations. If I plot the differential equation for the E1 field, I get: but this should be zero. I dont even know why there is a discontinuity at r~0.46. This doesnt match any other relevant scale in the problem. I tried increasing the working precision but I get the warning For the method FiniteElement, only machine-precision code is available. The original PDE still fails to be solved. I set $\lambda =0$ to decouple the PDE's, but the result is the same.

Update

So I managed to solve the system by transforming the differential operators into Inactive form. I had trouble with this initially since the Laplacian needs to be in spherical coordinates, but I managed to transform it into a totally Inactive form by

ULaplacian = (Cot[\[Theta]]  {0, 1, 0} . 
   Inactive[Grad][
    u[r, \[Theta], \[Phi]], {r, \[Theta], \[Phi]}])/r^2 + (
 2  {1, 0, 0} . 
   Inactive[Grad][u[r, \[Theta], \[Phi]], {r, \[Theta], \[Phi]}])/r + 
 Inactive[
   Div][(-{{-1, 0, 0}, {0, -(1/r^2), 0}, {0, 
       0, -(Csc[\[Theta]]^2/r^2)}} . 
     Inactive[Grad][
      u[r, \[Theta], \[Phi]], {r, \[Theta], \[Phi]}]), {r, \[Theta], \
\[Phi]}]

I checked this is indeed the operator I want by Simplify[ Activate[ULaplacian] == Laplacian[u[r, \[Theta], \[Phi]], {r, \[Theta], \[Phi]}, "Spherical"]]. Which returns True.

I then use some lower level functions in NDSolve to solve the system of equations using the code

EIC = NeumannValue[(-2 E1[r, \[Theta], \[Phi]])/r, True];
uIC = NeumannValue[(1 - u[r, \[Theta], \[Phi]])/r, True];
Block[{A = 1, r0 = 1, \[Mu] = 1, \[Lambda] = 0, M = 1, a = 0},
 {state} = NDSolve`ProcessEquations[{EQ1 == EIC, EQ2 == uIC},
    {E1, u},
   {r, \[Epsilon], rmax}, {\[Theta], 0, \[Pi]}, {\[Phi], 0, 
    2 \[Pi]}];
 NDSolve`Iterate[state];
 sols = NDSolve`ProcessSolutions[state];
 Ediff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[EQ1] /. sols]];
 Udiff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[EQ2] /. sols]];
 E1sol = E1 /. sols;
 Usol = u /. sols;
 ]

However, this still produces the unphysical discontinuities. Heres a plot of E1 .

I then found this post and this answer, which inspired me to create my own ElementGrid for the PDE solver.

I define the Functions

meshGrowth[x0_, xf_, n_, ratio_] := 
 Module[{k, fac, delta}, k = Log[ratio]/(n - 1);
  fac = Exp[k];
  delta = (xf - x0)/Sum[fac^(i - 1), {i, 1, n - 1}];
  N[{x0}~
    Join~(x0 + 
      delta Rest@
        FoldList[(#1 + #2) &, 0, 
         PowerRange[fac^0, fac^(n - 3), fac]])~Join~{xf}]]
PointsToMesh[data_] := 
 MeshRegion[Transpose[{data}], 
  Line@Table[{i, i + 1}, {i, Length[data] - 1}]]

and create my grid for NDSolve with

\[Epsilon] = 10^-6;
rmax = 10; \[Theta]max = \[Pi] - \[Epsilon]; \[Phi]max = 
 2 \[Pi] - \[Epsilon];
nrpoints = 200; n\[Theta]points = 20; n\[Phi]points = 20;
rratio = 7;
rpoints = meshGrowth[0.1, rmax, nrpoints, rratio];
\[Theta]points = Subdivide[0, \[Theta]max, n\[Theta]points];
\[Phi]points = Subdivide[0, \[Phi]max, n\[Phi]points];
rmesh = PointsToMesh[rpoints]
\[Theta]mesh = PointsToMesh[\[Theta]points]
\[Phi]mesh = PointsToMesh[\[Phi]points]
fullmesh = RegionProduct[rmesh, \[Theta]mesh, \[Phi]mesh];
HighlightMesh[fullmesh, Style[1, Brown]];
meshcoords = MeshCoordinates[fullmesh];
elementmesh = ToElementMesh[meshcoords];
groups = elementmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["TemperatureMap"][#] & /@ temp;
elementmesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]

This generates a uniform grid for $\theta, \phi$ and a nonlinear grid with increased density for small radius for the radial coordinates

When I execute the solver with

Block[{A = 1, r0 = 1, \[Mu] = 1, \[Lambda] = 0, M = 1, a = 0},
 {state} = NDSolve`ProcessEquations[{EQ1 == EIC, EQ2 == uIC},
    {E1, u},
   {r, \[Theta], \[Phi]} \[Element] elementmesh];
 NDSolve`Iterate[state];
 sols = NDSolve`ProcessSolutions[state];
 Ediff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[EQ1] /. sols]];
 Udiff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[EQ2] /. sols]];
 E1sol = E1 /. sols;
 Usol = u /. sols;
 ]

I get a much more reasonable solution. Heres a plot of E1

however, this still doesnt solve the original PDE satisfactorily. Heres a plot of the PDE for $E1$

It has pretty poor convergence for small radius.

How can I increase the convergence of the solver without running out of my laptops 16gb of RAM?

Update 2

I'm beginning to converge on something useful, though still not there yet. I ditched the ElementMesh approach above and instead relied on NDSolves internal Mesh generator. The reason I ditched the above approach is because the mesh I had to generate to get good convergence resulted in memory errors, since my laptop has only 16gb of ram.

I also make use the the InitialSeed Option in NDSOlve.

Choosing good initial seeds evidently is important. To do so, I make use of the knowledge of the physics and use the asymptote of the u function, which is equal to unity. This is reflected by the Neumann boundary condition.

I first solve the differential equation for E1 with the replacements u-> 1. Then compile the result for quick evaluation

EQInitEq = EQ1 /. u -> Function[{r, \[Theta], \[Phi]}, 1]
Esol = DSolve[EQInitEq == 0, E1, {r, \[Theta], \[Phi]}] /. 
     C[1][\[Theta], \[Phi]] -> 0 // First // Simplify;
E1InitGuess = 
  Compile[{r, \[Theta], \[Phi], \[Mu], M, A, r0, a, \[Lambda]}, 
   Evaluate[ E1[r, \[Theta], \[Phi]] /. Esol // Activate], 
   CompilationTarget -> "C"];

Evidently, the general solution for E1 contains an integral over the function u, so I cant use the general formula for E1 and instead rely on the asymptotic form. Together with $u(r, \theta, \phi)=1$, this constitutes my seed functions for NDSolve.

I compute the PDE's with

Block[{A = 1, r0 = 1, \[Mu] = 1, \[Lambda] = 0, M = 1, a = 0.},
 E1seed[r_?NumericQ, \[Theta]_?NumericQ, \[Phi]_?NumericQ] := 
  E1InitGuess[r, \[Theta], \[Phi], \[Mu], M, A, r0, a, \[Lambda]];
 useed[r_?NumericQ, \[Theta]_?NumericQ, \[Phi]_?NumericQ] := 1;
 
 {state} = NDSolve[{EQ1 == EIC, EQ2 == uIC},
    {E1, u},
   {r, \[Epsilon], 30}, {\[Theta], 0, \[Pi]}, {\[Phi], 0, 2 \[Pi]},
   Method -> {"FiniteElement"},
   InitialSeeding -> {E1[r, \[Theta], \[Phi]] == 
      E1seed[r, \[Theta], \[Phi]], 
     u[r, \[Theta], \[Phi]] == useed[r, \[Theta], \[Phi]]}
   ];
 Ediff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[EQ1] /. state]];
 Udiff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[EQ2] /. state]];
 E1sol = E1 /. state;
 Usol = u /. state;
 ]

and plot the result with

Block[{\[Theta]v = \[Pi]/2, \[Phi]v = \[Pi], pmin = \[Epsilon], 
  pmax = 30},
 Column[{
   Row[{LogLinearPlot[E1sol[r, \[Theta]v, \[Phi]v], {r, pmin, pmax}, 
      ImageSize -> Medium, PlotRange -> All], 
     LogLinearPlot[Usol[r, \[Theta]v, \[Phi]v], {r, pmin, pmax}, 
      ImageSize -> Medium, PlotRange -> All]}],
   Row[{LogLinearPlot[Ediff[r, \[Theta]v, \[Phi]v], {r, pmin, pmax}, 
      PlotRange -> All, ImageSize -> Medium], 
     LogLinearPlot[Udiff[r, \[Theta]v, \[Phi]v], {r, pmin, pmax}, 
      PlotRange -> All, ImageSize -> Medium]
     }]
   }]
 ]

This produces the results:

I get something even more reasonable than before, however the errors are still too large to be useful for my application. I need pretty good convergence on the order of 1e-3 or lower.

Verifying solution

I can verify my solution isn't actually a solution of the PDE's by using the method of manufactured solutions as detailed here.

I choose two arbitrary functions which agree with my boundary conditions: $E_{arb} = \frac{e^{-r^2}}{r^2}$ and $u_{arb} = 1 + \frac{1}{r}$.

I plug these expressions into the PDE's and since they dont actually solve the PDE's, there are residual terms left over, called source terms. I label these EQ1source and EQ2source. I subtract these from the original PDE's, EQ1 and EQ2, so that these arbitrary functions become actual solutions of these modified PDE's

arbEsol = Exp[-r^2]/r^2;
arbusol = 1 + 1/r;
EQ1source = 
  Activate[EQ1] /. {E1 -> Function[{r, \[Theta], \[Phi]}, arbEsol], 
     u -> Function[{r, \[Theta], \[Phi]}, arbusol]} // FullSimplify;
EQ2source =  
  Activate[EQ2] /. {E1 -> Function[{r, \[Theta], \[Phi]}, arbEsol], 
     u -> Function[{r, \[Theta], \[Phi]}, arbusol]} // Simplify;
newEQ1 = EQ1 - EQ1source;
newEQ2 = EQ2 - EQ2source;

I then use the same numerical procedure as above to numerically solve these PDE's

Block[{A = 1, r0 = 1, \[Mu] = 1, \[Lambda] = 0, M = 1, a = 0.},
 E1seed[r_?NumericQ, \[Theta]_?NumericQ, \[Phi]_?NumericQ] := 
  Evaluate[arbEsol];
 useed[r_?NumericQ, \[Theta]_?NumericQ, \[Phi]_?NumericQ] := 
  Evaluate[arbusol];
 
 {state} = NDSolve[{newEQ1 == EIC, newEQ2 == uIC},
    {E1, u},
   {r, \[Epsilon], Rmax}, {\[Theta], 0, \[Pi]}, {\[Phi], 0, 2 \[Pi]},
   Method -> {"FiniteElement", "BoundaryTolerance" -> None, 
     "MeshOptions" -> {"MeshOrder" -> 2}},
   InitialSeeding -> {E1[r, \[Theta], \[Phi]] == 
      E1seed[r, \[Theta], \[Phi]], 
     u[r, \[Theta], \[Phi]] == useed[r, \[Theta], \[Phi]]}
   ];
 Ediff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[newEQ1] /. state]];
 Udiff = 
  Function[{r, \[Theta], \[Phi]}, Evaluate[Activate[newEQ2] /. state]];
 E1sol = E1 /. state;
 Usol = u /. state;
 ]

Comparing the numerical solution to the known analytic solution gives us a measure of the error of the numerical scheme. I plot the numerical solution and analytic solution for u:

LogLinearPlot[{arbusol, state[[2, 2]][r, 0.5, 0.5]}, {r, 0.3, 10}, 
 PlotRange -> All]

Clearly, the solutions dont match. Meaning something is really happening with my desired solutions above. I'm not sure where though...