Is your code valgrind clean?
Yes, I also initialize all allocations with NaNs to be sure I'm not using anything uninitialized.
We can try and test this. Replace your MatMFFD with an actual matrix and run. Do you see any variability?
I think I did what you're asking. I have -snes_mf_operator set, and then SNESSetJacobian(snes, diag_ones, diag_ones, NULL, NULL) where diag_ones is a matrix with ones on the diagonal. Two runs below, still with differences but sometimes identical.
0 SNES Function norm 3.424003312857e+04
0 KSP Residual norm 3.424003312857e+04
1 KSP Residual norm 2.871734444536e+04
2 KSP Residual norm 2.490276930242e+04
3 KSP Residual norm 2.131675872968e+04
4 KSP Residual norm 1.973129814235e+04
5 KSP Residual norm 1.832377856317e+04
6 KSP Residual norm 1.716783617436e+04
7 KSP Residual norm 1.583963149542e+04
8 KSP Residual norm 1.482272170304e+04
9 KSP Residual norm 1.380312106742e+04
10 KSP Residual norm 1.297793480658e+04
11 KSP Residual norm 1.208599123244e+04
12 KSP Residual norm 1.137345655227e+04
13 KSP Residual norm 1.059676909366e+04
14 KSP Residual norm 1.003823862398e+04
15 KSP Residual norm 9.425879221354e+03
16 KSP Residual norm 8.954805890038e+03
17 KSP Residual norm 8.592372470456e+03
18 KSP Residual norm 8.060707175821e+03
19 KSP Residual norm 7.782057728723e+03
20 KSP Residual norm 7.449686095424e+03
Linear solve converged due to CONVERGED_ITS iterations 20
KSP Object: 1 MPI process
type: gmres
restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
happy breakdown tolerance 1e-30
maximum iterations=20, initial guess is zero
tolerances: relative=0.1, absolute=1e-15, divergence=10.
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI process
type: none
linear system matrix followed by preconditioner matrix:
Mat Object: 1 MPI process
type: mffd
rows=16384, cols=16384
Matrix-free approximation:
err=1.49012e-08 (relative error in function evaluation)
Using wp compute h routine
Does not compute normU
Mat Object: 1 MPI process
type: seqaij
rows=16384, cols=16384
total: nonzeros=16384, allocated nonzeros=16384
total number of mallocs used during MatSetValues calls=0
not using I-node routines
1 SNES Function norm 1.085015646971e+04
Nonlinear solve converged due to CONVERGED_ITS iterations 1
SNES Object: 1 MPI process
type: newtonls
maximum iterations=1, maximum function evaluations=-1
tolerances: relative=0.1, absolute=1e-15, solution=1e-15
total number of linear solver iterations=20
total number of function evaluations=23
norm schedule ALWAYS
Jacobian is never rebuilt
Jacobian is applied matrix-free with differencing
Preconditioning Jacobian is built using finite differences with coloring
SNESLineSearch Object: 1 MPI process
type: basic
maxstep=1.000000e+08, minlambda=1.000000e-12
tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
maximum iterations=40
KSP Object: 1 MPI process
type: gmres
restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
happy breakdown tolerance 1e-30
maximum iterations=20, initial guess is zero
tolerances: relative=0.1, absolute=1e-15, divergence=10.
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI process
type: none
linear system matrix followed by preconditioner matrix:
Mat Object: 1 MPI process
type: mffd
rows=16384, cols=16384
Matrix-free approximation:
err=1.49012e-08 (relative error in function evaluation)
Using wp compute h routine
Does not compute normU
Mat Object: 1 MPI process
type: seqaij
rows=16384, cols=16384
total: nonzeros=16384, allocated nonzeros=16384
total number of mallocs used during MatSetValues calls=0
not using I-node routines
0 SNES Function norm 3.424003312857e+04
0 KSP Residual norm 3.424003312857e+04
1 KSP Residual norm 2.871734444536e+04
2 KSP Residual norm 2.490276931041e+04
3 KSP Residual norm 2.131675873776e+04
4 KSP Residual norm 1.973129814908e+04
5 KSP Residual norm 1.832377852186e+04
6 KSP Residual norm 1.716783608174e+04
7 KSP Residual norm 1.583963128956e+04
8 KSP Residual norm 1.482272160069e+04
9 KSP Residual norm 1.380312087005e+04
10 KSP Residual norm 1.297793458796e+04
11 KSP Residual norm 1.208599115602e+04
12 KSP Residual norm 1.137345657533e+04
13 KSP Residual norm 1.059676906197e+04
14 KSP Residual norm 1.003823857515e+04
15 KSP Residual norm 9.425879177747e+03
16 KSP Residual norm 8.954805850825e+03
17 KSP Residual norm 8.592372413320e+03
18 KSP Residual norm 8.060706994110e+03
19 KSP Residual norm 7.782057560782e+03
20 KSP Residual norm 7.449686034356e+03
Linear solve converged due to CONVERGED_ITS iterations 20
KSP Object: 1 MPI process
type: gmres
restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
happy breakdown tolerance 1e-30
maximum iterations=20, initial guess is zero
tolerances: relative=0.1, absolute=1e-15, divergence=10.
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI process
type: none
linear system matrix followed by preconditioner matrix:
Mat Object: 1 MPI process
type: mffd
rows=16384, cols=16384
Matrix-free approximation:
err=1.49012e-08 (relative error in function evaluation)
Using wp compute h routine
Does not compute normU
Mat Object: 1 MPI process
type: seqaij
rows=16384, cols=16384
total: nonzeros=16384, allocated nonzeros=16384
total number of mallocs used during MatSetValues calls=0
not using I-node routines
1 SNES Function norm 1.085015821006e+04
Nonlinear solve converged due to CONVERGED_ITS iterations 1
SNES Object: 1 MPI process
type: newtonls
maximum iterations=1, maximum function evaluations=-1
tolerances: relative=0.1, absolute=1e-15, solution=1e-15
total number of linear solver iterations=20
total number of function evaluations=23
norm schedule ALWAYS
Jacobian is never rebuilt
Jacobian is applied matrix-free with differencing
Preconditioning Jacobian is built using finite differences with coloring
SNESLineSearch Object: 1 MPI process
type: basic
maxstep=1.000000e+08, minlambda=1.000000e-12
tolerances: relative=1.000000e-08, absolute=1.000000e-15, lambda=1.000000e-08
maximum iterations=40
KSP Object: 1 MPI process
type: gmres
restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement
happy breakdown tolerance 1e-30
maximum iterations=20, initial guess is zero
tolerances: relative=0.1, absolute=1e-15, divergence=10.
left preconditioning
using PRECONDITIONED norm type for convergence test
PC Object: 1 MPI process
type: none
linear system matrix followed by preconditioner matrix:
Mat Object: 1 MPI process
type: mffd
rows=16384, cols=16384
Matrix-free approximation:
err=1.49012e-08 (relative error in function evaluation)
Using wp compute h routine
Does not compute normU
Mat Object: 1 MPI process
type: seqaij
rows=16384, cols=16384
total: nonzeros=16384, allocated nonzeros=16384
total number of mallocs used during MatSetValues calls=0
not using I-node routines