NDSolve becomes very slow when the integration domain is 4 dimensional

Question

I am trying to solve a PDE problem on a 4d manifold (3 space + 1 time). However, NDSolve becomes very slow when I jump from 3 to 4 dimensions. Here I want to integrate a simple ndeps dimensional diffusion equation on a nindeps+1 dimensional manifold.

TimeNDSolve[nindeps_, ndeps_] := Module[
   {t, indeps, man, deps, [CapitalOmega], BCs, ICs, eqs},
   man = Table[Symbol["$x" <> ToString@i], {i, nindeps}];
   indeps = {t}~Join~man;
   deps = Table[Symbol["$y" <> ToString@i], {i, nindeps}];
   [CapitalOmega] = Cuboid@Table[0, {i, 1, Length@man}];
   BCs = {DirichletCondition[deps[[1]] @@ indeps == Sin[3 t], 
      indeps[[2]] == 0]};
   ICs = Flatten@Table[{q @@ (indeps /. t -> 0) == 0}, {q, deps}];
   eqs = Table[
     D[q @@ indeps, {t, 1}] == Laplacian[q @@ indeps, man], {q, 
      deps}]; AbsoluteTiming[
    nsol = NDSolve[{eqs, BCs, ICs} // Flatten, deps, 
      man [Element] [CapitalOmega], {t, 0, 1}]]];

For nindeps =1 the program runs basically constant as ndeps is increased (Try Table[TimeNDSolve[i, j][[1]], {i, 1, 1}, {j, 1, 100}]). Timing says this takes $<1$ second per function call. If we increase to nindeps =2 there is no significant slowdown.

However, if nindeps is increased to 3, suddenly the function takes ~170 seconds to run. I first thought that increasing dimensions would exponentially increase the number of simplices that any finite element solver would use to cover the base manifold, but this is more than exponential growth! It is absolutely ridiculous! I want to simulate a much more complicated PDE in this space which is currently unfeasible.

What is going on? Is there any way to speed things up? Thanks!

differential equations finite element method numerical integration

user21 · Accepted Answer

I think it's better to compare the number of degrees of freedom (DOF) for each dimension. (The DOF for a single PDE is roughly the number of coordinates in the mesh.) For that I have changed your code a bit:

Needs["NDSolve`FEM`"]
(*ord=1;*)
ord = 2;

Ω = Cuboid@Table[0, {i, 1, 1}];
m[1] = ToElementMesh[Ω, "MeshOrder" -> ord(*,
   MaxCellMeasure[Rule]0.00036*), MaxCellMeasure -> 0.0002];
Print[Length[m[1]["Coordinates"]]];

10001


Ω = Cuboid@Table[0, {i, 1, 2}];
m[2] = ToElementMesh[Ω, "MeshOrder" -> ord, 
   MaxCellMeasure -> 0.0003];
Print[Length[m[2]["Coordinates"]]];

10732

Ω = Cuboid@Table[0, {i, 1, 3}];
m[3] = ToElementMesh[Ω, "MeshOrder" -> ord, 
   "MeshElementType" -> 
    TetrahedronElement,(*MaxCellMeasure[Rule]0.00022,*)
   MaxCellMeasure -> 0.00041];
Print[Length[m[3]["Coordinates"]]];

10849

TimeNDSolve[nindeps_, ndeps_] := 
  Module[{t, indeps, man, deps, Ω, BCs, ICs, eqs}, 
   man = Table[Symbol["$x" <> ToString@i], {i, nindeps}];
   indeps = {t}~Join~man;
   deps = Table[Symbol["$y" <> ToString@i], {i, nindeps}];
   Ω = Cuboid@Table[0, {i, 1, Length@man}];
   BCs = {DirichletCondition[deps[[1]] @@ indeps == Sin[3 t], 
      indeps[[2]] == 0]};
   ICs = Flatten@Table[{q @@ (indeps /. t -> 0) == 0}, {q, deps}];
   eqs = Table[
     D[q @@ indeps, {t, 1}] == Laplacian[q @@ indeps, man], {q, deps}];
   AbsoluteTiming[
    Monitor[nsol = 
      NDSolve[{eqs, BCs, ICs} // Flatten, deps, 
       man ∈ m[nindeps], {t, 0, 1}, 
       EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])],
     monitor]]];

These are then the timings I see:

Table[TimeNDSolve[i, 0][[1]], {i, 1, 3}, {j, 1, 1}]

{{0.342661}, {1.9421}, {8.83981}}

As for speeding things up in 3D, you might have success with reducing the accuracy/precision goal; something reducing the BDF order is helpful, though probably not in this case given here.

NDSolve becomes very slow when the integration domain is 4 dimensional

One Answer

Add your own answers!

Ask a Question