Assuming that you calculate the Power Spectral Entropy as per the first link I've added to the question:
$${\rm PSE} = -\sum_{i=1}^N p_i\ln p_i $$
then the spectral flatness (using the same notation) would be
$$
{\rm Flatness} = \exp\left(\frac{1}{N}\sum_{i=1}^{N} \ln p_i\right)
$$
because the $p_i$ are normalized to unit sum, so the denominator in the second linked to question is one.
As you can see, the two are similar in form except that Flatness appears to a) exponentiated and b) assume that the $p_i$ are uniform.
As to which one to use, I'd probably opt for the PSE rather than Flatness. It appears to make more use of the spectral information.
To test my theory, I've written a comparison that starts off with a purely harmonic signal and goes to a purely white noise signal and calculates both measures for each blend.
There is very little difference between the two measures, as you can see from the closeness of the red and blue curves in the bottom plot.
Note that I've normalize both curves so they go between 0 and 1, just to make it easier to compare them.
Bottom line: Choose the simplest; there appears to be little difference between the two for this example.

R Code Below
#30534
T <- 1000
white_noise <- rnorm(T,0,1)
K <- array(seq(1,5),c(5,1)) %*% array(1,c(1,T))
omega <- 2*pi*0.053492384
t <- array(seq(0,T-1), dim=c(1,T))
t <- array(1,c(5,1)) %*% t
harmonic <- colSums(sin(omega*K*t))
spectral_entropy <- function(signal)
{
psd <- abs(fft(signal))^2
psd <- psd/sum(psd)
pse <- 0
for (ix in seq(1,length(psd)))
{
pse <- pse - psd[ix]*log(psd[ix])
}
return(pse)
}
flatness <- function(signal)
{
psd <- abs(fft(signal))^2
psd <- psd/sum(psd)
flatness <- 0
for (ix in seq(1,length(psd)))
{
flatness <- flatness + log(psd[ix])/length(psd)
}
return (exp(flatness))
}
P <- 100
pse <- rep(0,P)
flt <- rep(0,P)
for (alpha in seq(0,P))
{
signal <- white_noise*(alpha/P) + (1-alpha/P)*harmonic
pse[alpha+1] <- spectral_entropy(signal)
flt[alpha+1] <- flatness(signal)
}
normalize <- function(signal)
{
return( (signal - min(signal))/(max(signal) - min(signal)))
}
layout(matrix(c(1,2,3,3), 2, 2, byrow = TRUE))
plot(abs(fft(harmonic)), type="l", col="blue")
title("Harmonic Spectrum")
plot(abs(fft(white_noise)), type="l", col="blue")
title("White Noise Spectrum")
plot(seq(0,P)/P, normalize(pse), type="l", col="blue")
lines(seq(0,P)/P, normalize(flt), type="l", col="red")
title('Comparison of PSE vs Flatness')