Overhaul excision treatment - [merged]

[jmstone]

In GitLab by @pdmullen on Mar 28, 2022, 16:36

_Merges excisionv3 -> master

This PR overhauls the excision treatment.

Summary of Changes

• Excision criterion is now dependent on spherical KS radius r_ks (lower r), not capital R. Effectively, this means that the logically excised region is an ellipsoid for a spinning black hole.
• Metric evaluation switches from flooring R to flooring r_ks (i.e., r_ks floored to 1.0). Furthermore, we have eliminated the floor on z (see later To Do).
• Add LLF GR Hydro/MHD Riemann solvers (taken/modified from @jmstone216's FOFC branch).
• Introduce cell-centered and face-centered masks (cc_mask and fc_mask, respectively).
  • The cc_mask is used to check what zones need to have their primitive density, pressure, and velocities reset to <coord>/dexcise, <coord>/pexcise, and 0, respectively, in ConsToPrim. (NOTE: the excised region is no longer reset to initial condition values). The cc_mask criterion follows: if r_ks evaluated at this cell-center is <= 1, mask the cell.
  • The fc_mask signals to replace a face's flux to a first-order (donor cell) LLF flux. The fc_mask criterion follows: if r_ks evaluated at this face-center is <=1, or if any other portion of grid cells sharing this face (i.e., cell centers, cell faces, cell edges, cell nodes) is <=1, mask the cell (added complexity here as two neighboring cells share a face).
• Add UserDefinedBCs for the gr_torus problem. These boundary conditions are "quasi"-no-inflow---that is, they zero the primitive velocity if it is inwardly directed. The primitive velocity is not the advection velocity. @jmstone216, I'd appreciate another set of eyes on this part of the code in particular.
• Fixes the UserDefinedBCs in the gr_bondi problem and turns this pgen into a regression test. For now, this test is only performed on GPUs, otherwise, the regression tests take too long.

Results

The above changes are what I found to be required in order to run a 6--12 spin 0.9375 torus to 10,000 M with rk2+ppm (see above for a 3D rendering and poloidal slice for a 256^3 calculation at 10,000 M). Note that Athena++ does not have many of the changes described above...I believe they can successfully run these tori because of the more stringent floors/ceilings (e.g., sigma_max, beta_min, gamma_max, etc.) employed in their ConsToPrim algorithm. Our integration to 10,000 M with this MR employed only a density floor and pressure floor. Eventually, I think we will later need to better our floors/ceilings (in particular, I believe we should switch our density and pressure floors to be operating in the normal observer frame. Any thoughts on this, @emost?)

I have also tested the Bondi problem with our new excision treatment. Notably, this MR spoils the beautiful second order convergence we observed for the Bondi problem in !13 (above left: !13, above right: this MR). The differences in the errors is due to the behavior of cells about the excised region in this MR. In !13, we used an excision radius of R=2.5 and we reset the excised region to initial condition values. In this MR, we use an excision radius of r=1.0 (note r=R for spin a=0) and we reset the excised region to floor values. The plot on the right is used as a rubric for the Bondi test now included in CI.

To Do

1. Revamp floors. Potentially, a better treatment of floors would make integrating tori much more robust (or better yet, much less fragile).
2. Keep in the back of our minds: WhiskyMHD actually extrapolates state vectors from the innermost zones into the excised region (their excision region is cubical), rather than setting them directly to floors. Perhaps we should also consider something like this. As implemented here, reconstruction stencils that include zones inside the excision region will see a sharp discontinuity in primitive densities, pressures, and velocities.
3. I have removed a floor on z in the metric evaluation (I hope this won't cause troubles down the road). I believe that, in the future, we should consider eliminating the r_ks = 1.0 floor. If we do this, we will likely have to re-introduce the floor on z for r_ks < 1.
4. There are seemingly still issues with SMR. Things still crash with NaN, even after 45bad6f422cc381a95286d35288381476cdf2859.
5. Test @gnwong's monopole problem.
6. Put our 6--12 spin=0.9375 calculation through the EHT code comparison analysis (requires adding flux calciations; e.g., mass accretion rates, angular momentum fluxes, etc.).

[jmstone]

In GitLab by @pdmullen on Mar 28, 2022, 16:36

requested review from @jmstone216

[jmstone]

In GitLab by @pdmullen on Mar 28, 2022, 16:57

added 1 commit

• 88cdcb87 - Fix typos in function briefs

Compare with previous version

[jmstone]

In GitLab by @pdmullen on Mar 28, 2022, 17:13

added 1 commit

• 342a9b5f - Further separate excision treatment from rest of flux calculation

Compare with previous version

[jmstone]

In GitLab by @pdmullen on Mar 29, 2022, 14:43

marked this merge request as ready

[jmstone]

In GitLab by @jmstone216 on Mar 31, 2022, 13:50

approved this merge request

[jmstone]

In GitLab by @jmstone216 on Mar 31, 2022, 13:51

This all looks great to me. A huge amount of work, so thanks for taking the lead on it. It is a solid foundation on which to build better floors and add FOFC.

[jmstone]

In GitLab by @jmstone216 on Mar 31, 2022, 13:51

mentioned in commit 87e89858d3a53dadee23d6c92b096afd666adc6b