Thanks sir. I just tried to run tutorial mixed space code, it is confusing to see why running in parallel produces different result compared to running it normally. As the increasing the processor the values further changing. Could anyone please tell me, is there any specific reasons for the same. The example code is

```
import os
from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.mesh import CellType, create_unit_square
from dolfinx.nls.petsc import NewtonSolver
from ufl import dx, grad, inner
lmbda = 1.0e-02 # surface parameter
dt = 5.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicholson
msh = create_unit_square(MPI.COMM_WORLD, 4, 4, CellType.triangle)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))
q, v = ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
c, mu = ufl.split(u)
c0, mu0 = ufl.split(u0)
# Zero u
u.x.array[:] = 0.0
# Interpolate initial condition
np.random.seed(30)
u.sub(0).interpolate(lambda x: 0.25 + 0.01 * (0.5 - np.random.rand(x.shape[1])))
u.x.scatter_forward()
# Compute the chemical potential df/dc
c = ufl.variable(c)
f = 100 * c**2 * (1 - c)**2
dfdc = ufl.diff(f, c)
# mu_(n+theta)
mu_mid = (1.0 - theta) * mu0 + theta * mu
# Weak statement of the equations
F0 = inner(c, q) * dx - inner(c0, q) * dx + dt * inner(grad(mu_mid), grad(q)) * dx
F1 = inner(mu, v) * dx - inner(dfdc, v) * dx - lmbda * inner(grad(c), grad(v)) * dx
F = F0 + F1
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, u)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
ksp.setFromOptions()
# Step in time
t = 0.0
timesteps = [0]
# Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 3 * dt
else:
T = 1 * dt
# Get the sub-space for c and the corresponding dofs in the mixed space
# vector
V0, dofs = ME.sub(0).collapse()
c = u.sub(0)
u0.x.array[:] = u.x.array
while (t < T):
t += dt
r = solver.solve(u)
print(f"Step {int(t/dt)}: num iterations: {r[0]}")
u0.x.array[:] = u.x.array
print (u.x.array)
```

When I tried with running normal it produces values like

```
[ 0.25268046 0.25113871 0.25368124 18.6876801 18.70976266 18.66311557
0.24805195 18.78723264 0.25262463 18.6849591 0.25215371 18.69868973
0.25259528 18.69255675 0.25236828 18.69802514 0.24901909 18.77010902
0.24577314 18.84374189 0.24803946 18.79348743 0.24675552 18.82031573
0.24922422 18.76874523 0.25341396 18.67165108 0.25237948 18.69379297
0.24902894 18.77342398 0.2514791 18.72305914 0.24916687 18.77011863
0.25445935 18.63935928 0.24505966 18.86024589 0.24862446 18.78465413
0.25159698 18.71086092 0.24561697 18.84971027 0.24966788 18.75577913
0.25484011 18.6462961 ]
```

But running with processors 4, it produces different numbers, is there any difference? Even I constraint the randomness in code for initial concentration.

```
[0.25119405 0.25119405 0.2495566 0.2495566 0.2495566 0.2495566 ]
[0.24914252 0.24914252 0.24914252 0.2495566 0.2495566 0.2495566
0.25153556 0.25153556]
[0.25153556 0.25153556 0.25153556 0.25153556 0.2495566 0.24914252
0.24914252 0.24914252]
[0.2495566 0.2495566 0.2495566 0.2495566 0.2495566 0.25119405
0.25119405 0.24914252 0.24914252 0.24914252]
```

Thanks in advance for all your advices.