Skip to content

Commit

Permalink
linear log interp combine metric
Browse files Browse the repository at this point in the history
  • Loading branch information
Mike Park authored and Mike Park committed Jun 30, 2021
1 parent 525812b commit 1d0faf0
Show file tree
Hide file tree
Showing 3 changed files with 27 additions and 71 deletions.
51 changes: 0 additions & 51 deletions src/ref_metric.c
View file
Expand Up @@ -2049,57 +2049,6 @@ REF_STATUS ref_metric_local_scale(REF_DBL *metric, REF_DBL *weight,

return REF_SUCCESS;
}
REF_STATUS ref_metric_scale_combine(REF_DBL *hess1, REF_DBL *hess2,
REF_DBL *metric, REF_GRID ref_grid,
REF_INT p_norm) {
REF_NODE ref_node = ref_grid_node(ref_grid);
REF_INT i, node;
REF_INT dimension;
REF_DBL det1, det2, exponent;

dimension = 3;
if (ref_grid_twod(ref_grid)) {
dimension = 2;
}

if (ref_grid_twod(ref_grid)) {
each_ref_node_valid_node(ref_node, node) {
hess1[2 + 6 * node] = 0.0;
hess1[4 + 6 * node] = 0.0;
hess1[5 + 6 * node] = 1.0;
}
}
if (ref_grid_twod(ref_grid)) {
each_ref_node_valid_node(ref_node, node) {
hess2[2 + 6 * node] = 0.0;
hess2[4 + 6 * node] = 0.0;
hess2[5 + 6 * node] = 1.0;
}
}

/* local scaling */
exponent = -1.0 / ((REF_DBL)(2 * p_norm + dimension));
each_ref_node_valid_node(ref_node, node) {
RSS(ref_matrix_det_m(&(hess1[6 * node]), &det1), "det_m local hess scale");
RSS(ref_matrix_det_m(&(hess2[6 * node]), &det2), "det_m local hess scale");
if (det1 > 0.0 && det2 > 0.0) {
for (i = 0; i < 6; i++) {
metric[i + 6 * node] = hess1[i + 6 * node] * pow(det1, exponent) +
hess2[i + 6 * node] * pow(det2, exponent);
}
}
}

if (ref_grid_twod(ref_grid)) {
each_ref_node_valid_node(ref_node, node) {
metric[2 + 6 * node] = 0.0;
metric[4 + 6 * node] = 0.0;
metric[5 + 6 * node] = 1.0;
}
}

return REF_SUCCESS;
}

