FEM doesn't output exactly the same solution for exactly the same code?

Question

Consider this toy example:

sol := 
 NDSolveValue[{Laplacian[u[x, y], {x, y}] + 1 == 0, 
               DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}]

sol == sol
(* True *)

sol - sol
(* - InterpolatingFunction[…] + InterpolatingFunction[…] *)

sol === sol
(* False *)

Table[sol[0.6, 0.6] - sol[0.6, 0.6], {10}]
(* {3.05311*10^-16, -1.249*10^-16, -9.71445*10^-17, 
   -4.85723*10^-16, 9.71445*10^-17, 2.77556*10^-17, 
   -5.55112*10^-17, 0., 1.80411*10^-16, 2.498*10^-16}  *)

Tested in v12.3 and v13.2.

What's happening here?

Is it possible to force NDSolveValue to produce exactly the same output for every execution?

Answer 1

You are right to expected that a deterministic algorithm should produce the exact same result every time it is run. And solving the linear equations of this elliptic PDE should be deterministic. Or not?

In fact, there are several algorithms that use some randomness. Mesh generation is such a thing as user21 mentioned.

I'd like to point to another source of randomness: the linear solver that finally solves the linear equations. Or more precisely the algorithm that permutes rows and columns of the sparse system matrix. Here the aim is to find a permutation such that Gaussian elimination introduces as little fill-in as possible (and some other subtle side conditions, e.g. that the induced elimination tree is well balanced and not too deep ). This is a NP hard problem. So one has to settle with heuristics. I don't know which exact library is employed here. But I would not be surprised if it employs some randomization. That would imply that from run to run of the solver the Gauss algorithm eliminates the rows of the system matrix in different orders -- which would not matter in exact arithmetic but introduces some rounding errors in floating point arithmetic.

Answer 2

Inspired by Henrik's answer, I think I find the root of problem: "Pardiso" method of LinearSolve is the culprit. According to the tutorial NDSolve Options for Finite Elements:

The finite element method uses the efficient direct PARDISO solver as the default linear solver.

Let's test Pardiso:

n = 1000; m = 
 SparseArray[{{i_, i_} -> -2., {i_, j_} /; Abs[i - j] == 1 -> 1.}, {n, n}];
b = Table[0., {n}]; b[[99]] = 1.;

sollinear := LinearSolve[m, b, Method -> "Pardiso"]

sollinear - sollinear // Abs // Max
(* 2.13163*10^-12 *)

It's worth noting that n should be large enough. For small system e.g. n = 100, the randomness becomes 0.. This explains why Alex's method works. (MaxCellMeasure -> 10^-2 actually results in a mesh coarser than default. )

If we turn to following methods for LinearSolve, deterministic solution will be produced (in reasonable time):

solmethod[method_] := 
 NDSolveValue[{Laplacian[u[x, y], {x, y}] + 1 == 0, 
   DirichletCondition[u[x, y] == 0, True]}, u, {x, 0, 1}, {y, 0, 1}, 
  Method -> {"PDEDiscretization" -> {"FiniteElement", 
      "PDESolveOptions" -> {"LinearSolver" -> {Automatic, Method -> method}}}}]

Table[
   solmethod[#] === solmethod[#] & /@ {"Multifrontal", "Krylov", 
     "IterativeRefinement", "Banded"}, {10}] // Flatten // Union
(* {True} *)

Still, this isn't the end: what if the Pardiso method is necessary for certain equation solving? Is it possible to adjust option of Pardiso method to remove its randomness?

Answer 3

Contra example

Needs["NDSolve`FEM`"]
mesh = ToElementMesh[FullRegion[2], {{0, 1}, {0, 1}}, 
  MaxCellMeasure -> 10^-2]

(*Out[]= ElementMesh[{{0., 1.}, {0., 1.}}, {QuadElement["<" 100 ">"]}]*)

 Do[
 sol[i] = 
   NDSolveValue[{Laplacian[u[x, y], {x, y}] + 1 == 0, 
     DirichletCondition[u[x, y] == 0, True]}, u, 
    Element[{x, y}, mesh]];, {i, 10}]

Table[sol[1][6/10, 6/10] - sol[i][6/10, 6/10], {i, 10}]

{0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`, 0.`}