Divergence-preserving prolongation of coarse/fine shared face-centered fields for 3D spherical polar grids

[mariareneemf]

Description

I changed the `ProlongateSharedFieldX` functions in `mesh_refinement.cpp` to do the interpolation that defines the fine face-centered fields considering the areas of the overlapping fine and coarse cells for 3D spherical polar grids. To have access to the necessary coarse areas, I allocated memory for the `coord_area_x` arrays in the spherical polar file. These were previously skipped for coarse meshes with AMR. The coarse area arrays were defined in `mesh_refinement.hpp`. I modified the field_loop_poles.cpp problem generator to include a divergence check of the magnetic field. This was used for testing and may not need to be part of a merge. I included the input file `athinput.field_loop_poles_checkdivb` used for this test.

The current AMR implementation does not preserve the divergence of the magnetic field during refinement in spherical polar grids. Following a prescription suggested by @tomidakn, I have modified the interpolation that sets the fine shared face fields to enforce the following constraints:

1. The transverse gradients on the fine cells match the (slope-limited) ones on the coarse cells.
2. The total magnetic flux on the fine cells matches that on the coarse cell. With these changes, the divergence is preserved during refinement.

• Interpolation in ProlongateSharedFieldX functions now preserves the divergence and magnetic flux during refinement of 3D spherical polar grids.
• Added divergence check to field loop poles problem generator and an input file to test it.

Testing and validation

I ran the field loop poles test with AMR and checked the divergence of the magnetic field. To do this, I added a refinement condition in the problem generator to force refinement at a specific cycle in meshblocks with a magnetic field above a threshold and added a history output of the total divergence in the simulation. With the current prolongation scheme for shared face fields, the divergence increased by several orders of magnitude (~9) during refinement, as it does not preserve the magnetic flux. After the above-mentioned changes, the divergence change was less than an order of magnitude during refinement compared to the evolution of the divergence with no AMR. The divergence diagnostic was scaled with the number of cells used to compute it. This test can be run using the problem generator field_loop_poles.cpp with the input file athinput.field_loop_poles_checkdivb. Note that it is necessary to suppress the AMR-polar boundary error flag in bvals/utils/check_polar.cpp to run this test.

I ran the field loop poles test. I compiled the code with gcc 4.8.5, testing both, with openmpi 3.1.3 and no mpi.

The field loop poles test has a polar boundary, therefore, to run the field loop poles test it is necessary to comment L48-52 in bvals/utils/check_polar.cpp to suppress the AMR-polar boundary error. The test did not refine meshblocks adjacent to the pole so it was not a problem. I reverted this comment in check_polar.cpp as it cannot be generally applied.

To-do

The changes are currently only implemented for 3D but may be useful to implement in 2D and 1D and for other non-uniform coordinate systems. The changes should apply to any curvilinear coordinate system but are only triggered for spherical_polar coordinates since I have not tested them with other curvilinear coordinates. Additionally, it may be useful to modify the polar boundary or refinement scheme to support AMR in spherical polar grids with polar boundaries. Finally, only the r coordinate used in the prolongation of shared face fields is area-averaged (L125 src/coordinates/spherical_polar.cpp). Face-area averaged coordinates could be also used for X2 and X3, although the difference may be very small.

[buildbot-princeton]

Can one of the admins verify this patch?

[felker]

ok to test

[felker]

Thanks for contributing! Still looking through it so I can understand the changes.

With the current prolongation scheme for shared face fields, the divergence increased by several orders of magnitude (~9) during refinement, as it does not preserve the magnetic flux. After the above-mentioned changes, the divergence change was less than an order of magnitude during refinement compared to the evolution of the divergence with no AMR

Shouldn't the divergence be zero within machine precision after this fix? Can you check/ report in absolute values vs. relative order of magnitudes?

I had completely forgotten that AMR (and also SMR?) + curvilinear coordinates causes the divergence constraint to be violated. A bit of an edge case, although the Stone (2020) method paper has such a test in Figure 31 for the MHD blast wave with AMR in spherical polar coordinates. The violation of the constraint is not mentioned in the paper, nor anywhere in the Wiki. We should document this and any expecations (and possibly emit a warning when MHD is enabled in the code in such cases?).

Any though to turning your inputs/mhd/athinput.field_loop_poles_checkdivb and changes to src/pgen/field_loop_poles.cpp into a Python-based regression test?

[tomidakn]

And indeed, for better consistency, it is probably better to use the area-averaged centers for the other directions, although I agree with you that the difference will be small.

[felker]

I accidentally committed your commits directly to the master branch on the command line when fixing some style things. I cannot easily revert the pseudo-merge nor re-open this PR. Can you open another PR containing @tomidakn's requested fixes? We can continue the discussion there. Sorry!

[mariareneemf]

Shouldn't the divergence be zero within machine precision after this fix? Can you check/ report in absolute values vs. relative order of magnitudes?

Yes, this plot illustrates it better. I am plotting the total sum of the absolute value of the divergence normalized by dividing by the magnitude of the magnetic field and multiplying by the minimum cell size. This quantity is then divided by the total number of cells that are used in computing it. The divergence diagnostic is preserved with the new changes.

Any though to turning your inputs/mhd/athinput.field_loop_poles_checkdivb and changes to src/pgen/field_loop_poles.cpp into a Python-based regression test?

I am not sure if it is optimal since the running of this test requires the suppression of the AMR-polar boundary error flagged in src/bvals/utils/check_polar.cpp. This error flag can't be generally removed since it could result in meshblocks adjacent to the pole having different refinement levels, which raises another error.