Using ParallelTable (as OP suggested in comments):

sol = ParallelTable[
   Solve[{x == n1^2 + n2^2 + n3^2, x == n4^2 + n5^2 + n6^2, 
     x == n7^2 + n8^2 + n9^2, x == n1^2 + n4^2 + n7^2, 
     x == n2^2 + n5^2 + n8^2, x == n3^2 + n6^2 + n9^2, 
     x == n1^2 + n5^2 + n9^2}, {n1, n2, n3, n4, n5, n6, n7, n8, n9}, 
    PositiveIntegers], {x, 0, 2000}];

Map[Select[Length[#] == 9 &]@*DeleteDuplicates, 
  Values[sol], {2}] //. {} -> Nothing

We Map at 2nd level (at all sub-lists of sol): composition (@*) of deleting repeated values (DeleteDuplicates) and slecting sub-lists with length of 9 (all distinct values). Finally one can delete empty lists by repeatedly applying rule {} -> Nothing and see if something is left in the end.

Post here my comment as an answer (actually I'm not sure this is the right answer). Define table of x values (just as given by OP):

    xvals = ParallelTable[
   If[TrueQ[
     Length[Select[
        Select[PowersRepresentations[n, 3, 2], Times @@ # != 0 &], 
        Length[#] == Length[Union[#]] &]] >= 6], n, Nothing], {n, 0, 
    10000}];

Now we can use this to substitute for x when solving the given system of equations, and to speed-up the code we can use ParallelMap:

ParallelMap[
  Solve[{# == n1^2 + n2^2 + n3^2, # == n4^2 + n5^2 + n6^2, # == 
      n7^2 + n8^2 + n9^2, # == n1^2 + n4^2 + n7^2, # == 
      n2^2 + n5^2 + n8^2, # == n3^2 + n6^2 + n9^2}, {n1, n2, n3, n4, 
     n5, n6, n7, n8, n9}, PositiveIntegers] &, xvals[[;; 50]]];

This gives the solution in terms of positive integer numbers, and because of time and memory consumption one can use several ranges of xvals (here the first 50 values are used).

One can also solve for particular x:

x = 3051;
Solve[{x == n1^2 + n2^2 + n3^2, x == n4^2 + n5^2 + n6^2, 
  x == n7^2 + n8^2 + n9^2, x == n1^2 + n4^2 + n7^2, 
  x == n2^2 + n5^2 + n8^2, x == n3^2 + n6^2 + n9^2}, {n1, n2, n3, n4, 
  n5, n6, n7, n8, n9}, PositiveIntegers]

This gives many solutions and {n1 -> 29, n2 -> 41, n3 -> 23, n4 -> 1, n5 -> 37, n6 -> 41, n7 -> 47, n8 -> 1, n9 -> 29} is among them.

Post here my comment as an answer (actually I'm not sure this is the right answer). Define table of x values (just as given by OP):

    xvals = ParallelTable[
   If[TrueQ[
     Length[Select[
        Select[PowersRepresentations[n, 3, 2], Times @@ # != 0 &], 
        Length[#] == Length[Union[#]] &]] >= 6], n, Nothing], {n, 0, 
    10000}];

Now we can use this to substitute for x when solving the given system of equations, and to speed-up the code we can use ParallelMap:

ParallelMap[
  Solve[{# == n1^2 + n2^2 + n3^2, # == n4^2 + n5^2 + n6^2, # == 
      n7^2 + n8^2 + n9^2, # == n1^2 + n4^2 + n7^2, # == 
      n2^2 + n5^2 + n8^2, # == n3^2 + n6^2 + n9^2}, {n1, n2, n3, n4, 
     n5, n6, n7, n8, n9}, PositiveIntegers] &, xvals[[;; 50]]];

This gives the solution in terms of positive integer numbers, and because of time and memory consumption one can use several ranges of xvals (here the first 50 values are used).

One can also solve for particular x:

x = 3051;
sol=Solve[{x == n1^2 + n2^2 + n3^2, x == n4^2 + n5^2 + n6^2, 
  x == n7^2 + n8^2 + n9^2, x == n1^2 + n4^2 + n7^2, 
  x == n2^2 + n5^2 + n8^2, x == n3^2 + n6^2 + n9^2}, {n1, n2, n3, n4, 
  n5, n6, n7, n8, n9}, PositiveIntegers]

This gives many solutions and {n1 -> 29, n2 -> 41, n3 -> 23, n4 -> 1, n5 -> 37, n6 -> 41, n7 -> 47, n8 -> 1, n9 -> 29} is among them.

To check if all the solutions are unique we can do:

Select[DeleteDuplicates /@ Values[sol], Length[#] == 9 &]

In the case of x = 3051 we have repeated values (1, 29, 41) and this Select command gives empty list as an output {}, that means there is no unrepeated values.