HSE init in ghost zone

[pgrete]

@mfournier01 discovered an issue when applying the (initialized) density field as a weight to the scalar potential that initialized the magnetic field.

The MWE is the default cluster.in (just with periodic BC and perturbations): cluster.txt with these changes (so loop bounds are already adjusted to also init one layer into the ghost cells): https://github.com/parthenon-hpc-lab/athenapk/compare/main...pgrete/hse-ghost-mwe

Note, that using

u(IDN, k, j, i) = 1/(1+r);

instead of

u(IDN, k, j, i) = rho_r;

results in a correctly initialized field.

The magnetic field divergence for the initial dump looks like:

slc = yt.SlicePlot(ds,"x",("gas","magnetic_field_divergence_absolute"))
slc.annotate_grids()
slc.show()

so there are artifacts at the block boundaries (and more specifically, it looks like at boundaries where the coordinate do not cross 0). (Note, the artifacts at the outermost edge are related to yt domain boundary handling, so should not be a concern.)

We went through the hse code, but didn't spot anything that looks like being related to this kind of behavior. Any ideas, @forrestglines ?