# Total flux of the gradient of the numerical solution of a PDE through a surface

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;
<< NDSolveFEM
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}}]

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]

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}}]

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]

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.

• Mesh is too bad. May be you should restrict solution to x^2+y^2+z^2<=4&&x>=0&& y >= 0 && z <= -1 + p, with mesh of 50323 Tetrahedron Element (MaxCellMeasure -> .0001), check solution and then try to expand region. Commented May 28, 2021 at 16:10

The basic problem appears to be a convective-diffusive heat transfer problem of X-directed fluid flow across a heated spherical cap tip. To study this type of problem, it probably is easier to construct a virtual cuboid wind tunnel. When simulating virtual wind tunnels, the upstream section is typically much shorter than the downstream week region, so there really is not any spherical symmetry.

The following workflow will show how to construct a virtual wind tunnel and add refinement zones so that gradients may be captured near the object of interest without blowing up the total model size.

# Virtual wind tunnel construction with refinement zones.

As discussed here, a MeshRefinementFunction will not necessarily refine the surface mesh in 3D. The suggested workaround was to use BoundaryDiscretizeRegion to obtain a finely discretized surface mesh on the input geometry before applying the MeshRefinementFunction.

The following workflow creates a refined region with a finely discretized tip and a wind tunnel domain less refined region:

<< NDSolveFEM
Pe = 15;
tip = BoundaryDiscretizeRegion[Ball[], MaxCellMeasure -> .00125,
Axes -> True, AxesLabel -> {"X", "Y", "Z"}];
refCuboid =
BoundaryDiscretizeRegion[Cuboid[{-1.5, 0, -1.5}, {3.5, 1.5, -0.8}],
MaxCellMeasure -> .01, Axes -> True, AxesLabel -> {"X", "Y", "Z"}];
reftip = RegionDifference[refCuboid, tip, Axes -> True,
AxesLabel -> {"X", "Y", "Z"}]
domCuboid =
BoundaryDiscretizeRegion[Cuboid[{-2.5, 0, -3.8}, {7.5, 4, -0.8}],
Axes -> True, AxesLabel -> {"X", "Y", "Z"}];
domref = RegionDifference[domCuboid, refCuboid, Axes -> True,
AxesLabel -> {"X", "Y", "Z"}]

## Create volume mesh with refinement zones

The following workflow creates a refinement zone with a high level of surface discretization at the tip. The total number of elements is about 134,000, which does not take too long to solve.

(* Create Mesh Refinement Function *)
mrf = With[{rmf = RegionMember[reftip]},
Function[{vertices, volume},
Block[{x, y, z}, {x, y, z} = Mean[vertices];
If[rmf[{x, y, z}], volume > (0.07/1.5)^3, volume > (0.3)^3]]]];
(* Create and Display Volumetric Mesh *)
mesh = ToElementMesh[RegionUnion[domref, reftip],
MeshQualityGoal -> "Maximal", "MeshElementConstraint" -> 40,
"MaxBoundaryCellMeasure" -> {"Length" -> .1},
MeshRefinementFunction -> mrf]
mesh["Wireframe"[
"MeshElement" -> "MeshElements",
"ElementMeshDirective" -> Directive[EdgeForm[Black]],
PlotRange -> {{-2.5, 7.5}, {0.2,
4.}, {-3.8, -0.7999999999999998}}]]
groups = mesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp
mesh["Wireframe"["MeshElementStyle" -> FaceForm /@ colors]]

# Solution

We and solve and now compare our solution on the refined mesh versus the OP mesh.

sol = NDSolveValue[{D[T[x, y, z], x] ==
1/Pe Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
T[x, y, z] == 1., ElementMarker == 3],
DirichletCondition[T[x, y, z] == 0., ElementMarker == 4]}},
T, {x, y, z} ∈ mesh];
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]

Visually, the solution appears to be much smoother than the images shown in the OP.

## Other plots

The other plots appear much smoother with the refined mesh.

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}}]
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]

# Integration of total flux

## Using NIntegrate

The integration strategy given in the OP does solve albeit quite slowly due to slow convergence issues. When I evaluate the OP code, I obtained:

