Add gas halo to DiskGalaxy problem

[aditivijayan]

Description

The edit reads in initial rho, momr, eint, etot from tables. These tables were generated by the cooling flow code of Stern+19 and modified to include Agora vc and Grackle cooling. The modifications to the original code can be accessed here - https://bitbucket.org/aditivijayan/cooling_flow_stern-19/src/master/
The code generates rho, v, T profiles for Mdot=0.99 Msun/yr and rsonic of 0.38 kpc.

Related issues

Depends on:

• Fix 2D slices of derived vars #1595
• miscellaneous GPU fixes #1596
• add helper functions for SFR diagnostics #1607

Checklist

Before this pull request can be reviewed, all of these tasks should be completed. Denote completed tasks with an x inside the square brackets [ ] in the Markdown source below:

• I have added a description (see above).
• I have added a link to any related issues (if applicable; see above).
• I have read the Contributing Guide.
• I have added tests for any new physics that this PR adds to the code.
• (For quokka-astro org members) I have manually triggered the GPU tests with the magic comment /azp run.

[gemini-code-assist]

Summary of Changes

Hello @aditivijayan, I'm Gemini Code Assist¹! I'm currently reviewing this pull request and will post my feedback shortly. In the meantime, here's a summary to help you and other reviewers quickly get up to speed!

This pull request significantly enhances the AgoraGalaxy problem setup by introducing the capability to model a galactic halo. It achieves this by reading initial density, momentum, and energy profiles from a new input table and integrating these halo properties with the existing disk galaxy model. This provides a more comprehensive and astrophysically relevant initial state for simulations, moving towards a more complete representation of a 1.e12 galaxy.

Highlights

• Halo Property Integration: The simulation now reads and incorporates initial density (rho), radial momentum (momr), internal energy (eint), and total energy (etot) profiles for a galactic halo, enhancing the realism of the initial conditions.
• Input File Update: The configuration for initial conditions has been updated to use a new input file, 'agora_halo_IC.dat', which contains the necessary halo properties, replacing the previous 'vcirc.dat'.
• Enhanced Initial Condition Calculation: The preCalculateInitialConditions function has been extended to parse and store these new halo properties, and the setInitialConditionsOnGrid function now combines these halo properties with the existing disk properties to set the grid's initial state.
• Simplified Velocity Interpolation: The circular velocity interpolation logic has been refined by removing the vcirc_inner and R_table_min parameters, simplifying the boundary conditions for velocity calculations.

Using Gemini Code Assist

The full guide for Gemini Code Assist can be found on our documentation page, here are some quick tips.

Invoking Gemini

You can request assistance from Gemini at any point by creating a comment using either /gemini <command> or @gemini-code-assist <command>. Below is a summary of the supported commands on the current page.

Feature	Command	Description
Code Review	/gemini review	Performs a code review for the current pull request in its current state.
Pull Request Summary	/gemini summary	Provides a summary of the current pull request in its current state.
Comment	@gemini-code-assist	Responds in comments when explicitly tagged, both in pull request comments and review comments.
Help	/gemini help	Displays a list of available commands.

Customization

To customize Gemini Code Assist for GitHub experience, repository maintainers can create a configuration file and/or provide a custom code review style guide (such as PEP-8 for Python) by creating and adding files to a .gemini/ folder in the base of the repository. Detailed instructions can be found here.

Limitations & Feedback

Gemini Code Assist may make mistakes. Please leave feedback on any instances where its feedback is incorrect or counter productive. You can react with 👍 and 👎 on @gemini-code-assist comments. If you're interested in giving your feedback about your experience with Gemini Code Assist for Github and other Google products, sign up here.

You can also get AI-powered code generation, chat, as well as code reviews directly in the IDE at no cost with the Gemini Code Assist IDE Extension.

Footnotes

1. Review the Privacy Notices, Generative AI Prohibited Use Policy, Terms of Service, and learn how to configure Gemini Code Assist in GitHub here. Gemini can make mistakes, so double check it and use code with caution. ↩

[gemini-code-assist]

Code Review

