Skip to content

Commit

Permalink
combines two scalar field hessians into one metric
Browse files Browse the repository at this point in the history
  • Loading branch information
Mike Park authored and Mike Park committed Jun 23, 2021
1 parent 96742b2 commit eabaea1
Show file tree
Hide file tree
Showing 3 changed files with 149 additions and 0 deletions.
51 changes: 51 additions & 0 deletions src/ref_metric.c
View file
Expand Up @@ -1985,6 +1985,57 @@ 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: 3 additions & 0 deletions src/ref_metric.h
View file
Expand Up @@ -109,6 +109,9 @@ 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
95 changes: 95 additions & 0 deletions src/ref_metric_test.c
View file
Expand Up @@ -70,6 +70,7 @@ int main(int argc, char *argv[]) {
REF_INT explore_pos = REF_EMPTY;
REF_INT multigrad_pos = REF_EMPTY;
REF_INT lp_pos = REF_EMPTY;
REF_INT combine_pos = REF_EMPTY;
REF_INT opt_goal_pos = REF_EMPTY;
REF_INT no_goal_pos = REF_EMPTY;
REF_INT venditti_pos = REF_EMPTY;
Expand Down Expand Up @@ -102,6 +103,8 @@ int main(int argc, char *argv[]) {
RXS(ref_args_find(argc, argv, "--wlp", &wlp_pos), REF_NOT_FOUND,
"arg search");
RXS(ref_args_find(argc, argv, "--lp", &lp_pos), REF_NOT_FOUND, "arg search");
RXS(ref_args_find(argc, argv, "--combine", &combine_pos), REF_NOT_FOUND,
"arg search");
RXS(ref_args_find(argc, argv, "--multigrad", &multigrad_pos), REF_NOT_FOUND,
"arg search");
RXS(ref_args_find(argc, argv, "--moving", &moving_pos), REF_NOT_FOUND,
Expand Down Expand Up @@ -389,6 +392,98 @@ int main(int argc, char *argv[]) {
return 0;
}

if (combine_pos != REF_EMPTY) {
REF_GRID ref_grid;
REF_DBL *scalar1, *scalar2, *hess1, *hess2, *metric;
REF_INT p;
REF_DBL gradation, complexity, current_complexity;
REF_RECON_RECONSTRUCTION reconstruction = REF_RECON_L2PROJECTION;
REF_INT ldim;
REIS(1, combine_pos,
"required args: --combine grid.meshb scalar1.solb scalar2.solb p "
"gradation "
"complexity output-metric.solb");
if (8 > argc) {
printf(
"required args: --combine grid.meshb scalar1.solb scalar2.solb p "
"gradation "
"complexity output-metric.solb\n");
return REF_FAILURE;
}

p = atoi(argv[5]);
gradation = atof(argv[6]);
complexity = atof(argv[7]);
if (REF_EMPTY != kexact_pos) {
reconstruction = REF_RECON_KEXACT;
}
if (ref_mpi_once(ref_mpi)) {
printf("Lp=%d\n", p);
printf("gradation %f\n", gradation);
printf("complexity %f\n", complexity);
printf("reconstruction %d\n", (int)reconstruction);
}

if (ref_mpi_once(ref_mpi)) printf("reading grid %s\n", argv[2]);
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");

if (ref_mpi_once(ref_mpi)) printf("reading scalar1 %s\n", argv[3]);
RSS(ref_part_scalar(ref_grid, &ldim, &scalar1, argv[3]),
"unable to load scalar1 in position 3");
REIS(1, ldim, "expected one scalar1");
ref_mpi_stopwatch_stop(ref_mpi, "read scalar1");

if (ref_mpi_once(ref_mpi)) printf("reading scalar2 %s\n", argv[4]);
RSS(ref_part_scalar(ref_grid, &ldim, &scalar2, argv[4]),
"unable to load scalar2 in position 4");
REIS(1, ldim, "expected one scalar2");
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);

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_scale_combine(hess1, hess2, metric, ref_grid, p),
"local scale lp norm");
ref_mpi_stopwatch_stop(ref_mpi, "combine hess");

RSS(ref_metric_gradation_at_complexity(metric, ref_grid, gradation,
complexity),
"gradation at complexity");

ref_mpi_stopwatch_stop(ref_mpi, "metric gradation");

RSS(ref_metric_complexity(metric, ref_grid, &current_complexity), "cmp");
if (ref_mpi_once(ref_mpi))
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(scalar2);
ref_free(scalar1);

if (ref_mpi_once(ref_mpi)) printf("writing metric %s\n", argv[8]);
RSS(ref_gather_metric(ref_grid, argv[7]), "export curve limit metric");
ref_mpi_stopwatch_stop(ref_mpi, "write metric");

RSS(ref_grid_free(ref_grid), "free");
RSS(ref_mpi_free(ref_mpi), "free");
RSS(ref_mpi_stop(), "stop");
return 0;
}

if (multigrad_pos != REF_EMPTY) {
REF_GRID ref_grid;
REF_DBL *grad, *metric;
Expand Down

0 comments on commit eabaea1

Please sign in to comment.