NIntegrate[#, {x, y, z} ∈ SC1] & /@ (Grad[
sol[x, y, z], {x, y, z}] . {x, y, z})
NIntegrate[#, {x, y, z} ∈ SC2] & /@ (Grad[
sol[x, y, z], {x, y, z}] . {x, y, z})
(* -0.95414 *)
(* -2.96035 *)

## Summing up the discretized data

Using a variety of Mathematica functions, we can extract the normals, areas, and thermal gradients of each triangle on the mesh corresponding to the spherical cap. This should be enough information to estimate an integrated flux.

The following workflow shows how to extract the necessary information for the left and right-hand sides of the spherical cap.

(*Element info shortcuts*)
ebi = ElementIncidents[#["BoundaryElements"]][[1]] &;
ebm = ElementMarkers[#["BoundaryElements"]][[1]] &;
ebn = #["BoundaryNormals"][[1]] &;
ei = ElementIncidents[#["MeshElements"]][[1]] &;
em = ElementMarkers[#["MeshElements"]][[1]] &;
epi = ElementIncidents[#["PointElements"]][[1]] &;
epm = Flatten@ElementMarkers[#["PointElements"]] &;
(*extract boundary mesh from element mesh*)
bmesh = ToBoundaryMesh[mesh];
bcrd = bmesh["Coordinates"];
bi = ebi[bmesh];(*boundary element incidents*)
bm = ebm[bmesh];(*boundary element markers*)
bn = ebn[bmesh];(*boundary normals*)
(*find markers corresponding to the spherical cap*)
mrk3pos = Flatten@Position[bm, 3, 1];
(*generate necessary info to estimate surface integral*)
bn3 = bn[[mrk3pos]];
polys = Map[Polygon,
GetElementCoordinates[bcrd, #] & /@ bi[[mrk3pos]]];
area3 = Area /@ polys;
center3 = Map[Mean, GetElementCoordinates[bcrd, #] & /@ bi[[mrk3pos]]];
(*find positions of left and right side of spherical cap*)
posXids = Position[center3[[All, 1]], _?(# >= 0 &), 1] // Flatten;
negXids = Complement[Range[Length[mrk3pos]], posXids];
Show[{Graphics3D[{Red, polys[[posXids]]}],
Graphics3D[{Blue, polys[[negXids]]}]}]

There is a bit of jaggedness at the seam of the left and right-hand sides. Hopefully, the errors will average out. We could have put a seam at the interface with more elaborate model construction. Now, we can estimate the fluxes using the following:

gradT[x_, y_, z_] = {Derivative[1, 0, 0][sol][x, y, z],
Derivative[0, 1, 0][sol][x, y, z],
Derivative[0, 0, 1][sol][x, y, z]};
f = NDSolveFEMMapThreadDot[(gradT @@@ center3[[#]]), bn3[[#]]] .
area3[[#]] &;
f[posXids]
f[negXids]
f[posXids~Join~negXids]
(* -0.952311 *)
(* -2.96147 *)
(* -3.91378 *)

These results agree with the OP integral formulation quite well and are significantly faster.

# Comparison to another code

When possible, it is always good to compare the Mathematica FEM results with another FEM code. In this case, I will show that the Mathematica results compare favorably with the FEM code COMSOL. First, I will create a SliceContourPlot3D for comparison purposes.

surf = {{x^2 + y^2 + z^2 ==
1.001^2}, {"XStackedPlanes", {7.5}}, {"YStackedPlanes", {0}}, \
{"ZStackedPlanes", {-0.8}}, {"BackPlanes"}};
SliceContourPlot3D[sol[x, y, z], surf, {x, y, z} ∈ mesh,
Contours -> 11, PlotPoints -> 100, BoxRatios -> Automatic,
ColorFunction -> "ThermometerColors", PlotRange -> {-0.001, 1},
PlotLegends -> Automatic]

Here is the comparable COMSOL plot:

As you can see, the visualizations compare favorably.

We also can see that the integrations of gradT on the left and right hemispherical cap surfaces also agreed to within about 2%

Application Left Right
Mathematica -2.96147 -0.95414
COMSOL -3.0255 -0.952755
%Diff 2.16226 0.0465888

# Update in response to a comment about low Péclet numbers

The Péclet number is defined by:

$$Pe=\frac{Advective\ Transport\ Rate}{Diffusive\ Transport\ Rate}$$

So, as the Péclet number goes down the Advective component becomes less significant indicating that a "wind tunnel" is perhaps not the best model to study this problem. Although you may want to consider modifying the geometry, you can substantially mitigate the effects of low Péclet numbers by changing the default wall conditions to DirichletCondition's. Here is a logarithmic sweep from Pe numbers from 0.01 to 100 (note that this process is slow):

pfun = ParametricNDSolveValue[{D[T[x, y, z], x] ==
1/P Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
T[x, y, z] == 1., ElementMarker == 3],
DirichletCondition[T[x, y, z] == 0.,
Or @@ (ElementMarker == # & /@ {4, 6, 7})]}},
T, {x, y, z} ∈ mesh, {P}];
surf = {{x^2 + y^2 + z^2 ==
1.001^2}, {"XStackedPlanes", {7.5}}, {"YStackedPlanes", {0}}, \
{"ZStackedPlanes", {-0.8}}, {"BackPlanes"}};
frames = SliceContourPlot3D[pfun[#][x, y, z],
surf, {x, y, z} ∈ mesh, Contours -> 11,
PlotPoints -> 100, BoxRatios -> Automatic,
ColorFunction -> "ThermometerColors", PlotRange -> {-0.001, 1},
PlotLegends -> Automatic, PlotLabel -> N@#] & /@ (10^# & /@
Subdivide[-2, 2, 20]);
ListAnimate[frames]

• Is it possible that we can solve coupeld pde (e.g. FSI) using Mathametica Solver? Commented May 30, 2021 at 23:09
• @ABCDEMMM Of course, the subject of FSI is very broad and generally quite difficult. If you can assume that the deformations are so small that they do not affect the fluid flow (e.g., thermal stress due to transient transfer from a fluid), then you could do a "FSI" problem following this example. Significant deformation would require the ability for Mathematica to solve the mesh deformation equations on the fluid side, which I have yet to see an example. Consider making a feature request. Commented May 31, 2021 at 2:08
• thanks for your reply! Actually, I am looking for an example for Staggered Scheme for NDsolve coupled fields. or in other words, how to set the staggered solver in mathematica. Commented May 31, 2021 at 11:25
• @umby Thank you very much! I will see what I can do. Commented Jun 1, 2021 at 3:29
• @umby See update. Commented Jun 3, 2021 at 14:33

Here are a couple of different ideas for generating the meshes. This works in 12.3

Needs["NDSolveFEM"]
Pe = 15;
p = 2/10;
reg = RegionDifference[Cuboid[{-2.5, 0, -3.8}, {7.5, 4, -0.8}],
Ball[]];
bmesh = ToBoundaryMesh[reg, AccuracyGoal -> 3.5];

The key here is the AccuracyGoal. In version 12.3 the 3D default boundary mesh generator is OpenCascade. OpenCascade will not refine flat boundaries if the accuracy goal is increased.

bmesh["Wireframe"]

Visualize the boundary mesh and it's boundary markers:

groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
Show[
bmesh["Edgeframe"],
bmesh["Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementStyle" -> (Directive[EdgeForm[], FaceForm[#]] & /@
colors)]]
, Epilog -> Inset[LineLegend[colors, groups], Scaled[{0.85, 0.6}]]
]

The solution can then be found with

mesh = ToElementMesh[bmesh];
sol = NDSolveValue[{D[T[x, y, z], x] ==
1/Pe Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
T[x, y, z] == 1., ElementMarker == 7],
DirichletCondition[T[x, y, z] == 0., ElementMarker == 1]}},
T, {x, y, z} \[Element] mesh];

Another idea is to split the cap. This can be done with the OpenCascadeLink. The idea is to create a left and right part of the geometry, take the faces and merge them back together.

rleft = RegionDifference[Cuboid[{-2.5, 0, -3.8}, {0, 4, -0.8}],
Ball[]];
rright = RegionDifference[Cuboid[{0, 0, -3.8}, {7.5, 4, -0.8}],
Ball[]];

Note, that we have a split along the x=0 plane. Next, we generate a boundary mesh

"ShapeSurfaceMeshOptions" -> {(*"LinearDeflection"->0.00025*)
"AngularDeflection" -> 0.025
}])["Wireframe"]

You can tell OpenCascade to either use a linear deflection or an angular deflection (See documentation here)

Inspect the boundary mesh and it's markers

groups = bmesh["BoundaryElementMarkerUnion"];
temp = Most[Range[0, 1, 1/(Length[groups])]];
colors = ColorData["BrightBands"][#] & /@ temp;
Show[
bmesh["Edgeframe"],
bmesh["Wireframe"["MeshElement" -> "BoundaryElements",
"MeshElementStyle" -> (Directive[EdgeForm[], FaceForm[#]] & /@
colors)]]
, Epilog -> Inset[LineLegend[colors, groups], Scaled[{0.85, 0.6}]]
]

Inspect the detail of the cap:

bmesh["Wireframe"[ElementMarker == 7 || ElementMarker == 9,
"MeshElement" -> "BoundaryElements"]]

Generate the mesh:

mesh = ToElementMesh[bmesh];
mesh["Wireframe"]

Solve the PDE

sol = NDSolveValue[{D[T[x, y, z], x] ==
1/Pe Laplacian[T[x, y, z], {x, y, z}], {DirichletCondition[
T[x, y, z] == 1., ElementMarker == 7 || ElementMarker == 9],
DirichletCondition[T[x, y, z] == 0., ElementMarker == 1]}},
T, {x, y, z} \[Element] mesh];
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]

FEMNBoundaryIntegrate[
Grad[sol[x, y, z], {x, y, z}] . {x, y, z}, {x, y, z}, mesh,
ElementMarker == 7]

-2.98032

FEMNBoundaryIntegrate[
Grad[sol[x, y, z], {x, y, z}] . {x, y, z}, {x, y, z}, mesh,
ElementMarker == 9]

-0.946663

You can still use a MeshRefinementFunction to refine the mesh away from the cap.

As a side note, check out these new (12.3) OpenCascade tutorials: Book Shelf Bracket and Helical Bevel Gear

• Is FEMNBoundaryIntegrate a new addition? I like it. Also, is it still true that MeshRefinementFunction will only operate on the volume elements leaving a big jump in element size from the volume to the surface? I was trying to avoid that. I initially tried a combination of OpenCascade and FEMAddOns to union the bmesh's, but it failed. That is why I fell back to the  BoundaryDiscretizeRegion. Commented Jun 1, 2021 at 13:35
• @TimLaska, FEMNBoundaryIntegrate and FEMNIntegrate are not new, they are the back bone of NIntegrate for the cases when the method goes to FEM. I never really documented them because NIntegrate did a good job. However, they have one 'feature' that makes them valuable and that is the 4rth argument. Where you can specify a predicate. In principal one can have this predicate in the first argument like FEMNBoundaryIntegrate[If[pred, 1,0]*val,{x,y,..}, bmesh] but in some cases this gives the wrong result. So please be careful with that, probably better to use the 4 arg version. Commented Jun 2, 2021 at 4:46
• @TimLaska, concerning MeshRefinementFunction: this is a WL interface to a C language interface both Triangle and TetGen provide and both tools do not preform a mesh refinement on the boundary. So in theory what would need to happen is that the MeshRefinementFunction would need to be passed to the boundary mesh generator refine the boundary and then during the mesh generation that same function would then refine the full mesh. I'll think a bit about that, but please do not hold your breath for that - I am so consumed with structural mechancis right now that I can not foresee when I'll be Commented Jun 2, 2021 at 4:50
• @TimLaska, slightly tangential, I'd like to point you to another update, that I think, you might find useful: You can now use NeumannValue in NIntegrate. For example, if you want to compute the heat loss through a radiation boundary you can directly use the NeumannValue that specifies the radiation condition, like shown here. Hope this is useful. Commented Jun 2, 2021 at 5:12
• @TimLaska, I hope that before to long the anisotropic mesh generation you provided here will make it into the kernel, that should reduce the amount of programing needed to generate those kind of meshes; meshing is probably an area that one can always find new things to improve.... One other feature I'd really like to have is being able to specify 2D geometries in OpenCascade and then project them into 3D, for extrusion, lofting, etc. Commented Jun 3, 2021 at 3:41