# Postprocessing of FEM: extremely slow speed of NIntegrate

I've gone through the documentation for NDSolve FEM. There are quite a lot tips and techniques for accelerating NDSolve, and only quite few tips for post processing. Somehow I found out that the numerical integral of the solution from NDSolve takes a great amount of time, even much longer than NDSolve itself.

Take this fluid/solid domain as an example.

Needs["NDSolveFEM"]
{L, b, h, ri} = {0.5, 9.3, 1.25, 0.01}
hi = 2/3 h
m2D = ToElementMesh[
ImplicitRegion[-L/2 <= x <= L/2 && 0 <= y <= 1.25 && ! ((x^2 + (y - hi)^2 < ri^2)), {x, y}],
"RegionHoles" -> None, "MeshOrder" -> 1, MaxCellMeasure -> Infinity]
coords = m2D["Coordinates"];
c1 = Join[coords, ConstantArray[{0.}, {Length[coords]}], 2];
c2 = Join[coords, ConstantArray[{b}, {Length[coords]}], 2];
ele = m2D["MeshElements"];
max = Max[ElementIncidents[ele]]
tohigherIndex[type_[inci_, marker___], max_] :=
type[inci + max, marker]
sides = QuadElement[Map[Join[#1, Reverse[#1 + max]] &, #]] & /@
ElementIncidents[
Flatten[MeshElementSplitByMarker[m2D["BoundaryElements"]]]];
bm3D = ToBoundaryMesh["Coordinates" -> Join[c1, c2],
"BoundaryElements" ->
Join[ele, tohigherIndex[#, max] & /@ ele, sides]]
m3D = ToElementMesh[bm3D, MaxCellMeasure -> 1,
"MeshElementType" -> TetrahedronElement,
"MeshOrder" -> 1]


Here are the conservation equations and the boundary conditions.

w = 0.015
λs = 1.8;
ρs = 2104;
cps = 776;
te = 10;
ti = 20;

λf = 0.6;
ρf = 1000;
cpf = 4180;

eq = -Piecewise[{{w, x^2 + (y - hi)^2 < ri^2}}, 0] ρf cpf D[
t[x, y, z, tau], z] -
Piecewise[{{ρf cpf , x^2 + (y - hi)^2 < ri^2}}, ρs cps]/
3600/24 D[t[x, y, z, tau], tau] +
Inactivate[
Div[Piecewise[{{{{λf, 0, 0}, {0, λf, 0}, {0,
0, λf}}, x^2 + (y - hi)^2 < ri^2}}, {{λs,
0, 0}, {0, λs, 0}, {0, 0, λs}}].Grad[
t[x, y, z, tau], {x, y, z}], {x, y, z}], Div | Grad];
bc1 = NeumannValue[0, x == -L/2] + NeumannValue[0, x == L/2] +
NeumannValue[0.1 (t[x, y, z, tau] - te), z == 0 && x^2 + (y - hi)^2 >= ri^2] +
NeumannValue[0.1 (t[x, y, z, tau] - te), z == b && x^2 + (y - hi)^2 >= ri^2] +
NeumannValue[0, z == b && x^2 + (y - hi)^2 < ri^2] +
NeumannValue[0.1 (t[x, y, z, tau] - ti), y == 1.25] +
NeumannValue[0.1 (t[x, y, z, tau] - te), y == 0];
bc2 = DirichletCondition[t[x, y, z, tau] == 40, z == 0 && x^2 + (y - hi)^2 <= ri^2];


Solving this equation takes only 26 seconds on my laptop.

end = 45;
(tDyn = NDSolveValue[{
eq == bc1,
bc2,
t[x, y, z, 0] ==
Piecewise[{{40, z == 0 && x^2 + (y - hi)^2 <= ri^2}}, 15]}, t, {x, y, z} ∈ m3D, {tau, 0, end},
"ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False},
Method -> {"FiniteElement", "LinearSolveMethod" -> {"Pardiso"}},
EvaluationMonitor :> (currentTime = tau; percent = tau/end),
AccuracyGoal -> 4, PrecisionGoal -> 4];) // AbsoluteTiming
(*{26.1435, Null}*)


However, the integration of the solution takes more than a minute.

AbsoluteTiming[
NIntegrate[
tDyn[x, y, z, 0], {x, y, z} ∈ tDyn["ElementMesh"],
Method -> {Automatic, "SymbolicProcessing" -> 0},
AccuracyGoal -> 4, PrecisionGoal -> 4]/(L b h);]
(*{62.3974, Null}*)


I could imagine, it could last forever, if I want to plot the integration over the time variable tau. Is there any possibility to make the NIntegrate faster?

EDIT

Accelerating NDSolve by giving an initial condition from a steady state calculation.

(* initial condition from steady-state *)
(t0 = NDSolveValue[{
Inactivate[
Div[Piecewise[{{{{λf, 0, 0}, {0, λf, 0}, {0, 0, λf}},
x^2 + (y - hi)^2 <= ri^2}},
{{λs, 0, 0}, {0, λs, 0}, {0, 0, λs}}].
Grad[t[x, y, z], {x, y, z}], {x, y, z}], Div | Grad] ==

NeumannValue[0, x == -L/2]
+ NeumannValue[0, x == L/2]
+ NeumannValue[0.1 (t[x, y, z] - te), z == 0]
+ NeumannValue[0.1 (t[x, y, z] - te), z == b]
+ NeumannValue[0.1 (t[x, y, z] - ti), y == 1.25]
+ NeumannValue[0.1 (t[x, y, z] - te), y == 0]},
t, {x, y, z} \[Element] m3D,
"ExtrapolationHandler" -> {Automatic,
"WarningMessage" -> False}];) // AbsoluteTiming
(*{0.792426, Null}*)

end = 45;
(Monitor[tDyn =
NDSolveValue[{eq == bc1, bc2, t[x, y, z, 0] == Piecewise[{
{40, z == 0 && x^2 + (y - hi)^2 <= ri^2}}, t0[x, y, z]]},
t, {x, y, z} \[Element] m3D, {tau, 0, end},
"ExtrapolationHandler" -> {Automatic, "WarningMessage" -> False},
Method -> {"FiniteElement", "LinearSolveMethod" -> {"Pardiso"}},
EvaluationMonitor :> (currentTime = tau; percent = tau/end),
AccuracyGoal -> 4, PrecisionGoal -> 4];,
Row[{ProgressIndicator[(currentTime - 0)/end], "   ",
ToString[Floor[100 percent]] <> "%", "    ",
"tau = " <> ToString[currentTime]}]]) // AbsoluteTiming
(*{18.5054, Null}*)

AbsoluteTiming[
NIntegrate[
tDyn[x, y, z, 0], {x, y, z} \[Element] tDyn["ElementMesh"],
Method -> {Automatic, "SymbolicProcessing" -> 0},
AccuracyGoal -> 4, PrecisionGoal -> 4]/(L b h);]
(*{48.3468, Null}*)

Plot[NIntegrate[tDyn[x, y, z, tau], {x, y, z} \[Element] tDyn["ElementMesh"]]/(L b h), {tau, 0, end}]
(* this plot takes forever *)

• L is undefined. – bbgodfrey Jul 26 '18 at 1:24
• Oh sorry, now it's been given – 407PZ Jul 26 '18 at 6:15
• I cannot reproduce your solution tDyn, MMA ( 11.0.1 Windows ) takes forever. – Ulrich Neumann Jul 26 '18 at 7:07
• @UlrichNeumann thanks for pointing that out. I've given a new initial condition from steady-state calculation. It largely accelerates the calculation. It's also very interesting to see how NDSolve behaves under different initial conditions. – 407PZ Jul 26 '18 at 8:11
• @ 407Peezy: If I understand your question , you are looking vor the timedependig integral over region m3D? Perhaps you can put this integration inside the NDSolve-command as a discrete variable to run your code faster? – Ulrich Neumann Jul 26 '18 at 10:08