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

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?

• It is not so clear what do you try to test? Nor any numerical parameters are fixed. Can you compare sol on the same mesh? Jan 6 at 6:10
• @alex I'm testing the equivalence of different PDE systems. (To be more specific, I'm exploring this question: mathematica.stackexchange.com/q/278215/1871 ) I've tried generating mesh beforehand, the error is still there. Jan 6 at 6:20
• I think the problem comes from "SameQ". It is very specific what it considers to be True. Jan 6 at 8:40

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.

• @xzczd what your all think, would you expect this type of information in the documentation?, If yes, where would you place it? Maybe a possible issue for NDSolve? Jan 6 at 9:06
• This seems to be the correct direction. See my answer. Jan 6 at 9:19
• @user21 I think possible issue of NDSolve is the place to go. Jan 6 at 9:22
• @xzczd, I have added something to the LinearSolve and NDSolve possible issues sections. Jan 6 at 13:31
• @Henrik How can you explain my example? :) Jan 6 at 15:32

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?

• Aha. Yes. I recall that Pardiso had some option to turn on CNR mode (CNR = conditional numerical reproducibility). @user21 One can switch it on by modifying some of the infamous iparm parameters when initializing the Pardiso solver. Have a look at this page of Pardiso's documentation and watch out for "CNR". Jan 6 at 9:37
• Alas, that will switch off parallelization of computing the nested dissection ordering. And doubt that anybody really wants to do that... Jan 6 at 9:38
• @HenrikSchumacher, thanks, I had a look. However exposing that option is not going to happen. The amount of work and unexpected trouble this will produce is not warranted for such a 'cosmetic' improvement IMHO, Jan 6 at 13:29
• "For small system e.g. n = 100, the randomness becomes 0." That makes sense. The Pardiso doc page states that randomness is introduced by only for parallel computation of the matrix reordering. And parallelization is probably truended off for too small matrices. Jan 7 at 6:52
• By the way, I think that Pardiso uses Metis for computing the reordering. I skimmed a bit through the documentation of the multithreaded version MT-Metis. Could be that the "randomness" is entirely due to race conditions. (I was somewhat chocked to read that!) The Metis authors state that this were not an issue because they also make sure that the reordering is valid. Nonetheless, the race conditions make the library quite nondeterministic. Jan 7 at 6:55

Contra example

Needs["NDSolveFEM"]
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.}

• If MaxCellMeasure -> 10^-2 is removed, the random error is there again. Might be because the error is smaller than $MachineEpsilon for very fine mesh. Still, this is an answer to my last question. (+1) Jan 6 at 6:48 • I am not a 100% sure but I believe that the mesh generation has some degree of randomness to seed the algorithms. Jan 6 at 6:58 • "Might be because the error is smaller than$MachineEpsilon for very fine mesh." I retract this line, because with MaxCellMeasure -> .001 the randomness shows up again, so MaxCellMeasure -> 10^-2 is actually too coarse… Jan 7 at 2:39
• @xzczd Nevertheless we have same answer for every run. Jan 9 at 12:45
• You mean for MaxCellMeasure -> .001 the error is constant? It's not the case on my laptop (Win 10) and Wolfram Cloud (Linux): wolframcloud.com/obj/a4938367-256c-4772-a0d2-11baa90c095a Jan 9 at 14:01