Make GrMhd system compatible with 3D EoS #5479

isaaclegred · 2023-09-20T14:52:15Z

Proposed changes

This PR adds support for 3D thermodynamic EoSs in the evolution system. Very little is changed (another instantiation is added) except in primitive recovery, where some logic needed to be added.

Upgrade instructions

Code review checklist

The code is documented and the documentation renders correctly. Run
make doc to generate the documentation locally into BUILD_DIR/docs/html.
Then open index.html.
The code follows the stylistic and code quality guidelines listed in the
code review guide.
The PR lists upgrade instructions and is labeled bugfix or
new feature if appropriate.

Further comments

isaaclegred · 2023-09-20T14:53:47Z

@nilsvu, @nilsdeppe, (perhaps @ermost should take a look at the modifications to primitive recovery too)

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp

src/Evolution/Systems/GrMhd/ValenciaDivClean/PrimitiveFromConservative.cpp

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp

nilsdeppe · 2023-09-20T15:51:09Z

src/Evolution/Systems/GrMhd/ValenciaDivClean/NewmanHamlin.cpp

@@ -141,7 +141,12 @@ std::optional<PrimitiveRecoveryData> NewmanHamlin::apply(
        return std::nullopt;
      }
    } else if constexpr (ThermodynamicDim == 3) {
-      ERROR("3d EOS not implemented");
+      Scalar<double> local_internal_energy(


const Scalar<double>

nilsdeppe · 2023-09-20T21:43:45Z

src/Evolution/Systems/GrMhd/ValenciaDivClean/Characteristics.cpp

+            rest_mass_density, specific_internal_energy, electron_fraction);
+    get(sound_speed_squared) =
+        get(equation_of_state.sound_speed_squared_from_density_and_temperature(
+            rest_mass_density, temperature, electron_fraction))


Oh, CI is complaining about a missing semicolon on this line :)

ermost · 2023-09-21T02:39:34Z

Hi @isaaclegred,

this looks great! The Con2Prim changes in Newman and KastaunMHD look good. Do you also want to modify KastaunHydro?
(I think you could just copy the code you changed in the MHD routine.)

isaaclegred · 2023-09-25T19:13:07Z

this looks great! The Con2Prim changes in Newman and KastaunMHD look good. Do you also want to modify KastaunHydro? (I think you could just copy the code you changed in the MHD routine.)

Yep I can definitely do this, @ermost!

isaaclegred · 2023-09-25T19:18:39Z

@nilsdeppe should I rebase on develop and fix the conflict?

nilsdeppe · 2023-09-25T19:19:36Z

Yes, though maybe wait a bit because #5494 will probably also cause a conflict?

Edit: I'm currently rebasing that one

isaaclegred · 2023-09-27T17:37:35Z

@nilsdeppe and @nilsvu any other comments or should I squash?

nilsdeppe

Looks good, please squash!

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAlHydro.cpp

src/Evolution/Systems/GrMhd/ValenciaDivClean/PalenzuelaEtAl.cpp

src/Evolution/Systems/GrMhd/ValenciaDivClean/PrimitiveFromConservative.cpp

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp

nilsvu · 2023-09-27T18:16:23Z

And yes please squash

nilsvu

Changes look good, but it looks like none of the ThermodynamicDim=3 code is tested. That seems important since we'll mostly run with 3D EOSs once the other changes to the executables go through. Could you please add tests for the code added here?

ermost · 2023-10-03T22:07:11Z

For the 3D Primitive recovery test, would it be possible to check what happens when the evolved electron fraction goes out of bounds?

isaaclegred · 2023-10-04T16:02:09Z

For the 3D Primitive recovery test, would it be possible to check what happens when the evolved electron fraction goes out of bounds?

@ermost I can certainly add a test, what is the expected behavior in that case?

ermost · 2023-10-04T20:00:59Z

That would be great, thank you @isaaclegred!

It should be pushed back to the bounds of the EOS. so if you input ye_star < 0, you should get ye_min, and if you have ye_star > rho*ye_max, you should get ye_max.

The undertermined part is what happens if D <0, but that should be caught in the fix atmosphere, or fix conservative part. i.e., whether you keep ye, or just replace it by the atmospheric value.

isaaclegred · 2023-10-06T20:09:05Z

@ermost, @nilsdeppe, @nilsvu whenever you get a chance this is ready for review! I've also added some potentially unnecessary tests to GhMhd subcell related things with 3d EoSs. Let me know if there's anything else you'd like to see.

