Skip to content
Snippets Groups Projects
Commit 16d1bb08 authored by Sergio Ruiz-Bonilla's avatar Sergio Ruiz-Bonilla Committed by Matthieu Schaller
Browse files

Fix FPE from rare duplicate values in the SESAME planetary EoS tables

parent 35eaa463
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/equation_of_state/planetary/sesame.h
+ 58
− 18
View file @ 16d1bb08
...@@ -338,15 +338,35 @@ INLINE static float SESAME_pressure_from_internal_energy( ...@@ -338,15 +338,35 @@ INLINE static float SESAME_pressure_from_internal_energy(
idx_u_2 = mat->num_T - 2; idx_u_2 = mat->num_T - 2;
} }
intp_rho = (log_rho - mat->table_log_rho[idx_rho]) / // Check for duplicates in SESAME tables before interpolation
(mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]); if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) / intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
(mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] - (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]); } else {
intp_u_2 = intp_rho = 1.;
(log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) / }
(mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]); // Check for duplicates in SESAME tables before interpolation
if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] !=
mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) {
intp_u_1 =
(log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
(mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
} else {
intp_u_1 = 1.;
}
// Check for duplicates in SESAME tables before interpolation
if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
intp_u_2 =
(log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
(mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
} else {
intp_u_2 = 1.;
}
// Table values // Table values
P_1 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1]; P_1 = mat->table_P_rho_T[idx_rho * mat->num_T + idx_u_1];
...@@ -442,15 +462,35 @@ INLINE static float SESAME_soundspeed_from_internal_energy( ...@@ -442,15 +462,35 @@ INLINE static float SESAME_soundspeed_from_internal_energy(
idx_u_2 = mat->num_T - 2; idx_u_2 = mat->num_T - 2;
} }
intp_rho = (log_rho - mat->table_log_rho[idx_rho]) / // Check for duplicates in SESAME tables before interpolation
(mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]); if (mat->table_log_rho[idx_rho + 1] != mat->table_log_rho[idx_rho]) {
intp_u_1 = (log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) / intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
(mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] - (mat->table_log_rho[idx_rho + 1] - mat->table_log_rho[idx_rho]);
mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]); } else {
intp_u_2 = intp_rho = 1.;
(log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) / }
(mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]); // Check for duplicates in SESAME tables before interpolation
if (mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] !=
mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) {
intp_u_1 =
(log_u - mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]) /
(mat->table_log_u_rho_T[idx_rho * mat->num_T + (idx_u_1 + 1)] -
mat->table_log_u_rho_T[idx_rho * mat->num_T + idx_u_1]);
} else {
intp_u_1 = 1.;
}
// Check for duplicates in SESAME tables before interpolation
if (mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] !=
mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) {
intp_u_2 =
(log_u - mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]) /
(mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + (idx_u_2 + 1)] -
mat->table_log_u_rho_T[(idx_rho + 1) * mat->num_T + idx_u_2]);
} else {
intp_u_2 = 1.;
}
// Table values // Table values
c_1 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1]; c_1 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment