NIntegrate doesn't work on regional integration

Question

I use FEM to solve a PDE in a spherical region and I want to integrate the result within half of the spherical region. I wonder why NIntegrate doesn't work here. The integrant should not have singularity cause the solution is finite everywhere.

Clear[u];
reg = Ball[{0, 0, 0}, 1];
sol = NDSolveValue[{D[u[t, x, y, z], t] - \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\(u[t, x, y, 
       z]\)\) == NeumannValue[0, True], 
   u[0, x, y, z] == 0.5 Tanh[30 z] + 0.5}, 
  u, {t, 0, 1}, {x, y, z} \[Element] reg]
reg1 = ImplicitRegion[x^2 + y^2 + z^2 <= 1 && z > 0, {x, y, z}];
data = Table[
  NIntegrate[sol[t, x, y, z], {x, y, z} \[Element] reg1, 
   AccuracyGoal -> 6, PrecisionGoal -> 6], {t, .05, 1, .05}]
ListPlot[data]

It gives

If I use Method -> "MonteCarlo" it works but with low precision even when using a lot of sampling points. Is there a way to accurately evaluate this integration?

Answer 1

The problem is not with NIntegrate -- the sol function gives Indeterminate for some points of the region reg1:

SeedRandom[343];
lsRPoints = RandomPoint[reg1, 200];
Select[Association[
  Map[# -> sol[0.1, Sequence @@ #] &, lsRPoints]], ! NumberQ[#] &]

During evaluation of In[51]:= InterpolatingFunction::femdmval: Input value {-0.928406,-0.347651,0.0901622} lies outside the range of data in the interpolating function.

During evaluation of In[51]:= InterpolatingFunction::femdmval: Input value {-0.676274,-0.360619,0.637002} lies outside the range of data in the interpolating function.

During evaluation of In[51]:= InterpolatingFunction::femdmval: Input value {-0.165827,-0.415793,0.890016} lies outside the range of data in the interpolating function.

During evaluation of In[51]:= General::stop: Further output of InterpolatingFunction::femdmval will be suppressed during this calculation.

(* <|{-0.928406, -0.347651, 0.0901622} -> 
  Indeterminate, {-0.676274, -0.360619, 0.637002} -> 
  Indeterminate, {-0.165827, -0.415793, 0.890016} -> Indeterminate|> *)

Defining a function that, say, replaces Indeterminate with 0 would help NIntegrate to give numerical results:

Clear[sol2];
sol2[t_?NumericQ, x_?NumericQ, y_?NumericQ, z_?NumericQ] :=
  Block[{r},
   r = sol[t, x, y, z];
   If[NumberQ[r], r, 0]
  ];

AbsoluteTiming[
 data = Table[
   NIntegrate[sol2[t, x, y, z], {x, y, z} \[Element] reg1, 
    AccuracyGoal -> 6, PrecisionGoal -> 6, 
    MaxRecursion -> 4], {t, .05, 1, .05}]
 ]
(*Lots of warning messages*)

(*{280.363, {1.68084, 1.53627, 1.43274, 1.35325, 1.29053, 1.24077, 
  1.20077, 1.16844, 1.14222, 1.12142, 1.10546, 1.09276, 1.08204, 
  1.07332, 1.06635, 1.06092, 1.05636, 1.05272, 1.04979, 1.04745}}*)

Remark: Note that I am using MaxRecursion->4 since I did not want to wait for too long...

Remark: Examine the messages -- they indicate some convergence problems.

Here we plot the result:

ListPlot[data, PlotRange -> All]

Answer 2

There is two ways you can go about this. As Anton noted, NDSolve will return Indeterimate for solution function evaluations that are outside of the region. This is the correct behavior for the finite element method as generally speaking no information is available beyond the boundary condition. That behavior, however, can be changed. This is explained in the ExtrapolationHandler section of the NDSolve Finite Element Options tutorial.

Clear[u];
reg = Ball[{0, 0, 0}, 1];
sol = NDSolveValue[{D[u[t, x, y, z], t] - \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\(u[t, x, y, 
       z]\)\) == NeumannValue[0, True], 
      u[0, x, y, z] == 0.5 Tanh[30 z] + 0.5}, 
    u, {t, 0, 1}, {x, y, z} \[Element] reg, 
  "ExtrapolationHandler" -> {Automatic,
              "WarningMessage" -> False}]

Now, the function evaluation uses extrapolation and gives no message:

sol[1, 2, 2, 2]
0.46176609059001855`

The second approach would be to use the same mesh both for NDSolve and NIntegrate. First we create a boundary mesh that has an interface at z==0:

Needs["OpenCascadeLink`"]
s1 = OpenCascadeShape[Ball[]];
s2 = OpenCascadeShape[
   Polygon[{{-2, -2, 0}, {2, -2, 0}, {2, 2, 0}, {-2, 2, 0}}]];
surface = OpenCascadeShapeFaces[s1];
innerDisk = OpenCascadeShapeIntersection[s1, s2];
boundary = OpenCascadeShapeUnion[Join[surface, {innerDisk}]];
(bmesh = OpenCascadeShapeSurfaceMeshToBoundaryMesh[
    boundary])["Wireframe"]

Next, we generate the mesh:

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[bmesh, 
   "RegionMarker" -> {{{0, 0, -1/2}, 1}, {{0, 0, 1/2}, 2}}];

Let's visualize:

parts = Map[
  mesh["Wireframe"[ElementMarker == #[[1]], 
     "MeshElement" -> "MeshElements", 
     "ElementMeshDirective" -> 
      Directive[EdgeForm[], FaceForm[#[[2]]]]]] &, {{1, Gray}, {2, 
    Orange}}]

Rasterize[Show[parts, PlotRange -> {All, {-0.02, 2}, All}]]

Solve:

sol = NDSolveValue[{D[u[t, x, y, z], t] - \!\(
\*SubsuperscriptBox[\(\[Del]\), \({x, y, z}\), \(2\)]\(u[t, x, y, 
       z]\)\) == NeumannValue[0, True], 
      u[0, x, y, z] == 0.5 Tanh[30 z] + 0.5}, 
    u, {t, 0, 1}, {x, y, z} \[Element] mesh]

data = Table[
  NIntegrate[If[z > 0, sol[t, x, y, z], 0], {x, y, z} \[Element] mesh, 
      AccuracyGoal -> 6, PrecisionGoal -> 6], {t, .05, 1, .05}]
ListPlot[data]