nilsdeppe

I'm happy for you to squash once @ermost has taken a look :)

nilsdeppe · 2023-10-07T00:34:35Z

src/Evolution/Systems/GrMhd/ValenciaDivClean/NewmanHamlin.cpp

+            (rho_h_w_squared / square(current_lorentz_factor) -
+             current_pressure) /
+                current_rest_mass_density -
+            1.0;


LGTM, this has a catastrophic cancellation, but I don't know what to do about that other than to rework the entire recovery scheme (not worth it!)

Yeah, I would guess we would be using Kastaun most of the time anyway, if this becomes a serious limitation we can take a look at trying to patch it up.

I think it's an intrinsic limitation of the Newman Hamlin scheme in terms of iterating over pressure. For using tabulated EOS the above has always been fine for me, when I was using Newman-Hamlin.

Yea, I think the "solution" is "use a better recovery scheme" 😅 Kastaun isn't perfect in this regard, but it isn't obvious to me that a "perfect" scheme is even possible

I think something better would probably go towards a physicality preserving evolution scheme, but that is even harder :)

Yes, that would be the dream... But agreed, very, very difficult :(

ermost · 2023-10-10T18:20:22Z

Thanks a lot for additional test @isaaclegred . This looks good to me.

isaaclegred · 2023-10-10T20:57:05Z

Just rebased on develop to fix conflicts

nilsdeppe · 2023-10-10T22:14:40Z

Uh oh, something went wrong during the rebase :(

isaaclegred · 2023-10-11T16:51:29Z

Yeah Con2Prim is broken now (like after fixing things to make it compile), I'm investigating

isaaclegred · 2023-10-12T17:57:06Z

@nilsdeppe @ermost I think this is ready to be looked at again.

ermost · 2023-10-13T05:08:23Z

Looks good! Thanks for fixing this!

nilsdeppe · 2023-10-13T13:58:35Z

LGTM, please squash :D

isaaclegred · 2023-10-13T22:54:34Z

Changed some things to get rid of clang-tidy complaints, should be good to go when CI completes @nilsdeppe and @ermost!

nilsdeppe

@nilsvu please take a look again...

nilsdeppe · 2023-10-14T01:34:28Z

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp

@@ -297,15 +323,17 @@ Primitives FunctionOfMu<EnforcePhysicality, ThermodynamicDim>::primitives(
    p_hat = get(equation_of_state_.pressure_from_density_and_energy(
        Scalar<double>(rho_hat), Scalar<double>(epsilon_hat)));
  } else if constexpr (ThermodynamicDim == 3) {
-    ERROR("3d EOS not implemented");
+    p_hat = get(equation_of_state_.pressure_from_density_and_energy(
+        Scalar<double>(rho_hat), Scalar<double>(epsilon_hat),


If you wanted to fix these, I think you need to swap the parens for braces: Scalar<double>{...}

@nilsvu said to ignore these earlier I think, which is why I didn't change them, but I can change them if that's the way things are going.

Sounds good!

nilsvu

Merging as-is is fine, or squash in these things if you want to do them and merge then.

nilsvu · 2023-10-17T20:45:44Z

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp

+      double eps_min{};
+      if constexpr (ThermodynamicDim == 3) {
+        eps_min = equation_of_state_.specific_internal_energy_lower_bound(
+            rest_mass_density_times_lorentz_factor_ / lorentz_max,
+            electron_fraction_);
+      } else {
+        eps_min = equation_of_state_.specific_internal_energy_lower_bound(
+            rest_mass_density_times_lorentz_factor_ / lorentz_max);
+      }


(optional) move this code outside the if because it's repeated in both branches?

nilsvu · 2023-10-17T20:48:17Z

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp

+        eps_min = equation_of_state_.specific_internal_energy_lower_bound(
+            rest_mass_density_times_lorentz_factor_ / lorentz_max);
+      }
+      q_ = tau / rest_mass_density_times_lorentz_factor;


(optional) Move this line outside the if as well, and only do q_ = std::max(q_, eps_min) in the EnforcePhysicality branch? Seems easier to understand that way.

nilsvu · 2023-10-17T20:58:17Z

tests/Unit/Evolution/Systems/GrMhd/GhValenciaDivClean/Subcell/Test_PrimsAfterRollback.cpp

+              get<hydro::Tags::RestMassDensity<DataVector>>(subcell_prims));
+    }
+
+    else if constexpr (EquationOfStateType::thermodynamic_dim == 2) {


I recommend you build clang-format into your workflow, it makes editing code a lot easier. Try running git-clang-format and/or use your editor's formatting extensions.

clang-format for some reason doesn't remove this whitespace :(

Maybe the extra empty line tells it not to reformat

isaaclegred · 2023-10-18T23:52:52Z

@nilsvu, @nilsdeppe changes fixed and squashed

nilsvu requested changes Sep 20, 2023

View reviewed changes

src/Evolution/Systems/GrMhd/ValenciaDivClean/KastaunEtAl.cpp Show resolved Hide resolved

src/Evolution/Systems/GrMhd/ValenciaDivClean/PrimitiveFromConservative.cpp Outdated Show resolved Hide resolved

nilsdeppe requested changes Sep 20, 2023

View reviewed changes

nilsdeppe reviewed Sep 20, 2023

View reviewed changes

isaaclegred force-pushed the Add3DEoSSupport branch from 9d14485 to 71c44fd Compare September 27, 2023 17:36

nilsdeppe reviewed Sep 27, 2023

View reviewed changes

nilsvu requested changes Sep 27, 2023

View reviewed changes

isaaclegred force-pushed the Add3DEoSSupport branch 2 times, most recently from 29d4ee5 to b4c78a8 Compare September 28, 2023 21:23

nilsvu requested changes Sep 29, 2023

View reviewed changes

isaaclegred force-pushed the Add3DEoSSupport branch from 6b9c0df to 1132e8b Compare October 6, 2023 15:55

nilsdeppe reviewed Oct 7, 2023

View reviewed changes

isaaclegred force-pushed the Add3DEoSSupport branch from 7124155 to c709337 Compare October 10, 2023 20:56

isaaclegred force-pushed the Add3DEoSSupport branch from c709337 to eaf07b3 Compare October 10, 2023 20:57

nilsdeppe previously approved these changes Oct 13, 2023

View reviewed changes

isaaclegred force-pushed the Add3DEoSSupport branch from 81b3fa9 to 438c6c9 Compare October 13, 2023 16:11

isaaclegred dismissed nilsdeppe’s stale review via e7b4d9a October 13, 2023 22:53

isaaclegred force-pushed the Add3DEoSSupport branch from 438c6c9 to e7b4d9a Compare October 13, 2023 22:53

nilsdeppe previously approved these changes Oct 14, 2023

View reviewed changes

nilsvu previously approved these changes Oct 17, 2023

View reviewed changes

Make GrMhd system compatible with 3D EoS

dfe89fc

isaaclegred dismissed stale reviews from nilsvu and nilsdeppe via dfe89fc October 18, 2023 23:52

isaaclegred force-pushed the Add3DEoSSupport branch from e7b4d9a to dfe89fc Compare October 18, 2023 23:52

nilsdeppe approved these changes Oct 19, 2023

View reviewed changes

nilsvu approved these changes Oct 19, 2023

View reviewed changes

nilsvu merged commit f9cd99e into sxs-collaboration:develop Oct 19, 2023
20 of 22 checks passed

Make GrMhd system compatible with 3D EoS #5479

Make GrMhd system compatible with 3D EoS #5479

Conversation

isaaclegred commented Sep 20, 2023

Proposed changes

Upgrade instructions

Code review checklist

Further comments

isaaclegred commented Sep 20, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ermost commented Sep 21, 2023

isaaclegred commented Sep 25, 2023

isaaclegred commented Sep 25, 2023

nilsdeppe commented Sep 25, 2023 • edited Loading

isaaclegred commented Sep 27, 2023

nilsdeppe left a comment

Choose a reason for hiding this comment

nilsvu commented Sep 27, 2023

nilsvu left a comment

Choose a reason for hiding this comment

ermost commented Oct 3, 2023

isaaclegred commented Oct 4, 2023

ermost commented Oct 4, 2023

isaaclegred commented Oct 6, 2023

nilsdeppe left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ermost commented Oct 10, 2023

isaaclegred commented Oct 10, 2023

nilsdeppe commented Oct 10, 2023

isaaclegred commented Oct 11, 2023

isaaclegred commented Oct 12, 2023

ermost commented Oct 13, 2023

nilsdeppe commented Oct 13, 2023

isaaclegred commented Oct 13, 2023

nilsdeppe left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

nilsvu left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

isaaclegred commented Oct 18, 2023

nilsdeppe commented Sep 25, 2023 •

edited

Loading