Use hydro speeds in GRMHD HLL where magnetic field is zero #5788
Conversation
|
|
||
| for (size_t i = 0; i < get(rest_mass_density).size(); i++) { | ||
| if (get(rest_mass_density)[i] < 1.e-08) { | ||
| double alpha = 0.0; // get(rest_mass_density)[i]/1.e-08; |
There was a problem hiding this comment.
confusion... why is alpha 0 and then we are doing math with alpha?
There was a problem hiding this comment.
@ffoucart do we still need this?
There was a problem hiding this comment.
It used to be necessary for stability of the atmosphere. Effectively, it makes the code use the speed of light as characteristic speed at low density. I have not tried to remove that correction though, and it's probably not a great idea to keep the hard-coded 1.e-8 in the final version.
My preference, if possible, would be to get rid of the alpha, just set the speeds to the speed of light at low density, and have the density threshold as an specifiable input that we can test on merger simulations and maybe lower / remove at a later time if it proves unnecessary.
|
The only difference with what I used for BNS mergers seems to be that the sound speed is now calculated within the equation of state class. Assuming that this was already reviewed separately, I'm good with this PR (up to the discussion above on the treatment of low-density regions). |
|
If you want me to try some simulations with different thresholds for the low-density regions, let me know. |
|
Okay, great! I'm making a few other minor cleanup changes. I'll add a parameter for the density cutoff for atmosphere and will just use the light speeds there. Then we can test disabling it if we want but can re-enable later. |
|
I pushed fixups. @ffoucart I changed the calculation, it didn't agree with what was in the Rezzolla book as far as I could tell |
|
The formula we had in spec comes from Font. I don't have the Rezzola book,
but I can have a look if there is a disagreement.
…On Wed, Feb 21, 2024, 17:07 Nils Deppe ***@***.***> wrote:
I pushed fixups. @ffoucart <https://github.com/ffoucart> I changed the
calculation, it didn't agree with what was in the Rezzolla book as far as I
could tell
—
Reply to this email directly, view it on GitHub
<#5788 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/ABLFTIXZJT5WV3UF6RE2MG3YU2K5XAVCNFSM6AAAAABDR5QSEOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTSNJYGQ3TSMJQGA>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
Here's the tex for the equations I switched to: https://spectre-code.org/namespaceRelativisticEuler_1_1Valencia.html#a8a63b3bfe1b17683401ba9abd4d6f7f4 Could this be related to grid vs. inertial frame solves? |
| // Compute v_dot_normal, v^i n_i | ||
| dot_product(make_not_null(&v_dot_normal), spatial_velocity, | ||
| normal_covector); | ||
| get(v_dot_normal) *= get(lapse); |
There was a problem hiding this comment.
I think our two original formulae are the same, but you are having additional lapse factors by multiplying v_dot_normal here (giving you my v0), then still multiplying, but using this as v_n in the discriminant and still multiplying by the lapse in the final expression. In your formulation, I think you don't want L147. Then you should get an expression equivalent to mine, if I am not mistaken.
There was a problem hiding this comment.
You're right, I've removed this!
nilsdeppe
force-pushed
the
mhd_hll_sound_speed
branch
from
6a00ff9
to
face522
Compare
There was a problem hiding this comment.
Looks good to me. As a sanity check, you could run the first few hours of evolution of the lowest resolution BNS system, e.g. on wheeler:
/panfs/ds09/sxs/ffoucart/Runs/Spectre/BNS_Gamma2/ConvergenceTest/5Points
The maximum density in particular is quite sensitive to the HLL solver, so checking that you get the same result as me would be a good way to verify consistency of the new implementation.
| // Compute v^2=v^i v_i | ||
| dot_product(make_not_null(&v_squared), spatial_velocity, | ||
| spatial_velocity_one_form); | ||
| get(v_squared) = clamp(get(v_squared), 0.0, 1.0 - 1.0e8); |
| } | ||
| } | ||
|
|
||
| // Correct for mesh velocity | ||
| if (normal_dot_mesh_velocity.has_value()) { | ||
| get(*packaged_largest_outgoing_char_speed) = |
There was a problem hiding this comment.
is there a reason not to use -=?
nilsdeppe
force-pushed
the
mhd_hll_sound_speed
branch
from
face522
to
04cda3d
Compare
|
@nilsdeppe is there a reason this didn't get merged? If not, it needs a rebase... |
|
@nilsdeppe will look into this |
nilsdeppe
force-pushed
the
mhd_hll_sound_speed
branch
from
04cda3d
to
17d661e
Compare
nilsdeppe
force-pushed
the
mhd_hll_sound_speed
branch
2 times, most recently
from
311de44
to
9d68366
Compare
nilsdeppe
force-pushed
the
mhd_hll_sound_speed
branch
2 times, most recently
from
2a84901
to
cd1c011
Compare
nilsdeppe
force-pushed
the
mhd_hll_sound_speed
branch
from
cd1c011
to
a76cf6d
Compare
|
I've rebased this and made sure it agrees with the paper. Here's an image of the inspiral trajectories from a few different runs I've done trying to improve the domain decomposition and tuning other parameters. The point is they all agree. |
Proposed changes
With taggedtuple support in pypp the boundary condition python files can be cleaned up, but I don't have the energy right now to do that huge refactor...
Upgrade instructions
Code review checklist
make docto generate the documentation locally intoBUILD_DIR/docs/html.Then open
index.html.code review guide.
bugfixornew featureif appropriate.Further comments