Skip to content

Conversation

@bendudson
Copy link
Collaborator

A relatively simple fix that seems to greatly help reduce unphysical behavior near density floors.

Defines the internal energy E as E = Cv * Nlim * T rather than the usual Cv * N * T. This effectively puts a floor on the heat capacity of the fluid. Then solve the pressure equation consistent with this equation of state.

Changes made to pressure equations in evolve_pressure, neutral_mixed and neutral_full_velocity. Unfortunately I ran clang-format on the files so the changes seem numerous. The actual change is quite small: In the pressure advection replace (5/3)P with Nlim * T + (2/3) * P.

This is an example from the midplane of a DIII-D simulation using neutral_full_velocity. The density floor is 1e-5, and previously produced quite unphysical behavior.

image image

bendudson added 3 commits June 8, 2025 20:50
Consistently treat internal energy as Nlim * T, so that
advection of thermal energy takes into account the density floor.

Appears to reduce artefacts where the density falls below the floor.
Improve handling of low density regions by defining the fluid
internal energy as Nlim * T. Modifies the pressure equation.
The internal energy of a fluid is modified to be Nlim * T.  This
changes the pressure equation to be consistent with the application of
a density floor.
@mikekryjak
Copy link
Collaborator

@bendudson this is awesome. Looking at evolve_pressure only for now: am I understanding correctly that this in effect makes the density floor in use everywhere consistently, whereas before it was only used where we divided by zero? And P_solver is the pressure without floors that the solver sees, while P is the pressure calculated using floored density?

What about flooring pressure itself? I see that now it's only flooring the density.

I raised an issue to put clang-format in the CI: #360

@bendudson
Copy link
Collaborator Author

Hey @mikekryjak . Thanks! The solver now evolves a modified internal energy P_solver = N_lim T, which is divided by N_lim to get temperature, and then the pressure is P = N T. This is why neutral pressure can continue to fall into the core, even when temperature is increasing and the density has hit the floor: the internal energy P_solver is going up into the core, while pressure P is going down.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants