creates a variable aspect ratio limit

nasa · Jan 17, 2022 · 04b65d0 · 04b65d0
1 parent 6d5dda4
commit 04b65d0
Show file tree

Hide file tree

Showing 2 changed files with 50 additions and 0 deletions.
diff --git a/src/ref_metric.c b/src/ref_metric.c
@@ -1870,6 +1870,53 @@ REF_STATUS ref_metric_limit_aspect_ratio(REF_DBL *metric, REF_GRID ref_grid,
   return REF_SUCCESS;
 }
 
+REF_STATUS ref_metric_limit_aspect_ratio_field(REF_DBL *metric,
+                                               REF_GRID ref_grid,
+                                               REF_DBL *aspect_ratio_field) {
+  REF_DBL diag_system[12];
+  REF_DBL max_eig, limit_eig;
+  REF_INT node;
+  REF_DBL aspect_ratio;
+  if (ref_grid_twod(ref_grid)) {
+    each_ref_node_valid_node(ref_grid_node(ref_grid), node) {
+      aspect_ratio = aspect_ratio_field[node];
+      RSS(ref_matrix_diag_m(&(metric[6 * node]), diag_system), "eigen decomp");
+      RSS(ref_matrix_ascending_eig_twod(diag_system), "2D eig sort");
+      max_eig = ref_matrix_eig(diag_system, 0);
+      max_eig = MAX(ref_matrix_eig(diag_system, 1), max_eig);
+      if (ref_math_divisible(max_eig, (aspect_ratio * aspect_ratio))) {
+        limit_eig = max_eig / (aspect_ratio * aspect_ratio);
+        ref_matrix_eig(diag_system, 0) =
+            MAX(ref_matrix_eig(diag_system, 0), limit_eig);
+        ref_matrix_eig(diag_system, 1) =
+            MAX(ref_matrix_eig(diag_system, 1), limit_eig);
+        RSS(ref_matrix_form_m(diag_system, &(metric[6 * node])), "reform m");
+        RSS(ref_matrix_twod_m(&(metric[6 * node])),
+            "enforce 2D m for roundoff");
+      }
+    }
+  } else {
+    each_ref_node_valid_node(ref_grid_node(ref_grid), node) {
+      aspect_ratio = aspect_ratio_field[node];
+      RSS(ref_matrix_diag_m(&(metric[6 * node]), diag_system), "eigen decomp");
+      max_eig = ref_matrix_eig(diag_system, 0);
+      max_eig = MAX(ref_matrix_eig(diag_system, 1), max_eig);
+      max_eig = MAX(ref_matrix_eig(diag_system, 2), max_eig);
+      if (ref_math_divisible(max_eig, (aspect_ratio * aspect_ratio))) {
+        limit_eig = max_eig / (aspect_ratio * aspect_ratio);
+        ref_matrix_eig(diag_system, 0) =
+            MAX(ref_matrix_eig(diag_system, 0), limit_eig);
+        ref_matrix_eig(diag_system, 1) =
+            MAX(ref_matrix_eig(diag_system, 1), limit_eig);
+        ref_matrix_eig(diag_system, 2) =
+            MAX(ref_matrix_eig(diag_system, 2), limit_eig);
+        RSS(ref_matrix_form_m(diag_system, &(metric[6 * node])), "reform m");
+      }
+    }
+  }
+  return REF_SUCCESS;
+}
+
 REF_STATUS ref_metric_limit_h(REF_DBL *metric, REF_GRID ref_grid, REF_DBL hmin,
                               REF_DBL hmax) {
   REF_DBL diag_system[12];

diff --git a/src/ref_metric.h b/src/ref_metric.h
@@ -97,6 +97,9 @@ REF_STATUS ref_metric_set_complexity(REF_DBL *metric, REF_GRID ref_grid,
                                      REF_DBL target_complexity);
 REF_STATUS ref_metric_limit_aspect_ratio(REF_DBL *metric, REF_GRID ref_grid,
                                          REF_DBL aspect_ratio);
+REF_STATUS ref_metric_limit_aspect_ratio_field(REF_DBL *metric,
+                                               REF_GRID ref_grid,
+                                               REF_DBL *aspect_ratio_field);
 REF_STATUS ref_metric_limit_h(REF_DBL *metric, REF_GRID ref_grid, REF_DBL hmin,
                               REF_DBL hmax);
 REF_STATUS ref_metric_limit_h_at_complexity(REF_DBL *metric, REF_GRID ref_grid,