Skip to content
Snippets Groups Projects
Commit 6f78a621 authored by Jacob Kegerreis's avatar Jacob Kegerreis Committed by Matthieu Schaller
Browse files

Revert "Store the gravity acceleration in the xpart so that we don't access p->gpart in drift_part"

This reverts commit 5e1ae7ca.
parent 8febf700
Branches
Tags
No related merge requests found
...@@ -425,6 +425,18 @@ elif test "$fixed_boundary_particles" != "no"; then ...@@ -425,6 +425,18 @@ elif test "$fixed_boundary_particles" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_FIXED_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.]) AC_DEFINE_UNQUOTED([SWIFT_FIXED_BOUNDARY_PARTICLES], [$enableval] ,[Particles with smaller ID than this will be considered as boundaries.])
fi fi
# Check if fixed entropy is on for settling planetary initial conditions
AC_ARG_ENABLE([planetary-fixed-entropy],
[AS_HELP_STRING([--enable-planetary-fixed-entropy],
[Force entropies to stay fixed for settling planetary initial conditions @<:@yes/no@:>@]
)],
[planetary_fixed_entropy="$enableval"],
[planetary_fixed_entropy="no"]
)
if test "$planetary_fixed_entropy" = "yes"; then
AC_DEFINE([PLANETARY_FIXED_ENTROPY],1,[Enable planetary fixed entropy])
fi
# Check whether we have any of the ARM v8.1 tick timers # Check whether we have any of the ARM v8.1 tick timers
AX_ASM_ARM_PMCCNTR AX_ASM_ARM_PMCCNTR
AX_ASM_ARM_CNTVCT AX_ASM_ARM_CNTVCT
...@@ -2468,6 +2480,7 @@ AC_MSG_RESULT([ ...@@ -2468,6 +2480,7 @@ AC_MSG_RESULT([
Custom icbrtf : $enable_custom_icbrtf Custom icbrtf : $enable_custom_icbrtf
Boundary particles : $boundary_particles Boundary particles : $boundary_particles
Fixed boundary particles : $fixed_boundary_particles Fixed boundary particles : $fixed_boundary_particles
Planetary fixed entropy : $planetary_fixed_entropy
Particle Logger : $with_logger Particle Logger : $with_logger
Python enabled : $have_python Python enabled : $have_python
......
...@@ -26,7 +26,7 @@ These allow for multiple materials to be used, ...@@ -26,7 +26,7 @@ These allow for multiple materials to be used,
chosen from the several available equations of state. chosen from the several available equations of state.
.. toctree:: .. toctree::
:maxdepth: 2 :maxdepth: 1
:caption: More information: :caption: More information:
Hydro Scheme <hydro_scheme> Hydro Scheme <hydro_scheme>
......
...@@ -20,4 +20,19 @@ Ruiz-Bonilla et al. (2020). ...@@ -20,4 +20,19 @@ Ruiz-Bonilla et al. (2020).
They are available with documentation and examples at They are available with documentation and examples at
https://github.com/srbonilla/WoMa and https://github.com/jkeger/seagen, https://github.com/srbonilla/WoMa and https://github.com/jkeger/seagen,
or can be installed directly with ``pip`` or can be installed directly with ``pip``
(https://pypi.org/project/woma/, https://pypi.org/project/seagen/). (https://pypi.org/project/woma/, https://pypi.org/project/seagen/).
\ No newline at end of file
Settling initial conditions with fixed entropies
------------------------------------------------
If the particles' equations of state include specific entropies,
and the initial conditions file includes specific entropies for each particle
(in ``PartType0/Entropies``),
then configuring SWIFT with ``--enable-planetary-fixed-entropy``
will override the internal energy of each particle each step such that its
specific entropy remains constant.
This should be used with caution, but may be a convenient way to maintain an
entropy profile while initial conditions settle to equilibrium with their
slightly different SPH densities.
...@@ -47,9 +47,9 @@ struct SESAME_params { ...@@ -47,9 +47,9 @@ struct SESAME_params {
float *table_log_u_rho_T; float *table_log_u_rho_T;
float *table_P_rho_T; float *table_P_rho_T;
float *table_c_rho_T; float *table_c_rho_T;
float *table_s_rho_T; float *table_log_s_rho_T;
int date, num_rho, num_T; int date, num_rho, num_T;
float P_tiny, c_tiny; float u_tiny, P_tiny, c_tiny, s_tiny;
enum eos_planetary_material_id mat_id; enum eos_planetary_material_id mat_id;
}; };
...@@ -138,7 +138,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat, ...@@ -138,7 +138,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
(float *)malloc(mat->num_rho * mat->num_T * sizeof(float)); (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
mat->table_c_rho_T = mat->table_c_rho_T =
(float *)malloc(mat->num_rho * mat->num_T * sizeof(float)); (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
mat->table_s_rho_T = mat->table_log_s_rho_T =
(float *)malloc(mat->num_rho * mat->num_T * sizeof(float)); (float *)malloc(mat->num_rho * mat->num_T * sizeof(float));
// Densities (not log yet) // Densities (not log yet)
...@@ -159,7 +159,8 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat, ...@@ -159,7 +159,8 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
if (c != 1) error("Failed to read the SESAME EoS table %s", table_file); if (c != 1) error("Failed to read the SESAME EoS table %s", table_file);
} }
// Sp. int. energies (not log yet), pressures, sound speeds, and entropies // Sp. int. energies (not log yet), pressures, sound speeds, and sp.
// entropies (not log yet)
for (int i_T = -1; i_T < mat->num_T; i_T++) { for (int i_T = -1; i_T < mat->num_T; i_T++) {
for (int i_rho = -1; i_rho < mat->num_rho; i_rho++) { for (int i_rho = -1; i_rho < mat->num_rho; i_rho++) {
// Ignore the first elements of rho = 0, T = 0 // Ignore the first elements of rho = 0, T = 0
...@@ -171,7 +172,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat, ...@@ -171,7 +172,7 @@ INLINE static void load_table_SESAME(struct SESAME_params *mat,
&mat->table_log_u_rho_T[i_rho * mat->num_T + i_T], &mat->table_log_u_rho_T[i_rho * mat->num_T + i_T],
&mat->table_P_rho_T[i_rho * mat->num_T + i_T], &mat->table_P_rho_T[i_rho * mat->num_T + i_T],
&mat->table_c_rho_T[i_rho * mat->num_T + i_T], &mat->table_c_rho_T[i_rho * mat->num_T + i_T],
&mat->table_s_rho_T[i_rho * mat->num_T + i_T]); &mat->table_log_s_rho_T[i_rho * mat->num_T + i_T]);
if (c != 4) error("Failed to read the SESAME EoS table %s", table_file); if (c != 4) error("Failed to read the SESAME EoS table %s", table_file);
} }
} }
...@@ -188,22 +189,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) { ...@@ -188,22 +189,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
mat->table_log_rho[i_rho] = logf(mat->table_log_rho[i_rho]); mat->table_log_rho[i_rho] = logf(mat->table_log_rho[i_rho]);
} }
// Convert sp. int. energies to log(sp. int. energy) // Initialise tiny values
for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) { mat->u_tiny = FLT_MAX;
for (int i_T = 0; i_T < mat->num_T; i_T++) {
// If not positive then set very small for the log
if (mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] <= 0) {
mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] = 1.f;
}
mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] =
logf(mat->table_log_u_rho_T[i_rho * mat->num_T + i_T]);
}
}
// Initialise tiny pressure and soundspeed
mat->P_tiny = FLT_MAX; mat->P_tiny = FLT_MAX;
mat->c_tiny = FLT_MAX; mat->c_tiny = FLT_MAX;
mat->s_tiny = FLT_MAX;
// Enforce that the 1D arrays of u (at each rho) are monotonic // Enforce that the 1D arrays of u (at each rho) are monotonic
// This is necessary because, for some high-density u slices at very low T, // This is necessary because, for some high-density u slices at very low T,
...@@ -223,7 +213,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) { ...@@ -223,7 +213,11 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
break; break;
} }
// Smallest positive pressure and sound speed // Smallest positive values
if ((mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] < mat->u_tiny) &&
(mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] > 0)) {
mat->u_tiny = mat->table_log_u_rho_T[i_rho * mat->num_T + i_T];
}
if ((mat->table_P_rho_T[i_rho * mat->num_T + i_T] < mat->P_tiny) && if ((mat->table_P_rho_T[i_rho * mat->num_T + i_T] < mat->P_tiny) &&
(mat->table_P_rho_T[i_rho * mat->num_T + i_T] > 0)) { (mat->table_P_rho_T[i_rho * mat->num_T + i_T] > 0)) {
mat->P_tiny = mat->table_P_rho_T[i_rho * mat->num_T + i_T]; mat->P_tiny = mat->table_P_rho_T[i_rho * mat->num_T + i_T];
...@@ -232,12 +226,38 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) { ...@@ -232,12 +226,38 @@ INLINE static void prepare_table_SESAME(struct SESAME_params *mat) {
(mat->table_c_rho_T[i_rho * mat->num_T + i_T] > 0)) { (mat->table_c_rho_T[i_rho * mat->num_T + i_T] > 0)) {
mat->c_tiny = mat->table_c_rho_T[i_rho * mat->num_T + i_T]; mat->c_tiny = mat->table_c_rho_T[i_rho * mat->num_T + i_T];
} }
if ((mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] < mat->s_tiny) &&
(mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] > 0)) {
mat->s_tiny = mat->table_log_s_rho_T[i_rho * mat->num_T + i_T];
}
} }
} }
// Tiny pressure to allow interpolation near non-positive values // Tiny values to allow interpolation near non-positive values
mat->u_tiny *= 1e-3f;
mat->P_tiny *= 1e-3f; mat->P_tiny *= 1e-3f;
mat->c_tiny *= 1e-3f; mat->c_tiny *= 1e-3f;
mat->s_tiny *= 1e-3f;
// Convert sp. int. energies to log(sp. int. energy), same for sp. entropies
for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
for (int i_T = 0; i_T < mat->num_T; i_T++) {
// If not positive then set very small for the log
if (mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] <= 0) {
mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] = mat->u_tiny;
}
mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] =
logf(mat->table_log_u_rho_T[i_rho * mat->num_T + i_T]);
if (mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] <= 0) {
mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] = mat->s_tiny;
}
mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] =
logf(mat->table_log_s_rho_T[i_rho * mat->num_T + i_T]);
}
}
} }
// Convert to internal units // Convert to internal units
...@@ -255,7 +275,7 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat, ...@@ -255,7 +275,7 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
units_cgs_conversion_factor(us, UNIT_CONV_DENSITY)); units_cgs_conversion_factor(us, UNIT_CONV_DENSITY));
} }
// Sp. Int. Energies (log), pressures, and sound speeds // Sp. int. energies (log), pressures, sound speeds, and sp. entropies
for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) { for (int i_rho = 0; i_rho < mat->num_rho; i_rho++) {
for (int i_T = 0; i_T < mat->num_T; i_T++) { for (int i_T = 0; i_T < mat->num_T; i_T++) {
mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] += logf( mat->table_log_u_rho_T[i_rho * mat->num_T + i_T] += logf(
...@@ -267,26 +287,118 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat, ...@@ -267,26 +287,118 @@ INLINE static void convert_units_SESAME(struct SESAME_params *mat,
mat->table_c_rho_T[i_rho * mat->num_T + i_T] *= mat->table_c_rho_T[i_rho * mat->num_T + i_T] *=
1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) / 1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
units_cgs_conversion_factor(us, UNIT_CONV_SPEED); units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
mat->table_s_rho_T[i_rho * mat->num_T + i_T] *= mat->table_log_s_rho_T[i_rho * mat->num_T + i_T] +=
units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) / logf(units_cgs_conversion_factor(
units_cgs_conversion_factor(us, UNIT_CONV_ENTROPY); &si, UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS) /
units_cgs_conversion_factor(
us, UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS));
} }
} }
// Tiny pressure and sound speed // Tiny values
mat->u_tiny *=
units_cgs_conversion_factor(&si, UNIT_CONV_ENERGY_PER_UNIT_MASS) /
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
mat->P_tiny *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) / mat->P_tiny *= units_cgs_conversion_factor(&si, UNIT_CONV_PRESSURE) /
units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE); units_cgs_conversion_factor(us, UNIT_CONV_PRESSURE);
mat->c_tiny *= 1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) / mat->c_tiny *= 1e3f * units_cgs_conversion_factor(&si, UNIT_CONV_SPEED) /
units_cgs_conversion_factor(us, UNIT_CONV_SPEED); units_cgs_conversion_factor(us, UNIT_CONV_SPEED);
mat->s_tiny *=
units_cgs_conversion_factor(&si,
UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS) /
units_cgs_conversion_factor(us, UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS);
} }
// gas_internal_energy_from_entropy // gas_internal_energy_from_entropy
INLINE static float SESAME_internal_energy_from_entropy( INLINE static float SESAME_internal_energy_from_entropy(
float density, float entropy, const struct SESAME_params *mat) { float density, float entropy, const struct SESAME_params *mat) {
error("This EOS function is not yet implemented!"); float u, log_u_1, log_u_2, log_u_3, log_u_4;
return 0.f; if (entropy <= 0.f) {
return 0.f;
}
int idx_rho, idx_s_1, idx_s_2;
float intp_rho, intp_s_1, intp_s_2;
const float log_rho = logf(density);
const float log_s = logf(entropy);
// 2D interpolation (bilinear with log(rho), log(s)) to find u(rho, s))
// Density index
idx_rho =
find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
// Sp. entropy at this and the next density (in relevant slice of s array)
idx_s_1 = find_value_in_monot_incr_array(
log_s, mat->table_log_s_rho_T + idx_rho * mat->num_T, mat->num_T);
idx_s_2 = find_value_in_monot_incr_array(
log_s, mat->table_log_s_rho_T + (idx_rho + 1) * mat->num_T, mat->num_T);
// If outside the table then extrapolate from the edge and edge-but-one values
if (idx_rho <= -1) {
idx_rho = 0;
} else if (idx_rho >= mat->num_rho) {
idx_rho = mat->num_rho - 2;
}
if (idx_s_1 <= -1) {
idx_s_1 = 0;
} else if (idx_s_1 >= mat->num_T) {
idx_s_1 = mat->num_T - 2;
}
if (idx_s_2 <= -1) {
idx_s_2 = 0;
} else if (idx_s_2 >= mat->num_T) {
idx_s_2 = mat->num_T - 2;
}
// Check for duplicates in SESAME tables before interpolation
if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
(mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
} else {
intp_rho = 1.f;
}
if (mat->table_log_s_rho_T[idx_rho * mat->num_T + (idx_s_1 + 1)] !=
mat->table_log_s_rho_T[idx_rho * mat->num_T + idx_s_1]) {
intp_s_1 =
(log_s - mat->table_log_s_rho_T[idx_rho * mat->num_T + idx_s_1]) /
(mat->table_log_s_rho_T[idx_rho * mat->num_T + (idx_s_1 + 1)] -
mat->table_log_s_rho_T[idx_rho * mat->num_T + idx_s_1]);
} else {
intp_s_1 = 1.f;
}
if (mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + (idx_s_2 + 1)] !=
mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2]) {
intp_s_2 =
(log_s - mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2]) /
(mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + (idx_s_2 + 1)] -
mat->table_log_s_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2]);
} else {
intp_s_2 = 1.f;
}
// Table values
log_u_1 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1];
log_u_2 = mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_s_1 + 1];
log_u_3 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2];
log_u_4 = mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_s_2 + 1];
// If below the minimum s at this rho then just use the lowest table values
if ((idx_rho > 0.f) && ((intp_s_1 < 0.f) || (intp_s_2 < 0.f) ||
(log_u_1 > log_u_2) || (log_u_3 > log_u_4))) {
intp_s_1 = 0;
intp_s_2 = 0;
}
// Interpolate with the log values
u = (1.f - intp_rho) * ((1.f - intp_s_1) * log_u_1 + intp_s_1 * log_u_2) +
intp_rho * ((1.f - intp_s_2) * log_u_3 + intp_s_2 * log_u_4);
// Convert back from log
u = expf(u);
return u;
} }
// gas_pressure_from_entropy // gas_pressure_from_entropy
...@@ -338,7 +450,7 @@ INLINE static float SESAME_pressure_from_internal_energy( ...@@ -338,7 +450,7 @@ INLINE static float SESAME_pressure_from_internal_energy(
const float log_rho = logf(density); const float log_rho = logf(density);
const float log_u = logf(u); const float log_u = logf(u);
// 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u) // 2D interpolation (bilinear with log(rho), log(u)) to find P(rho, u))
// Density index // Density index
idx_rho = idx_rho =
find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho); find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
...@@ -463,7 +575,7 @@ INLINE static float SESAME_soundspeed_from_internal_energy( ...@@ -463,7 +575,7 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
const float log_rho = logf(density); const float log_rho = logf(density);
const float log_u = logf(u); const float log_u = logf(u);
// 2D interpolation (bilinear with log(rho), log(u)) to find c(rho, u) // 2D interpolation (bilinear with log(rho), log(u)) to find c(rho, u))
// Density index // Density index
idx_rho = idx_rho =
find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho); find_value_in_monot_incr_array(log_rho, mat->table_log_rho, mat->num_rho);
......
...@@ -25,7 +25,7 @@ ...@@ -25,7 +25,7 @@
* equations) * equations)
* *
* The thermal variable is the internal energy (u). Simple constant * The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction * viscosity term with the Balsara (1995) switch. No thermal conduction
* term is implemented. * term is implemented.
* *
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
......
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
* @brief Minimal conservative implementation of SPH (Debugging routines) * @brief Minimal conservative implementation of SPH (Debugging routines)
* *
* The thermal variable is the internal energy (u). Simple constant * The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction * viscosity term with the Balsara (1995) switch. No thermal conduction
* term is implemented. * term is implemented.
* *
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
......
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
* @brief Minimal conservative implementation of SPH (Neighbour loop equations) * @brief Minimal conservative implementation of SPH (Neighbour loop equations)
* *
* The thermal variable is the internal energy (u). Simple constant * The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction * viscosity term with the Balsara (1995) switch. No thermal conduction
* term is implemented. * term is implemented.
* *
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
......
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
* @brief Minimal conservative implementation of SPH (i/o routines) * @brief Minimal conservative implementation of SPH (i/o routines)
* *
* The thermal variable is the internal energy (u). Simple constant * The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction * viscosity term with the Balsara (1995) switch. No thermal conduction
* term is implemented. * term is implemented.
* *
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
...@@ -225,8 +225,9 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) { ...@@ -225,8 +225,9 @@ INLINE static void hydro_write_flavour(hid_t h_grpsph) {
/* Viscosity and thermal conduction */ /* Viscosity and thermal conduction */
/* Nothing in this minimal model... */ /* Nothing in this minimal model... */
io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment"); io_write_attribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
io_write_attribute_s(h_grpsph, "Viscosity Model", io_write_attribute_s(
"Minimal treatment as in Monaghan (1992)"); h_grpsph, "Viscosity Model",
"as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
} }
/** /**
......
...@@ -24,7 +24,7 @@ ...@@ -24,7 +24,7 @@
* @brief Minimal conservative implementation of SPH (Particle definition) * @brief Minimal conservative implementation of SPH (Particle definition)
* *
* The thermal variable is the internal energy (u). Simple constant * The thermal variable is the internal energy (u). Simple constant
* viscosity term without switches is implemented. No thermal conduction * viscosity term with the Balsara (1995) switch. No thermal conduction
* term is implemented. * term is implemented.
* *
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with * This corresponds to equations (43), (44), (45), (101), (103) and (104) with
......
...@@ -392,10 +392,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p, ...@@ -392,10 +392,15 @@ hydro_set_drifted_physical_internal_energy(struct part *p,
/* Now recompute the extra quantities */ /* Now recompute the extra quantities */
/* Compute the sound speed */ /* Compute the sound speed */
const float pressure =
gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
const float soundspeed = hydro_get_comoving_soundspeed(p); const float soundspeed = hydro_get_comoving_soundspeed(p);
/* Update variables. */ /* Update variables. */
p->force.pressure = pressure;
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
} }
/** /**
...@@ -472,6 +477,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part( ...@@ -472,6 +477,10 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount_dh = 0.f; p->density.wcount_dh = 0.f;
p->rho = 0.f; p->rho = 0.f;
p->density.rho_dh = 0.f; p->density.rho_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
} }
/** /**
...@@ -507,6 +516,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( ...@@ -507,6 +516,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
p->density.rho_dh *= h_inv_dim_plus_one; p->density.rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= h_inv_dim; p->density.wcount *= h_inv_dim;
p->density.wcount_dh *= h_inv_dim_plus_one; p->density.wcount_dh *= h_inv_dim_plus_one;
const float rho_inv = 1.f / p->rho;
const float a_inv2 = cosmo->a2_inv;
/* Finish calculation of the (physical) velocity curl components */
p->density.rot_v[0] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
p->density.rot_v[1] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
p->density.rot_v[2] *= h_inv_dim_plus_one * a_inv2 * rho_inv;
/* Finish calculation of the (physical) velocity divergence */
p->density.div_v *= h_inv_dim_plus_one * a_inv2 * rho_inv;
} }
/** /**
...@@ -534,6 +554,10 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours( ...@@ -534,6 +554,10 @@ __attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
p->density.wcount = kernel_root * h_inv_dim; p->density.wcount = kernel_root * h_inv_dim;
p->density.rho_dh = 0.f; p->density.rho_dh = 0.f;
p->density.wcount_dh = 0.f; p->density.wcount_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
p->density.rot_v[2] = 0.f;
} }
/** /**
...@@ -558,15 +582,22 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -558,15 +582,22 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const struct cosmology *cosmo, const struct hydro_props *hydro_props, const struct cosmology *cosmo, const struct hydro_props *hydro_props,
const float dt_alpha) { const float dt_alpha) {
const float fac_mu = cosmo->a_factor_mu; const float fac_Balsara_eps = cosmo->a_factor_Balsara_eps;
/* Compute the norm of the curl */ /* Compute the norm of the curl */
const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] + const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
p->density.rot_v[1] * p->density.rot_v[1] + p->density.rot_v[1] * p->density.rot_v[1] +
p->density.rot_v[2] * p->density.rot_v[2]); p->density.rot_v[2] * p->density.rot_v[2]);
/* Compute the norm of div v */ /* Compute the norm of div v including the Hubble flow term */
const float abs_div_v = fabsf(p->density.div_v); const float div_physical_v = p->density.div_v + hydro_dimension * cosmo->H;
const float abs_div_physical_v = fabsf(div_physical_v);
#ifdef PLANETARY_FIXED_ENTROPY
/* Override the internal energy to satisfy the fixed entropy */
p->u = gas_internal_energy_from_entropy(p->rho, p->s_fixed, p->mat_id);
xp->u_full = p->u;
#endif
/* Compute the pressure */ /* Compute the pressure */
const float pressure = const float pressure =
...@@ -576,23 +607,30 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( ...@@ -576,23 +607,30 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float soundspeed = const float soundspeed =
gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id); gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id);
/* Compute the "grad h" term */ /* Compute the "grad h" term - Note here that we have \tilde{x}
const float rho_inv = 1.f / p->rho; * as 1 as we use the local number density to find neighbours. This
float rho_dh = p->density.rho_dh; * introduces a j-component that is considered in the force loop,
* meaning that this cached grad_h_term gives:
*
* f_ij = 1.f - grad_h_term_i / m_j */
const float common_factor = p->h * hydro_dimension_inv / p->density.wcount;
float grad_h_term = common_factor * p->density.rho_dh /
(1.f + common_factor * p->density.wcount_dh);
/* Ignore changing-kernel effects when h ~= h_max */ /* Ignore changing-kernel effects when h ~= h_max */
if (p->h > 0.9999f * hydro_props->h_max) { if (p->h > 0.9999f * hydro_props->h_max) {
rho_dh = 0.f; grad_h_term = 0.f;
} }
const float grad_h_term =
1.f / (1.f + hydro_dimension_inv * p->h * rho_dh * rho_inv);
/* Compute the Balsara switch */ /* Compute the Balsara switch */
#ifdef PLANETARY_SPH_NO_BALSARA #ifdef PLANETARY_SPH_NO_BALSARA
const float balsara = hydro_props->viscosity.alpha; const float balsara = hydro_props->viscosity.alpha;
#else #else
const float balsara = /* Pre-multiply in the AV factor; hydro_props are not passed to the iact
hydro_props->viscosity.alpha * abs_div_v / * functions */
(abs_div_v + curl_v + 0.0001f * fac_mu * soundspeed / p->h); const float balsara = hydro_props->viscosity.alpha * abs_div_physical_v /
(abs_div_physical_v + curl_v +
0.0001f * fac_Balsara_eps * soundspeed / p->h);
#endif #endif
/* Update variables. */ /* Update variables. */
...@@ -641,12 +679,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( ...@@ -641,12 +679,12 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->v[1] = xp->v_full[1]; p->v[1] = xp->v_full[1];
p->v[2] = xp->v_full[2]; p->v[2] = xp->v_full[2];
/* Re-set the entropy */ /* Re-set the internal energy */
p->u = xp->u_full; p->u = xp->u_full;
/* Compute the pressure */ /* Compute the pressure */
const float pressure = const float pressure =
gas_pressure_from_internal_energy(p->rho, xp->u_full, p->mat_id); gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id);
/* Compute the sound speed */ /* Compute the sound speed */
const float soundspeed = const float soundspeed =
...@@ -654,6 +692,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( ...@@ -654,6 +692,8 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
p->force.pressure = pressure; p->force.pressure = pressure;
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
} }
/** /**
...@@ -683,7 +723,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -683,7 +723,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->u += p->u_dt * dt_therm; p->u += p->u_dt * dt_therm;
/* Check against absolute minimum */ /* Check against absolute minimum */
const float min_u = hydro_props->minimal_internal_energy; const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
p->u = max(p->u, min_u); p->u = max(p->u, min_u);
...@@ -713,6 +754,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -713,6 +754,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->force.pressure = pressure; p->force.pressure = pressure;
p->force.soundspeed = soundspeed; p->force.soundspeed = soundspeed;
p->force.v_sig = max(p->force.v_sig, 2.f * soundspeed);
} }
/** /**
...@@ -759,7 +802,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -759,7 +802,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full);
/* Check against absolute minimum */ /* Check against absolute minimum */
const float min_u = hydro_props->minimal_internal_energy; const float min_u =
hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy;
if (xp->u_full < min_u) { if (xp->u_full < min_u) {
xp->u_full = min_u; xp->u_full = min_u;
......
...@@ -56,6 +56,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -56,6 +56,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
float wi, wj, wi_dx, wj_dx; float wi, wj, wi_dx, wj_dx;
#ifdef SWIFT_DEBUG_CHECKS
if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
/* Get r. */ /* Get r. */
const float r_inv = 1.0f / sqrtf(r2); const float r_inv = 1.0f / sqrtf(r2);
const float r = r2 * r_inv; const float r = r2 * r_inv;
...@@ -83,6 +90,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density( ...@@ -83,6 +90,33 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx); pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
pj->density.wcount += wj; pj->density.wcount += wj;
pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx); pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
/* Compute dv dot r */
float dv[3], curlvr[3];
const float faci = mj * wi_dx * r_inv;
const float facj = mi * wj_dx * r_inv;
dv[0] = pi->v[0] - pj->v[0];
dv[1] = pi->v[1] - pj->v[1];
dv[2] = pi->v[2] - pj->v[2];
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
pi->density.div_v -= faci * dvdr;
pj->density.div_v -= facj * dvdr;
/* Compute dv cross r */
curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
pi->density.rot_v[0] += faci * curlvr[0];
pi->density.rot_v[1] += faci * curlvr[1];
pi->density.rot_v[2] += faci * curlvr[2];
pj->density.rot_v[0] += facj * curlvr[0];
pj->density.rot_v[1] += facj * curlvr[1];
pj->density.rot_v[2] += facj * curlvr[2];
} }
/** /**
...@@ -103,6 +137,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( ...@@ -103,6 +137,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float wi, wi_dx; float wi, wi_dx;
#ifdef SWIFT_DEBUG_CHECKS
if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
/* Get the masses. */ /* Get the masses. */
const float mj = pj->mass; const float mj = pj->mass;
...@@ -118,6 +159,27 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( ...@@ -118,6 +159,27 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx); pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
pi->density.wcount += wi; pi->density.wcount += wi;
pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx); pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
/* Compute dv dot r */
float dv[3], curlvr[3];
const float faci = mj * wi_dx * r_inv;
dv[0] = pi->v[0] - pj->v[0];
dv[1] = pi->v[1] - pj->v[1];
dv[2] = pi->v[2] - pj->v[2];
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
pi->density.div_v -= faci * dvdr;
/* Compute dv cross r */
curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
pi->density.rot_v[0] += faci * curlvr[0];
pi->density.rot_v[1] += faci * curlvr[1];
pi->density.rot_v[2] += faci * curlvr[2];
} }
/** /**
...@@ -136,6 +198,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -136,6 +198,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
float r2, const float *dx, float hi, float hj, struct part *restrict pi, float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, float a, float H) { struct part *restrict pj, float a, float H) {
#ifdef SWIFT_DEBUG_CHECKS
if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
/* Cosmological factors entering the EoMs */ /* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a); const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H; const float a2_Hubble = a * a * H;
...@@ -168,9 +237,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -168,9 +237,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
kernel_deval(xj, &wj, &wj_dx); kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx; const float wj_dr = hjd_inv * wj_dx;
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
/* Compute gradient terms */ /* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f; const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* Compute dv dot r. */ /* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
...@@ -195,7 +268,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( ...@@ -195,7 +268,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
/* Convolve with the kernel */ /* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; const float visc_acc_term =
0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
/* SPH acceleration term */ /* SPH acceleration term */
const float sph_acc_term = const float sph_acc_term =
...@@ -253,6 +327,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -253,6 +327,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
float r2, const float *dx, float hi, float hj, struct part *restrict pi, float r2, const float *dx, float hi, float hj, struct part *restrict pi,
const struct part *restrict pj, float a, float H) { const struct part *restrict pj, float a, float H) {
#ifdef SWIFT_DEBUG_CHECKS
if (pi->time_bin >= time_bin_inhibited)
error("Inhibited pi in interaction function!");
if (pj->time_bin >= time_bin_inhibited)
error("Inhibited pj in interaction function!");
#endif
/* Cosmological factors entering the EoMs */ /* Cosmological factors entering the EoMs */
const float fac_mu = pow_three_gamma_minus_five_over_two(a); const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H; const float a2_Hubble = a * a * H;
...@@ -262,7 +343,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -262,7 +343,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float r = r2 * r_inv; const float r = r2 * r_inv;
/* Recover some data */ /* Recover some data */
// const float mi = pi->mass; const float mi = pi->mass;
const float mj = pj->mass; const float mj = pj->mass;
const float rhoi = pi->rho; const float rhoi = pi->rho;
const float rhoj = pj->rho; const float rhoj = pj->rho;
...@@ -285,9 +366,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -285,9 +366,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
kernel_deval(xj, &wj, &wj_dx); kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx; const float wj_dr = hjd_inv * wj_dx;
/* Variable smoothing length term */
const float f_ij = 1.f - pi->force.f / mj;
const float f_ji = 1.f - pj->force.f / mi;
/* Compute gradient terms */ /* Compute gradient terms */
const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f; const float P_over_rho2_i = pressurei / (rhoi * rhoi) * f_ij;
const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f; const float P_over_rho2_j = pressurej / (rhoj * rhoj) * f_ji;
/* Compute dv dot r. */ /* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
...@@ -314,7 +399,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( ...@@ -314,7 +399,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij; const float visc = -0.25f * v_sig * mu_ij * (balsara_i + balsara_j) / rho_ij;
/* Convolve with the kernel */ /* Convolve with the kernel */
const float visc_acc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv; const float visc_acc_term =
0.5f * visc * (wi_dr * f_ij + wj_dr * f_ji) * r_inv;
/* SPH acceleration term */ /* SPH acceleration term */
const float sph_acc_term = const float sph_acc_term =
......
...@@ -51,7 +51,11 @@ INLINE static void hydro_read_particles(struct part* parts, ...@@ -51,7 +51,11 @@ INLINE static void hydro_read_particles(struct part* parts,
struct io_props* list, struct io_props* list,
int* num_fields) { int* num_fields) {
#ifdef PLANETARY_FIXED_ENTROPY
*num_fields = 10;
#else
*num_fields = 9; *num_fields = 9;
#endif
/* List what we want to read */ /* List what we want to read */
list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY, list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
...@@ -72,6 +76,11 @@ INLINE static void hydro_read_particles(struct part* parts, ...@@ -72,6 +76,11 @@ INLINE static void hydro_read_particles(struct part* parts,
UNIT_CONV_DENSITY, parts, rho); UNIT_CONV_DENSITY, parts, rho);
list[8] = io_make_input_field("MaterialIDs", INT, 1, COMPULSORY, list[8] = io_make_input_field("MaterialIDs", INT, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, parts, mat_id); UNIT_CONV_NO_UNITS, parts, mat_id);
#ifdef PLANETARY_FIXED_ENTROPY
list[9] = io_make_input_field("Entropies", FLOAT, 1, COMPULSORY,
UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS, parts,
s_fixed);
#endif
} }
INLINE static void convert_S(const struct engine* e, const struct part* p, INLINE static void convert_S(const struct engine* e, const struct part* p,
......
...@@ -208,6 +208,11 @@ struct part { ...@@ -208,6 +208,11 @@ struct part {
#endif #endif
#ifdef PLANETARY_FIXED_ENTROPY
/* Fixed specific entropy */
float s_fixed;
#endif
} SWIFT_STRUCT_ALIGN; } SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_PLANETARY_HYDRO_PART_H */ #endif /* SWIFT_PLANETARY_HYDRO_PART_H */
...@@ -329,6 +329,12 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5], ...@@ -329,6 +329,12 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5],
baseUnitsExp[UNIT_TIME] = -2.f; baseUnitsExp[UNIT_TIME] = -2.f;
break; break;
case UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS:
baseUnitsExp[UNIT_LENGTH] = 2.f;
baseUnitsExp[UNIT_TIME] = -2.f;
baseUnitsExp[UNIT_TEMPERATURE] = -1.f;
break;
case UNIT_CONV_POWER: case UNIT_CONV_POWER:
baseUnitsExp[UNIT_MASS] = 1.f; baseUnitsExp[UNIT_MASS] = 1.f;
baseUnitsExp[UNIT_LENGTH] = 2.f; baseUnitsExp[UNIT_LENGTH] = 2.f;
......
...@@ -84,6 +84,7 @@ enum unit_conversion_factor { ...@@ -84,6 +84,7 @@ enum unit_conversion_factor {
UNIT_CONV_ENERGY_PER_UNIT_MASS, UNIT_CONV_ENERGY_PER_UNIT_MASS,
UNIT_CONV_ENTROPY, UNIT_CONV_ENTROPY,
UNIT_CONV_ENTROPY_PER_UNIT_MASS, UNIT_CONV_ENTROPY_PER_UNIT_MASS,
UNIT_CONV_PHYSICAL_ENTROPY_PER_UNIT_MASS,
UNIT_CONV_POWER, UNIT_CONV_POWER,
UNIT_CONV_PRESSURE, UNIT_CONV_PRESSURE,
UNIT_CONV_FREQUENCY, UNIT_CONV_FREQUENCY,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment