Trying to solve the following PDE with BC `T==1` on a spherical cap of a unit sphere and `T==0` at infinity (approximated as `r==(x^2 + y^2 + z^2)^0.5==40^0.5`) and the flux over the remaining surfaces taken to be zero (only half domains has been specified due to symmetry reasons):

 p = 0.2;
 Pe = 20;
 << NDSolve`FEM`
 boundaries = {-x^2 - y^2 - z^2 + 1, x^2 + y^2 + z^2 - 40, -y, 
 z - p + 1};
 \[CapitalOmega] = 
 ToElementMesh[
 ImplicitRegion[And @@ (# <= 0 & /@ boundaries), {x, y, z}]];
 Show[\[CapitalOmega][
 "Wireframe"["MeshElement" -> "MeshElements", Boxed -> True]], 
 Axes -> True, AxesLabel -> {"x", "y", "z"}, 
 PlotRange -> {{-7, 7}, {-0.5, 7}, {-7, 1}}]
 Show[\[CapitalOmega]["Wireframe"], Axes -> True, 
 AxesLabel -> {"x", "y", "z"}, 
 PlotRange -> {{-7, 7}, {-0.5, 7}, {-7, 1}}]

[![dominio][1]][1]

The free-surface is located at `z==-0.8` in the chosen reference system.
The differential equation is:

 sol = NDSolveValue[{D[T[x, y, z], x] == 
 1/Pe Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
 T[x, y, z] == 1., boundaries[[1]] == 0.],
 DirichletCondition[T[x, y, z] == 0., boundaries[[2]] == 0.]}}, 
 T, {x, y, z} \[Element] \[CapitalOmega]]
I noticed a non-smooth behavior of the solution `sol`, as confirmed by the density plots:

 z1 = -0.8;
 DensityPlot[sol[x, y, z1], {x, -4, 4}, {y, 0, 2}, PlotRange -> All, 
 PlotPoints -> 100, AspectRatio -> 1/2]
 DensityPlot[sol[x, 0, z], {x, -4, 4}, {z, -0.8, -2}, PlotRange -> All,
 PlotPoints -> 100, AspectRatio -> 1/2]
[![fig1][2]][2]
[![fig2][3]][3]

and by some plots of the solution

 Plot[sol[x, 0, -0.8], {x, 0.6, 6.1}, Frame -> True, 
 PlotRange -> {{-0.1, 7}, {-0.1, 1.2}}]
 Plot[sol[x, 0, -0.8], {x, -6.1, -0.6}, Frame -> True, 
 PlotRange -> {{-7, 0.1}, {-0.1, 1.2}}]

[![fig3][4]][4]
[![fig4][5]][5]


The derivatives of the solution show an even worse behavior. I report here as example just the derivative with respect to `x`:

 Dr[x_, y_, z_] = D[sol[x, y, z], x]
 Plot[Dr[x, 0, -0.8], {x, 0.6, 8}, Frame -> True]
 Plot[Dr[x, 0, -0.8], {x, -8, -0.6}, Frame -> True]
[![fig5][6]][6]
[![fig6][7]][7]


I am interested in the derivative because I am trying to calculate the total flux of the gradient of `sol` through the two portion of the spherical cap, `SC1` and `SC2`, which border the domain. As known, this is proportional to the heat flow through `SC1` and `SC2`. The definition of `SC1` and `SC2` is:

 SC1 = ImplicitRegion[
 x^2 + y^2 + z^2 == 1 && z <= p - 1 && y >= 0 && x >= 0, {x, y, z}];
 SC2 = ImplicitRegion[
 x^2 + y^2 + z^2 == 1 && z <= p - 1 && y >= 0 && x <= 0, {x, y, z}];
 Show[DiscretizeRegion[SC1, {{-5, 5}, {-5, 5}, {-3, 3}}, 
 MaxCellMeasure -> 0.0001], Axes -> True, 
 AxesLabel -> {"x", "y", "z"}, 
 PlotRange -> {{-1, 1}, {-0.2, 1}, {-0.8, -1}}]
 Show[DiscretizeRegion[SC2, {{-5, 5}, {-5, 5}, {-3, 3}}, 
 MaxCellMeasure -> 0.0001], Axes -> True, 
 AxesLabel -> {"x", "y", "z"}, 
 PlotRange -> {{-1, 1}, {-0.2, 1}, {-0.8, -1}}]

and the heat flow through `SC1` and `SC2` should be:

 NIntegrate[#, {x, y, z} \[Element] SC1] & /@ (Grad[
 sol[x, y, z], {x, y, z}].{x, y, z})
 NIntegrate[#, {x, y, z} \[Element] SC2] & /@ (Grad[
 sol[x, y, z], {x, y, z}].{x, y, z})

Unfortunately, I get a lot of errors from the last calculation, I don't know if this happens because of the non-smooth behavior of the solution and of the derivative or for other mistakes that I am doing. Thank you in advance.


 [1]: https://i.sstatic.net/VrfMr.gif
 [2]: https://i.sstatic.net/S4c63.gif
 [3]: https://i.sstatic.net/5P7Dm.gif
 [4]: https://i.sstatic.net/NRs8c.gif
 [5]: https://i.sstatic.net/Sw4XW.gif
 [6]: https://i.sstatic.net/I0bTo.gif
 [7]: https://i.sstatic.net/6EKBK.gif