Strange meshes with `compute_anisotropic_dwr_metric`

[ddundo]

I mentioned on the meeting yesterday that I am getting strange meshes when using RiemannianMetric.compute_anisotropic_dwr_metric and I haven't been able to get to the bottom of it.

I have 2 prognostic fields: thickness h and velocity u, and corresponding error indicators h_ind and u_ind. Computing compute_isotropic_dwr_metric(indicator) and adapting the mesh gives expected results. But when I compute hessians of h and u, and intersect them to compute compute_anisotropic_dwr_metric(indicator, intersected_hessian), I get an adapted mesh which looks as if the hessian acts against the indicator. Fine resolution gets prescribed in the right-hand side of the domain, where the hessian metric density is lowest.

Here is the code used to get the above:

# load msh, h, u, h_ind, u_ind from a checkpoint file

mp = {"dm_plex_metric": 
    {
        "target_complexity": 1600,
        "h_min": 1.0,
        "h_max": 50e3,
        "p": 2.0,
        "a_max": 1e30,
    }
}

P1_ten = TensorFunctionSpace(msh, "CG", 1)
metric = RiemannianMetric(P1_ten)
metric.set_parameters(mp)

# isotropic dwr metric from h_ind
metric_hiso = metric.copy(deepcopy=True)
metric_hiso.compute_isotropic_dwr_metric(h_ind)
metric_hiso.normalise()
amsh_hiso = adapt(msh, metric_hiso)

# isotropic dwr metric from u_ind
metric_uiso = metric.copy(deepcopy=True)
metric_uiso.compute_isotropic_dwr_metric(u_ind)
metric_uiso.normalise()
amsh_uiso = adapt(msh, metric_uiso)

# Compute hessians
hessian_h = metric.copy(deepcopy=True)
hessian_h.compute_hessian(h)
hessian_h.enforce_spd(restrict_sizes=True, restrict_anisotropy=True)

hessian_u = metric.copy(deepcopy=True)
hessian_u1 = metric.copy(deepcopy=True)
for i, hessian in enumerate([hessian_u, hessian_u1]):
    hessian.compute_hessian(u[i])
    hessian.enforce_spd(restrict_sizes=True, restrict_anisotropy=True)
hessian_u.intersect(hessian_u1)
hessian_u.enforce_spd(restrict_sizes=True, restrict_anisotropy=True)

intersected_hessian = hessian_h.copy(deepcopy=True).intersect(hessian_u)
# metric density of the intersected hessian
rho, _, _ = intersected_hessian.density_and_quotients()

# anisotropic dwr metric from h_ind and intersected_hessian
metric_haniso = metric.copy(deepcopy=True)
metric_haniso.compute_anisotropic_dwr_metric(h_ind, intersected_hessian)
amsh_haniso = adapt(msh, metric_haniso)

# anisotropic dwr metric from u_ind and intersected_hessian
metric_uaniso = metric.copy(deepcopy=True)
metric_uaniso.compute_anisotropic_dwr_metric(u_ind, intersected_hessian)
amsh_uaniso = adapt(msh, metric_uaniso)

Could you please take a look if I did something wrong? :)

[ddundo]

By the way, compute_weighted_hessian_metric seems to give reasonable meshes:

weighted_hessian_metric = metric.copy(deepcopy=True)
weighted_hessian_metric.compute_weighted_hessian_metric([h_ind, u_ind], [hessian_h, hessian_u])
weighted_hessian_metric.normalise()
amsh_weighted = adapt(msh, weighted_hessian_metric)

[jwallwork23]

Hi @ddundo, thanks for reporting this. Sorry to hear you're having issues. The code looks fine to me.

I'm afraid RiemannianMetric.compute_anisotropic_dwr_metric is not at all well-tested. Currently, the only tests we have are (a) error handling, (b) a test in the case of a uniform indicator, and (c) any demos it appears in. Sadly, these are not enough to check it is doing everything we expect. It's on the to do list but I'm afraid it won't be addressed any time soon.

[jwallwork23]

By the way, compute_weighted_hessian_metric seems to give reasonable meshes:

Note that the weighted Hessian metric is meant to be used a bit differently. The error indicator parts don't mean the same thing. See (Power et al., 2006) or (Wallwork et al., 2020) for details.

[ddundo]

Thanks @jwallwork23!