reports interpolation error estimate of multiscale

nasa · Jun 30, 2021 · f3cd9af · f3cd9af
1 parent 6f296ae
commit f3cd9af
Show file tree

Hide file tree

Showing 2 changed files with 23 additions and 0 deletions.
diff --git a/src/ref_metric.h b/src/ref_metric.h
@@ -82,6 +82,9 @@ REF_STATUS ref_metric_interpolation_error(REF_DBL *metric, REF_DBL *hess,
 REF_STATUS ref_metric_integrate_error(REF_GRID ref_grid,
                                       REF_DBL *interpolation_error,
                                       REF_DBL *total_error);
+REF_STATUS ref_metric_report_error(REF_DBL *metric, REF_GRID ref_grid,
+                                   REF_DBL *scalar,
+                                   REF_RECON_RECONSTRUCTION reconstruction);
 REF_STATUS ref_metric_complexity(REF_DBL *metric, REF_GRID ref_grid,
                                  REF_DBL *complexity);
 REF_STATUS ref_metric_set_complexity(REF_DBL *metric, REF_GRID ref_grid,

diff --git a/src/ref_subcommand.c b/src/ref_subcommand.c
@@ -2217,6 +2217,24 @@ static REF_STATUS parse_p(int argc, char *argv[], REF_INT *p) {
   return REF_SUCCESS;
 }
 
+static REF_STATUS ref_subcommand_report_error(
+    REF_DBL *metric, REF_GRID ref_grid, REF_DBL *scalar,
+    REF_RECON_RECONSTRUCTION reconstruction) {
+  REF_DBL *hess, *error;
+  REF_DBL total_error;
+  ref_malloc(hess, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
+  ref_malloc(error, ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
+
+  RSS(ref_recon_hessian(ref_grid, scalar, hess, reconstruction), "hess");
+  RSS(ref_metric_interpolation_error(metric, hess, ref_grid, error), "error")
+  RSS(ref_metric_integrate_error(ref_grid, error, &total_error), "int")
+  ref_free(error);
+  ref_free(hess);
+  if (ref_mpi_once(ref_grid_mpi(ref_grid)))
+    printf("interpolation error %e\n", total_error);
+  return REF_SUCCESS;
+}
+
 static REF_STATUS loop(REF_MPI ref_mpi_orig, int argc, char *argv[]) {
   char *in_project = NULL;
   char *out_project = NULL;
@@ -2666,6 +2684,8 @@ static REF_STATUS loop(REF_MPI ref_mpi_orig, int argc, char *argv[]) {
                         gradation, complexity),
           "lp norm");
       ref_mpi_stopwatch_stop(ref_mpi, "multiscale metric");
+      RSS(ref_subcommand_report_error(metric, ref_grid, scalar, reconstruction),
+          "report error");
     }
     ref_free(scalar);
   }