This pull request adds support for reading halo properties from a table. The changes are mostly in testDiskGalaxy.cpp. I've found a critical bug related to variable reuse that could lead to incorrect density calculations. There are also several high-severity issues related to potential out-of-bounds memory access and missing unit conversions. I've also included several medium-severity suggestions to improve code quality, such as removing duplicated code, improving performance of math functions, and adding assertions.

I am having trouble creating individual review comments. Click here to see my feedback.

src/problems/DiskGalaxy/testDiskGalaxy.cpp (396)

There appears to be a critical bug here due to variable shadowing and reuse.

1. The variable rho_halo is calculated from the input table on line 348.
2. This same rho_halo variable is then used in the condition if (rho_halo * T_halo > rho_disk * T_disk) on line 364 to decide if a cell is part of the background halo or the disk.
3. If the condition is true, rho is set to rho_halo on line 365.
4. Finally, on this line, the density is set to rho + rho_halo. If the cell was determined to be halo, this results in rho_halo + rho_halo, effectively doubling the halo density from the table.

The intention seems to be to distinguish between a gaseous disk and a background medium, and then add the new halo component from the table to either of them.

To fix this, you should use distinct variables for the background halo density (calculated from ndens_halo) and the new halo density (from the table).

A suggested fix:

1. Uncomment line 196 and rename the variable to avoid confusion, e.g., rho_ambient.

   // in setInitialConditionsOnGrid, around line 196
   const double rho_ambient = ndens_halo * quokka::EOS_Traits<AgoraGalaxy>::mean_molecular_weight;

2. Use rho_ambient in the if condition on line 364.

   // around line 364
   if (rho_ambient * T_halo > rho_disk * T_disk) {
       rho = rho_ambient;
       // ...
   }

3. The final state update on lines 396-401 will then correctly superimpose the halo from the table onto either the disk or the ambient medium.

src/problems/DiskGalaxy/testDiskGalaxy.cpp (135-138)

The code converts radius_h and vcirc_h from problem-specific units (kpc, km/s) to CGS. However, the new quantities rho_h, momr_h, eint_h, and etot_h are assigned directly without unit conversion. If these values are not already in CGS units in the input file, this will lead to incorrect physical results. Please add comments to clarify the units of the input data and apply conversions if necessary for correctness and clarity.

src/problems/DiskGalaxy/testDiskGalaxy.cpp (238-245)

The logic for handling radii smaller than the minimum radius in the table (R_table_min) has been removed. The previous implementation used a constant vcirc_inner for R < R_table_min. The new implementation will call interpolate_value for this case, which will cause an assertion failure if R is less than the first radius value in the table. This is a potential regression.

Please restore robust handling for small radii. A simple way to do this is to clamp to the first value in the table. You can achieve this by using interpolate_value<quokka::BoundaryPolicy::Clamp>.

auto vcirc_exact = [R_table_max, R_table, vcirc_outer, vcirc_table, len_table](const amrex::Real R) {
    double vcirc = NAN;
    if ((R <= R_table_max)) { // also changed to <= for consistency
        vcirc = interpolate_value<quokka::BoundaryPolicy::Clamp>(R, R_table, vcirc_table, len_table);
    } else {
        vcirc = vcirc_outer;
    }
    return vcirc;
};

src/problems/DiskGalaxy/testDiskGalaxy.cpp (125)

This line is a duplicate of the one above and can be removed.

src/problems/DiskGalaxy/testDiskGalaxy.cpp (261-300)

These four new lambda functions (rhoHalo, momHalo, eintHalo, etotHalo) are nearly identical. This code duplication can be avoided by creating a helper function or a generic lambda factory. This would also be a good place to fix the missing handling for radii smaller than the table's minimum radius, as mentioned in another comment.

For example, you could create a generic lambda factory that also handles out-of-bounds cases by clamping to the table boundaries:

auto make_halo_interpolator = [R_table, R_table_max, len_table](auto const *table, auto outer_val) {
    return [=](const amrex::Real R) {
        double value = NAN;
        if (R <= R_table_max) {
            // Using BoundaryPolicy::Clamp handles R < R_table_min by clamping to the first value.
            value = interpolate_value<quokka::BoundaryPolicy::Clamp>(R, R_table, table, len_table);
        } else {
            value = outer_val;
        }
        return value;
    };
};

