Troubling artifacts appear when strong mean magnetic field is aligned with grid

[dfielding14]

Hi everyone

I have been running some driven turbulence isothermal MHD simulations and have been finding some pretty bad artifacts when the simulations have a strong mean field aligned with the grid, and I wanted to see if anyone else has seen anything similar or has any thoughts.

For reference the simulations have Mach = 1 and beta = 0.1—more specifically rho = 1, dedt = 1, and |B| = sqrt(2/0.1). The magnetic field is initialized uniform in some direction. The problems also appear in cases with weaker driving and in simulations with somewhat weaker B. The simulations are 256^3, but the problems do not appear to go away at higher resolution as I have seen them in sims with resolutions as high as 2048^3.

Here is an example that shows 2D slices along the mean B direction (first movie) and perpendicular to it (second movie). In this case the mean B is perfectly aligned with the z-axis. The panels show density, mach number, velocity, electric current, plasma beta, and magnetic field. You can see the artifacts most clearly in the current, but they are present in the velocity and density too if you look closely. (note these sims also have an added small tangled field component in the ICs, but i've tried removing that and it makes no difference)

multislice.2d.yaxis.IsoMHD_256_mean_beta01.webm

multislice.2d.zaxis.IsoMHD_256_mean_beta01.webm

I have tried using xorder=2, 2c, 3, and 3c, and integrator = vl2, rk3, and rk4, and nothing changes. I have also tried aligning the B field with the x and y axis and again nothing changes. Tilting the B field in a plane so that it is aligned with (0,1,1), or tilting it by some small amount so it is aligned with (1, 0.1, 0.1) also doesn't help. The only thing that seems to help is having the mean B point in the (1,1,1) direction, which you can see in the following example

multislice.2d.yaxis.IsoMHD_256_mean_beta01_diag_rk3.webm

I am fairly worried about this and I know that others (e.g., philipp kempski and @ltlancas although I think Lachlans simulation are actually with athena4.2) have seen the same issue, so I would love to hear what people think.

[pdmullen]

Does the issue persist with an LLF Riemann solver? I may be off base here, but what you are describing sounds reminiscent of Carbuncle... I wonder if something like Athena C-Version's H-correction could help here?

[ltlancas]

I've been meaning to ask a wider audience about this so thanks @dfielding14 for bringing it up! To expand on the problems that I was having, I was running stellar wind-driven bubble simulations where the background medium has a grid-aligned mean magnetic field with beta ~ 0.01 that is then turbulently driven before the wind-driven bubble is started, so there isn't necessarily a strong grid-aligned field everywhere when the bubble starts expanding. That said, since the field is draped around the bubble as it expands, the bubble does tend to create a strong, grid-aligned field in certain places as it expands, and these spots are preferentially picked up when taking slices through the mid-plane, so very well could be the same thing.

At first I thought that this could have something to do with my method for injecting wind energy, though I tried a bunch of different things along that front and convinced myself that that isn't the problem. Happy to go in to more detail on why if needed.

My default run uses the Roe Solver, the VL integrator, and H-correction (relevant to @pdmullen 's comment). I tried the CTU integrator and HLLD solver and always got these same artifacts. The artifacts didn't come up in simulations I ran without cooling, which I would guess is because the lack of cooling prevents the collapse of the shell around the bubble and therefore prevents an amplification of the field that comes with that. I haven't tried an LLF solver.

[c-white]

I do wonder if this is an issue with constructing cell-centered fields. For example, for there is the idea that one should interpolate B^2 from faces to cells, in order to accurately subtract magnetic energy from total energy before variable inversion. There are papers on this, especially for relativity, probably out of Valencia, but I can't find them. This also wouldn't be easy to test.

If not for the effects on velocity, I wouldn't be convinced that there is a problem, given the trickiness of computing j entirely accurately. It should probably be calculated on edges from face-centered B, if I had to guess.

Given that the (1,1,1) orientation helps when (0,1,1) doesn't, perhaps having many Riemann problems with weak parallel B but strong perpendicular B (in at least one direction) is the problem.

[dfielding14]

I just tried the aligned with the grid case using HLLE and xorder=2 & 3, and in both cases I don't see any issues.

[c-white]

When all else fails, add diffusion I guess.

[pdmullen]

The solution to diffusion is resolution.
The resolution to the solution is diffusion.