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): 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 ?

[forrestglines]

I'm trying to understand the setup here

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.
You mean that this removes the grid artifacts? But it's not the setup you want to use, you would like to use u(IDN, k, j, i) = 1/(1+r); but you get the grid artifacts above, right?

I'm wondering if this is due to how the integral for HSE is constructed and interpolated. I chose to do the integration separately for each mesh block so that the HSE curve for small inner mesh blocks and large outer mesh blocks are resolved. However, since the r-coordinates of the integral for each mesh block is different, I think this could create discontinuities at mesh block boundaries.

Maybe the discontinuity disappears when the integral lines up with an axis, which smells like an interpolation issue. Is the plot through the origin? If you offset from x=0 do you see artifacts at all boundaries, indicating only the axes are free from the issue?

If you raise the parameter r_sampling do the artifacts diminish? Try 100 or even 1000.

I could also try hacking the HSE integral to be the same for all mesh blocks. If it works better, we can switch main to that scheme.

[pgrete]

I'm trying to understand the setup here

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.
You mean that this removes the grid artifacts? But it's not the setup you want to use, you would like to use u(IDN, k, j, i) = 1/(1+r); but you get the grid artifacts above, right?

I'm wondering if this is due to how the integral for HSE is constructed and interpolated. I chose to do the integration separately for each mesh block so that the HSE curve for small inner mesh blocks and large outer mesh blocks are resolved. However, since the r-coordinates of the integral for each mesh block is different, I think this could create discontinuities at mesh block boundaries.

Maybe the discontinuity disappears when the integral lines up with an axis, which smells like an interpolation issue. Is the plot through the origin? If you offset from x=0 do you see artifacts at all boundaries, indicating only the axes are free from the issue?

Yes, is centered on the origin and we haven't checked yet what happens if it's slightly misaligned. @mfournier01 do you have data readily available to check what a slice looks like at x!=0?

If you raise the parameter r_sampling do the artifacts diminish? Try 100 or even 1000.

We've tried to increase it, but not that dramatically. @mfournier01 would you mind giving this a shot, too?

I could also try hacking the HSE integral to be the same for all mesh blocks. If it works better, we can switch main to that scheme.

I'm not sure this is really necessary. So far it's been working nicely (and it still is within blocks). The main motivation might be to be able to use MeshData init versus individual MeshBlockData.

[mfournier01]

Hi and thanks for replying,

Yes, is centered on the origin and we haven't checked yet what happens if it's slightly misaligned. @mfournier01 do you have data readily available to check what a slice looks like at x!=0?

I just checked that and the grid looks quite the same everywhere along the x-axis.

If you raise the parameter r_sampling do the artifacts diminish? Try 100 or even 1000.

We've tried to increase it, but not that dramatically. @mfournier01 would you mind giving this a shot, too?

Yes I tried it already up to 1000 and did not noticed any significant change.

[mfournier01]

Hey, I don't know whether this is connected to the problem mentioned above but I noticed similar patterns between meshblocks when turning static refinement on (notice that for this example I am not applying any weight to the magnetic potential, the initialization routines are main branch ones, untouched) :

The value for div B * cell_size / B is roughly 1 within these patterns, and << 1 outside. Have you ever encountered this ?

[pgrete]

The value for div B * cell_size / B is roughly 1 within these patterns, and << 1 outside. Have you ever encountered this ?

Yes, this is a plotting artifact from yt as interpolation/smoothing between levels for derivatives is currently not supported for Parthenon data (I'm not even sure it's supported in general), so there'll always be artifacts across level boundaries.

For example, here's the quantify plotted with yt based on the magnetic fields in the output file

and here based on calculating the quantity inside AthenaPK directly

[pgrete]

I also reviewed the original issue.
It seems that there's a fairly easy fix: Using a "global" profile, see e.g., c83c6aa

So for your pgen you might want to go that route for now, i.e., not generating a profile relative to the local block indices, but use some global bounds (for smallest and largest radius and number of sample points -- adjusted to the dimensions in the input file).