auto rhoHalo = make_halo_interpolator(rhoH_table, rho_outer);
auto momHalo = make_halo_interpolator(momr_table, momr_outer);
auto eintHalo = make_halo_interpolator(eint_table, eint_outer);
auto etotHalo = make_halo_interpolator(etot_table, etot_outer);

This refactoring would make the code more concise and robust.

src/problems/DiskGalaxy/testDiskGalaxy.cpp (303-342)

There are a few small improvements that can be made in this block of lambda functions:

• std::pow(x, 2) is used for squaring. It's generally more efficient to use x*x.
• For calculating Euclidean norms, std::hypot is preferred over std::sqrt(x*x + ...) as it provides better numerical stability by avoiding intermediate overflow or underflow. For example, std::sqrt(std::pow(x, 2) + std::pow(y, 2) + std::pow(z, 2)) can be replaced with std::hypot(x, y, z).
• There is a stray semicolon ;; on lines 323, 332, and 340.

src/problems/DiskGalaxy/testDiskGalaxy.cpp (348-353)

It's good practice to add assertions to check for NaN values for the newly computed halo quantities (rho_halo, momx_halo, etc.), similar to the existing assertion for rho_disk. This will help catch potential issues early during development and debugging.

[aditivijayan]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[BenWibking]

The cooling rates in the Python code are not quite right. It's necessary to add the heating rates to the cooling rates to get a net cooling rate, and also add PE heating and Compton cooling by hand.

You can copy the code from here to do that:

quokka/extern/cooling/grackle_tables.py

Line 88 in 5fbfea8

def cooling_rate(nH, T, zmet, redshift=0., tables=None, include_pe=True):

[aditivijayan]

The cooling rates in the Python code are not quite right. It's necessary to add the heating rates to the cooling rates to get a net cooling rate, and also add PE heating and Compton cooling by hand.

You can copy the code from here to do that:

quokka/extern/cooling/grackle_tables.py

Line 88 in 5fbfea8

def cooling_rate(nH, T, zmet, redshift=0., tables=None, include_pe=True):

The net cooling rate is negative for some combination of T and nH. How do you propose I estimate dlnLambda_dlnrho and dlnLambda_dlnT? Having negative values of lambda just throws up nan when estimating the gradients.

[markkrumholz]

Physically the net cooling rate can indeed be negative for some combinations of n and T, because those correspond to places where the gas is colder than the equilibrium temperature given the radiation field to which it is being subjected. But does the solution actually pass through any of those points? I would guess not. If you need to supply a numerical value there anyway just to avoid numerical or computer issues, I would suggest that you instead just supply d ln(|Lambda|) / d ln T instead of d ln Lambda / d ln T, i.e., just take the absolute value of the cooling rate.

[BenWibking]

You can't take the logarithm when it's negative, but you should be able to rewrite it with the chain rule (by physicist math):
$\frac{d ln \Lambda}{d ln T} =\frac{T}{\Lambda} \frac{d \Lambda}{d T}$

Mathematically, this should be the analytic continuation of the function when $\Lambda < 0$.

[aditivijayan]

I've pushed the changes to bitbucket. Now profiles include PE and Compton. This is what they look like.

[aditivijayan]

I've pushed the changes to bitbucket. Now profiles include PE and Compton. This is what they look like.

Profile regenerated for Z=Zsun to match with the metallicity of the cooling table.

[BenWibking]

I don't think it makes sense to continue the solution inward of the sonic point. Can you remove that part of the solution?

[BenWibking]

It looks like it's going to hit another sonic point just outside of R200. If you want me to try to debug this, I have some time tomorrow.

[markkrumholz]

Stern et al. do in fact continue their solution inward of the sonic point, so I don't think it is necessarily a problem to do so -- and I think we do want to do so, because otherwise we are going to have an artificial density jump in our initial condition. There's no mass in this gas in any event, so I see no harm in continuing the solution inward. I think the second sonic point is just an artifact of the cooling function we're using, which has a temperature that starts to drop after some density rather than rising continuously as Stern et al. assume. Again, I don't think this is necessarily a problem or an error.

[aditivijayan]

I don't think it makes sense to continue the solution inward of the sonic point. Can you remove that part of the solution?