REF_STATUS ref_metric_opt_goal(REF_DBL *metric, REF_GRID ref_grid,
REF_INT nequations, REF_DBL *solution,
Expand Down
3 changes: 0 additions & 3 deletions src/ref_metric.h
View file
Expand Up @@ -115,9 +115,6 @@ REF_STATUS ref_metric_eig_bal(REF_DBL *metric, REF_GRID ref_grid,
REF_DBL complexity);
REF_STATUS ref_metric_local_scale(REF_DBL *metric, REF_DBL *weight,
REF_GRID ref_grid, REF_INT p_norm);
REF_STATUS ref_metric_scale_combine(REF_DBL *hess1, REF_DBL *hess2,
REF_DBL *metric, REF_GRID ref_grid,
REF_INT p_norm);
REF_STATUS ref_metric_opt_goal(REF_DBL *metric, REF_GRID ref_grid,
REF_INT nequations, REF_DBL *solution,
REF_RECON_RECONSTRUCTION reconstruction,
Expand Down
44 changes: 27 additions & 17 deletions src/ref_metric_test.c
View file
Expand Up @@ -400,11 +400,12 @@ int main(int argc, char *argv[]) {

if (combine_pos != REF_EMPTY) {
REF_GRID ref_grid;
REF_DBL *scalar1, *scalar2, *hess1, *hess2, *metric;
REF_NODE ref_node;
REF_DBL *scalar1, *scalar2, *metric1, *metric2, *metric;
REF_INT p;
REF_DBL gradation, complexity, current_complexity;
REF_RECON_RECONSTRUCTION reconstruction = REF_RECON_L2PROJECTION;
REF_INT ldim;
REF_INT ldim, i, node;
REIS(1, combine_pos,
"required args: --combine grid.meshb scalar1.solb scalar2.solb p "
"gradation "
Expand Down Expand Up @@ -434,6 +435,7 @@ int main(int argc, char *argv[]) {
RSS(ref_part_by_extension(&ref_grid, ref_mpi, argv[2]),
"unable to load target grid in position 2");
ref_mpi_stopwatch_stop(ref_mpi, "read grid");
ref_node = ref_grid_node(ref_grid);

if (ref_mpi_once(ref_mpi)) printf("reading scalar1 %s\n", argv[3]);
RSS(ref_part_scalar(ref_grid, &ldim, &scalar1, argv[3]),
Expand All @@ -448,21 +450,29 @@ int main(int argc, char *argv[]) {
ref_mpi_stopwatch_stop(ref_mpi, "read scalar2");

ref_malloc(metric, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
ref_malloc(hess1, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
ref_malloc(hess2, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
ref_malloc(metric1, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);
ref_malloc(metric2, 6 * ref_node_max(ref_grid_node(ref_grid)), REF_DBL);

RSS(ref_recon_hessian(ref_grid, scalar1, hess1, reconstruction), "recon");
RSS(ref_recon_hessian(ref_grid, scalar2, hess2, reconstruction), "recon");
RSS(ref_recon_roundoff_limit(hess1, ref_grid),
"floor metric eigenvalues based on grid size and solution jitter");
ref_mpi_stopwatch_stop(ref_mpi, "recon hess1");
RSS(ref_recon_roundoff_limit(hess2, ref_grid),
"floor metric eigenvalues based on grid size and solution jitter");
ref_mpi_stopwatch_stop(ref_mpi, "recon hess2");
RSS(ref_metric_lp(metric1, ref_grid, scalar1, NULL, reconstruction, p,
gradation, complexity),
"lp norm");
ref_mpi_stopwatch_stop(ref_mpi, "multiscale metric1");
RSS(ref_metric_lp(metric2, ref_grid, scalar2, NULL, reconstruction, p,
gradation, complexity),
"lp norm");
ref_mpi_stopwatch_stop(ref_mpi, "multiscale metric2");

RSS(ref_metric_scale_combine(hess1, hess2, metric, ref_grid, p),
"local scale lp norm");
ref_mpi_stopwatch_stop(ref_mpi, "combine hess");
each_ref_node_valid_node(ref_grid_node(ref_grid), node) {
REF_DBL log_m1[6];
REF_DBL log_m2[6];
REF_DBL log_m[6];
REF_DBL s;
s = 0.5;
RSS(ref_matrix_log_m(&(metric1[6 * node]), log_m1), "log");
RSS(ref_matrix_log_m(&(metric2[6 * node]), log_m2), "log");
for (i = 0; i < 6; i++) log_m[i] = (1.0 - s) * log_m1[i] + s * log_m2[i];
RSS(ref_matrix_exp_m(log_m, &(metric[6 * node])), "exp");
}

RSS(ref_metric_gradation_at_complexity(metric, ref_grid, gradation,
complexity),
Expand All @@ -475,8 +485,8 @@ int main(int argc, char *argv[]) {
printf("actual complexity %e\n", current_complexity);
RSS(ref_metric_to_node(metric, ref_grid_node(ref_grid)), "set node");
ref_free(metric);
ref_free(hess2);
ref_free(hess1);
ref_free(metric2);
ref_free(metric1);
ref_free(scalar2);
ref_free(scalar1);

Expand Down

0 comments on commit 1d0faf0

Please sign in to comment.