Issue with boundary Integration of FEM numerical solution (interpolation function)

Question

I am having issues in integrating the numerical solution obtained with FEM method and NDSolve, over a boundary. Specifically, i am integrating a specific combination of the gradient over a boundary obtaining errors which are much larger than the specified tolerances. I would like to have a way to avoid the errors and to have a faster integration.

Let me give you the full story:

I am solving the equation $\boldsymbol{\nabla}^2 f(x,y) -(0.001)^2 f(x,y)=0$ in two dimensions and in cartesian coordinates.

My domain consists in a very big outer circle of radius $L=50000$ and two small circles $C_1$ and $C_2$ of radius $R=1$ whose $x$ and $y$ center's coordinates are $\left(- \frac{dist}{2},0 \right)$ and $\left(\frac{dist}{2},0 \right)$, respectively. $dist$ is a parameter that I vary to fix the circles center to center distance.

On the external circle i have no flux boundary condition (Neumann boundary condition). Whereas on the two circles i have Dirchelet boundary condition:

$f(x,y)=1 +(x-\frac{dist}{2})^2-y^2$ on $C_2$

$f(x,y)=1 +(x+\frac{dist}{2})^2-y^2$ on $C_1$

This is the code that i use to generate the mesh and solve the problem (I specify the mesh elements on the boundaries):

L = 50000;
dist = 5.0;
nptspart = 2000;
nptsout = 60;
ptspart = 
 Table[{Cos[2.0*Pi*i/nptspart] - dist/2, Sin[2.0 Pi*i/nptspart]}, {i, 0, nptspart - 1}];
elempart = 
 Table[If[i < Length[ptspart], {i, i + 1}, {Length[ptspart], 1}], {i, 1, Length[ptspart]}];
ptspart1 = 
 Table[{Cos[2.0*Pi*i/nptspart] + dist/2, Sin[2.0 Pi*i/nptspart]}, {i, 0, nptspart - 1}];
elempart1 = 
 Table[If[i < Length[ptspart1], {nptspart + i, 
    nptspart + i + 1}, {nptspart + Length[ptspart1], 
    nptspart + 1}], {i, 1, Length[ptspart1]}];
ptsout = Table[{L Cos[2.0*Pi*i/nptsout], L Sin[2.0 Pi*i/nptsout]}, {i,
    0, nptsout - 1}];
elemout = 
 Table[If[i < Length[ptsout], {2 nptspart + i, 
    2 nptspart + i + 1}, {2 nptspart + Length[ptsout], 
    2 nptspart + 1}], {i, 1, Length[ptsout]}];
bmesh = ToBoundaryMesh[
  "Coordinates" -> Join[ptspart, ptspart1, ptsout], 
  "BoundaryElements" -> {LineElement[
     Join[elempart, elempart1, elemout]]}, "MeshOrder" -> 2];
mesh1 = ToElementMesh[bmesh, "MeshOrder" -> 2, 
  "RegionHoles" -> {{-dist/2, 0}, {dist/2, 0}}, 
  MeshQualityGoal -> "Maximal"];
ufun1 = NDSolveValue[{\!\(
\*SubsuperscriptBox[\(∇\), \({x, y}\), \(2\)]\(f[x, 
       y]\)\) - (0.001)^2 f[x, y] == 
    NeumannValue[0.0, x^2 + y^2 == L^2], 
   DirichletCondition[
    f[x, y] == 
     1 + ((x - dist/2)^2 - y^2), (x - dist/2)^2 + y^2 - 1^2 == 0], 
   DirichletCondition[
    f[x, y] == 
     1 + ((x + dist/2)^2 - y^2), (x + dist/2)^2 + y^2 - 1^2 == 0]}, 
  f, {x, y} ∈ mesh1, 
  Method -> {"PDEDiscretization" -> "FiniteElement"}];

The parameter nptspart defines how many points i have on the circles boundaries $C_1$ and $C_2$, whereas nptsout defines the number of points on the bigger circle of radius L.

I have now to perform an integral on the boundary of each circle:

$int_1=\int_{C_1} \left( \boldsymbol{\nabla} f \boldsymbol{\nabla} f -\frac{1}{2} \boldsymbol{\nabla} f \cdot \boldsymbol{\nabla} f \boldsymbol{I} \right) \cdot \boldsymbol{n} \; dl$

$int_2=\int_{C_2} \left( \boldsymbol{\nabla} f \boldsymbol{\nabla} f -\frac{1}{2} \boldsymbol{\nabla} f \cdot \boldsymbol{\nabla} f \boldsymbol{I} \right) \cdot \boldsymbol{n} \; dl$