And just copy the solution from r=rsonic? Sure.

[aditivijayan]

Stern et al. do in fact continue their solution inward of the sonic point, so I don't think it is necessarily a problem to do so -- and I think we do want to do so, because otherwise we are going to have an artificial density jump in our initial condition. There's no mass in this gas in any event, so I see no harm in continuing the solution inward. I think the second sonic point is just an artifact of the cooling function we're using, which has a temperature that starts to drop after some density rather than rising continuously as Stern et al. assume. Again, I don't think this is necessarily a problem or an error.

This is correct.

[BenWibking]

I don't think it makes sense to continue the solution inward of the sonic point. Can you remove that part of the solution?

And just copy the solution from r=rsonic? Sure.

This would be my preferred solution. It remains continuous in this case (although not differentiable).

I think it's not any less meaningful than the continued solution, since the smoothly-continued solution won't be realized in practice (there will be a shock, it will become time-dependent, etc.).

[BenWibking]

Stern et al. do in fact continue their solution inward of the sonic point, so I don't think it is necessarily a problem to do so -- and I think we do want to do so, because otherwise we are going to have an artificial density jump in our initial condition. There's no mass in this gas in any event, so I see no harm in continuing the solution inward. I think the second sonic point is just an artifact of the cooling function we're using, which has a temperature that starts to drop after some density rather than rising continuously as Stern et al. assume. Again, I don't think this is necessarily a problem or an error.

My mistake. That makes sense.

In terms of what physical solution to choose, I think it might be preferable to choose the sonic point location (or mass accretion rate) to maximize observational agreement. Essentially all the cooling flow solutions give temperatures ~1e6 K, so there's nothing to calibrate there. For density, the observational constraints are weak, but I think one that is as good as any other is the Miller & Bregman (https://ui.adsabs.harvard.edu/abs/2013ApJ...770..118M/abstract) electron density profile, which goes at $n_e^{-1.7}$ (see their paper for the normalization). This is consistent with all the more recent measurements I'm aware of within the (very large) observational uncertainties.

[aditivijayan]

Stern et al. do in fact continue their solution inward of the sonic point, so I don't think it is necessarily a problem to do so -- and I think we do want to do so, because otherwise we are going to have an artificial density jump in our initial condition. There's no mass in this gas in any event, so I see no harm in continuing the solution inward. I think the second sonic point is just an artifact of the cooling function we're using, which has a temperature that starts to drop after some density rather than rising continuously as Stern et al. assume. Again, I don't think this is necessarily a problem or an error.

My mistake. That makes sense.

In terms of what physical solution to choose, I think it might be preferable to choose the sonic point location (or mass accretion rate) to maximize observational agreement. Essentially all the cooling flow solutions give temperatures ~1e6 K, so there's nothing to calibrate there. For density, the observational constraints are weak, but I think one that is as good as any other is the Miller & Bregman (https://ui.adsabs.harvard.edu/abs/2013ApJ...770..118M/abstract) electron density profile, which goes at n e − 1.7 (see their paper for the normalization). This is consistent with all the more recent measurements I'm aware of within the (very large) observational uncertainties.

The normalisation I am getting for n_e (0.024) is within the error range (0.048(+0.085/-0.037). Their accretion rate is Msun/yr for large r.

[BenWibking]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[BenWibking]

@markkrumholz Can you review this? I've made some changes, so it would be good to have a third person look at it all.

[markkrumholz]

This generally looks fine except that that changes in the infrastructure to compute the SFR in PhysicsParticles seem to have been duplicated in this PR. Once they are removed, I'm happy to approve.

[BenWibking]

This generally looks fine except that that changes in the infrastructure to compute the SFR in PhysicsParticles seem to have been duplicated in this PR. Once they are removed, I'm happy to approve.

@markkrumholz I've removed those changes.

[sonarqubecloud]

Quality Gate passed

Issues
3 New issues
0 Accepted issues

Measures
0 Security Hotspots
0.0% Coverage on New Code
9.5% Duplication on New Code

See analysis details on SonarQube Cloud

[BenWibking]

/azp run

[azure-pipelines]

Azure Pipelines successfully started running 2 pipeline(s).

[BenWibking]

@markkrumholz This PR is ready to be approved.