Merge branch 'iss642' into 'master'

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

Closes #642

See merge request !1005

src/equation_of_state/planetary/sesame.h

@@ -338,15 +338,35 @@ INLINE static float SESAME_pressure_from_internal_energy(
     idx_u_2 = mat->num_T - 2;
   }
 
-  intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
-             (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]) /
-             (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_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]);
+  // 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.;
+  }
+
+  // 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
   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(
     idx_u_2 = mat->num_T - 2;
   }
 
-  intp_rho = (log_rho - mat->table_log_rho[idx_rho]) /
-             (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]) /
-             (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_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]);
+  // 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.;
+  }
+
+  // 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
   c_1 = mat->table_c_rho_T[idx_rho * mat->num_T + idx_u_1];