Where $\boldsymbol{I}$ is the identity matrix in two dimensions, $\boldsymbol{n}$ the outwardly directed normal to the circles boundary, $\nabla f \nabla f$ is the dyadic product (also called Kronecker product), and $dl$ is the differential element on the circle line.

In orer to perform the integration, on the circle $C_2$ for example, i use the following code:

NIntegrate[(KroneckerProduct[Grad[ufun1[x, y], {x, y}], 
     Grad[ufun1[x, y], {x, y}]] - (1/
       2) (Grad[ufun1[x, y], {x, y}].Grad[ufun1[x, y], {x, y}])*
     IdentityMatrix[2]).{x - dist/2, y}, {x, y} ∈ 
  Circle[{dist/2, 0}, 1], AccuracyGoal -> 6, PrecisionGoal -> 6, 
 MaxRecursion -> 100, Method -> "InterpolationPointsSubdivision"]

However not only i get a lot of warning messages like:

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 400 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.06632359368254395and 0.0005234490597411253 for the integral and error estimates.

but also the Integration is very slow, which is something I don't want as i will have to perform the integration, in my future research, for, say, 200 circles. therefore i need the integration to be fast on each boundary.

If i plot the results of the two integrals, for $dist=5$, as a function of the number of points for the circle boundaries (nptspart), I get something like this.

Where you can clearly see an error of the order of $\pm 0.001$ which is way bigger than the specified precision (6 digits) i specified in the integration.

Furthermore the error doesn't seem to reduce appreciably as the mesh on the circles gets refined.

POSSIBLE REASON: I think this has something to do with the interpolation function being evaluated in points which are not mesh points and are specified in the integration as:

  {x, y} ∈ Circle[{dist/2, 0}, 1]

Indeed if i let mathematica do a plot over the boundary $C_1$ (where the value of the function is specified by the Dirchelet boundary condition) i get weird results when i reduce the number of points on the boundary.

If i perform a parametric plot of the function $f(x(\theta),y(\theta))$ on $C_1$ (parametrised with the angle $\theta$ as shown in figure) for two different meshes

i get the following weird result:

In the above example i used directly the interpolation function obtained as result from the solution of NDSolve. An it is clear that is doing weird things on the boundary. So i guess that the gradient would be even worse.

So The question is. Is there a way to obtain the values of the function, and the gradient of the function, in the mesh nodes? Or if not, is it possible to perform an integration with the interpolation function that gives accurate results and is performed faster than it is at the moment? Do you have any tips?

Answer 1

More an extended comment than an answer but may help to illustrate my last comment.

If you plot the solution normalised by the boundary condition, you can see that making the mesh finer creates a high frequency noise. If you differentiate this, say with respect to the angle, will have very large values on the nodes of the elements. This is for nptspart=2000

Plot[ufun1[dist/2 + Cos[θ], 
   Sin[θ]]/(1 + (Cos[θ])^2 - 
    Sin[θ]^2), {θ, -.01, .02 π}]

I don't know the internals of how the derivative of an interpolated function is computed but I guess this noise is what makes the integral so time consuming.

As a workaround I found that interpolating in fewer points gives you a smoother solution. This is the code I use, which is not directly applicable to your problem but will give you an idea:

First I compute the gradient

gradC[r_, z_] := Evaluate[D[concentration[r, z], {{r, z}}]];

then the tangential component of the gradient on each particle (mine are in +- z0)

vSlip[r_, z_] := 
  Module[{ϕ = ArcTan[r, Abs[z] - z0]}, 
   gradC[r, 
     z] - {r, z - Sign[z] z0} ({r, z - Sign[z] z0}.gradC[r, z])];

And with that create an interpolated function fixing the values so it goes to zero on the axis (this was were my problem was)

fInterpol1 = 
  Interpolation[
   Table[{ϕ, 
     If[Abs[ϕ] >= π/2, 0, 
      vSlip[Cos[ϕ], z0 + Sin[ϕ]][[
       1]]]}, {ϕ, -π, π, π/5}], Method -> "Spline"];
fInterpol2 = 
  Interpolation[
   Table[{ϕ, 
     If[Abs[ϕ] >= π/2, 0, 
      vSlip[Cos[ϕ], z0 + Sin[ϕ]][[
       2]]]}, {ϕ, -π, π, π/5}], Method -> "Spline"];

with this, I can use a more smooth function that's fast to integrate.

vSlip2[r_, z_] := 
  Module[{ϕ = ArcTan[r, Abs[z] - z0]}, {fInterpol1[ϕ], 
    fInterpol2[ϕ]}];

In your case you should interpolate in 2D as I think yours is the full gradient that is used.

However, we were interested on how the integral of vslip changed with the distance between the particles and due to the errors we could not find anything conclusive, thus we are using BEM now.