accesses error est

nasa · Jun 29, 2021 · c61401f · c61401f
1 parent 454639c
commit c61401f
Showing 1 changed file with 71 additions and 0 deletions.
diff --git a/src/ref_metric_test.c b/src/ref_metric_test.c
@@ -89,6 +89,7 @@ int main(int argc, char *argv[]) {
   REF_INT decompose_pos = REF_EMPTY;
   REF_INT imply_pos = REF_EMPTY;
   REF_INT eigs_pos = REF_EMPTY;
+  REF_INT error_pos = REF_EMPTY;
 
   REF_MPI ref_mpi;
   RSS(ref_mpi_start(argc, argv), "start");
@@ -148,6 +149,8 @@ int main(int argc, char *argv[]) {
       "arg search");
   RXS(ref_args_find(argc, argv, "--eigs", &eigs_pos), REF_NOT_FOUND,
       "arg search");
+  RXS(ref_args_find(argc, argv, "--error", &error_pos), REF_NOT_FOUND,
+      "arg search");
 
   if (curve_limit_pos != REF_EMPTY) {
     REF_GRID ref_grid;
@@ -1962,6 +1965,74 @@ int main(int argc, char *argv[]) {
     return 0;
   }
 
+  if (error_pos != REF_EMPTY) {
+    REF_GRID ref_grid;
+    REF_DBL *field, *hess, *metric, *error;
+    REF_INT ldim;
+
+    REIS(1, error_pos,
+         "required args: --error grid.ext scalar.solb metric.solb");
+    REIS(5, argc, "required args: --error grid.ext scalar.solb metric.solb");
+
+    if (ref_mpi_once(ref_mpi)) printf("reading grid %s\n", argv[2]);
+    if (ref_mpi_para(ref_mpi)) {
+      if (ref_mpi_once(ref_mpi)) printf("part %s\n", argv[2]);
+      RSS(ref_part_by_extension(&ref_grid, ref_mpi, argv[2]), "part");
+      ref_mpi_stopwatch_stop(ref_mpi, "part mesh");
+    } else {
+      if (ref_mpi_once(ref_mpi)) printf("import %s\n", argv[2]);
+      RSS(ref_import_by_extension(&ref_grid, ref_mpi, argv[2]), "import");
+      ref_mpi_stopwatch_stop(ref_mpi, "import mesh");
+    }
+    ref_mpi_stopwatch_stop(ref_mpi, "read grid");
+
+    if (ref_mpi_once(ref_mpi)) printf("reading scalar %s\n", argv[3]);
+    RSS(ref_part_scalar(ref_grid, &ldim, &field, argv[3]),
+        "unable to load solution in position 3");
+    if (ref_mpi_once(ref_mpi)) printf("ldim %d\n", ldim);
+    ref_mpi_stopwatch_stop(ref_mpi, "read scalar");
+    REIS(1, ldim, "expect scalar");
+
+    if (ref_mpi_once(ref_mpi)) printf("reading metric %s\n", argv[4]);
+    RSS(ref_part_metric(ref_grid_node(ref_grid), argv[4]),
+        "unable to load metric in position 4");
+    ref_mpi_stopwatch_stop(ref_mpi, "read metric");
+
+    ref_malloc(metric, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
+    RSS(ref_metric_from_node(metric, ref_grid_node(ref_grid)), "from");
+
+    ref_malloc(hess, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
+    ldim = 1;
+    ref_malloc(error, ldim * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
+
+    RSS(ref_recon_hessian(ref_grid, field, hess, REF_RECON_L2PROJECTION),
+        "hess");
+
+    RSS(ref_metric_interpolation_error(metric, hess, ref_grid, error), "error")
+
+    /*
+    each_ref_node_valid_node(ref_grid_node(ref_grid), node) {
+      REF_DBL diag[12];
+      eigs[1 + ldim * node] = field[node];
+     }
+    */
+
+    RSS(ref_gather_scalar_by_extension(ref_grid, ldim, error, NULL,
+                                       "ref_metric_test_error.plt"),
+        "export eigs");
+
+    ref_free(error);
+    ref_free(hess);
+    ref_free(metric);
+    ref_free(field);
+
+    RSS(ref_grid_free(ref_grid), "free");
+    ref_mpi_stopwatch_stop(ref_mpi, "done.");
+    RSS(ref_mpi_free(ref_mpi), "free");
+    RSS(ref_mpi_stop(), "stop");
+    return 0;
+  }
+
   if (argc == 3) {
     REF_GRID ref_grid;