Merge branch 'diag-oct-metric' into 'main'

diagnose roundoff oct

See merge request fun3d-developers/refine!963

src/ref_oct_test.c

@@ -76,7 +76,7 @@ int main(int argc, char *argv[]) {
     REF_GRID ref_grid;
     REF_OCT ref_oct;
     REF_DBL h;
-    REF_DBL bbox[]={0.0,1.0,0.0,1.0,0,0.03};
+    REF_DBL bbox[] = {0.0, 1.0, 0.0, 1.0, 0, 0.03};
     REF_INT nleaf;
     h = atof(argv[pos + 1]);
     RSS(ref_grid_create(&ref_grid, ref_mpi), "make grid");
@@ -162,8 +162,8 @@ int main(int argc, char *argv[]) {
   RXS(ref_args_find(argc, argv, "--adapt", &pos), REF_NOT_FOUND, "arg search");
   if (pos != REF_EMPTY && pos + 3 < argc) {
     REF_OCT ref_oct;
-    char tec[1024];
     REF_INT nleaf;
+    REF_BOOL debug = REF_FALSE;
 
     RSS(ref_oct_create(&ref_oct), "make oct");
     {
@@ -190,15 +190,24 @@ int main(int argc, char *argv[]) {
       }
       RSS(ref_grid_free(ref_grid), "free grid");
     }
-    snprintf(tec, 1024, "%s-raw.tec", argv[pos + 3]);
+
     RSS(ref_oct_nleaf(ref_oct, &nleaf), "count leaves");
-    printf("writing %d vox to %s from %s\n", nleaf, tec, argv[pos + 3]);
-    RSS(ref_oct_tec(ref_oct, tec), "tec");
+    if (debug) {
+      char tec[1024];
+      snprintf(tec, 1024, "%s-raw.tec", argv[pos + 3]);
+      printf("writing %d vox to %s from %s\n", nleaf, tec, argv[pos + 3]);
+      RSS(ref_oct_tec(ref_oct, tec), "tec");
+    } else {
+      printf("raw vox %d\n", nleaf);
+    }
     RSS(ref_oct_gradation(ref_oct), "grad");
-    snprintf(tec, 1024, "%s-grad.tec", argv[pos + 3]);
     RSS(ref_oct_nleaf(ref_oct, &nleaf), "count leaves");
-    printf("writing %d vox to %s from %s\n", nleaf, tec, argv[pos + 3]);
-    RSS(ref_oct_tec(ref_oct, tec), "tec");
+    if (debug) {
+      char tec[1024];
+      snprintf(tec, 1024, "%s-grad.tec", argv[pos + 3]);
+      printf("writing %d vox to %s from %s\n", nleaf, tec, argv[pos + 3]);
+      RSS(ref_oct_tec(ref_oct, tec), "tec");
+    }
     printf("writing %d vox to %s\n", nleaf, argv[pos + 3]);
     {
       REF_GRID ref_grid;

src/ref_recon.c

@@ -1265,7 +1265,11 @@ REF_FCN REF_STATUS ref_recon_roundoff_limit(REF_DBL *recon, REF_GRID ref_grid) {
         "element with zero edge length or missing element");
     eig_floor = (4 * round_off_jitter) / (radius[node] * radius[node]);
 
-    RSS(ref_matrix_diag_m(&(recon[6 * node]), diag_system), "eigen decomp");
+    RSB(ref_matrix_diag_m(&(recon[6 * node]), diag_system), "eigen decomp", {
+      ref_matrix_show_m(&(recon[6 * node]));
+      ref_node_location(ref_node, node);
+      printf("eig_floor %e radius %e\n", eig_floor, radius[node]);
+    });
     ref_matrix_eig(diag_system, 0) =
         MAX(ref_matrix_eig(diag_system, 0), eig_floor);
     ref_matrix_eig(diag_system, 1) =

src/ref_subcommand.c

@@ -2612,6 +2612,13 @@ static REF_STATUS fixed_point_metric(
   REF_DBL inv_total;
   REF_INT im, node;
   REF_INT fixed_point_ldim;
+  REF_BOOL ensure_finite = REF_TRUE;
+
+  each_ref_node_valid_node(ref_grid_node(ref_grid), node) {
+    for (im = 0; im < 6; im++) {
+      metric[im + 6 * node] = 0.0; /* initialize */
+    }
+  }
 
   ref_malloc(hess, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
   total_timesteps = 0;
@@ -2624,19 +2631,29 @@ static REF_STATUS fixed_point_metric(
     RSS(ref_part_scalar(ref_grid, &fixed_point_ldim, &scalar, solb_filename),
         "unable to load scalar");
     REIS(1, fixed_point_ldim, "expected one scalar");
+    if (ensure_finite)
+      RSS(ref_validation_finite(ref_grid, fixed_point_ldim, scalar),
+          "input scalar");
     if (strong_sensor_bc) {
       RSS(ref_phys_strong_sensor_bc(ref_grid, scalar, strong_value,
                                     ref_dict_bcs),
           "apply strong sensor bc");
+      if (ensure_finite)
+        RSS(ref_validation_finite(ref_grid, fixed_point_ldim, scalar),
+            "strong scalar");
     }
     RSS(ref_recon_hessian(ref_grid, scalar, hess, reconstruction), "hess");
+    if (ensure_finite)
+      RSS(ref_validation_finite(ref_grid, 6, hess), "recon hess");
     ref_free(scalar);
     total_timesteps++;
     each_ref_node_valid_node(ref_grid_node(ref_grid), node) {
       for (im = 0; im < 6; im++) {
         metric[im + 6 * node] += hess[im + 6 * node];
       }
     }
+    if (ensure_finite)
+      RSS(ref_validation_finite(ref_grid, 6, metric), "metric sum");
   }
   free(hess);
   ref_mpi_stopwatch_stop(ref_mpi, "all timesteps processed");
@@ -2648,6 +2665,8 @@ static REF_STATUS fixed_point_metric(
       metric[im + 6 * node] *= inv_total;
     }
   }
+  if (ensure_finite)
+    RSS(ref_validation_finite(ref_grid, 6, metric), "metric avg");
 
   RSS(ref_recon_roundoff_limit(metric, ref_grid),
       "floor metric eigenvalues based on grid size and solution jitter");