3
$\begingroup$

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["NDSolve`FEM`"]
{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 *)
$\endgroup$
7
  • $\begingroup$ L is undefined. $\endgroup$ – bbgodfrey Jul 26 '18 at 1:24
  • $\begingroup$ Oh sorry, now it's been given $\endgroup$ – 407PZ Jul 26 '18 at 6:15
  • $\begingroup$ I cannot reproduce your solution tDyn, MMA ( 11.0.1 Windows ) takes forever. $\endgroup$ – Ulrich Neumann Jul 26 '18 at 7:07
  • $\begingroup$ @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. $\endgroup$ – 407PZ Jul 26 '18 at 8:11
  • 2
    $\begingroup$ @ 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? $\endgroup$ – Ulrich Neumann Jul 26 '18 at 10:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.