Getting very different meshes with anisotropic dwr metric

[ddundo]

Hey! In my glacier experiments, in my get_solver function I solve for the thickness h and then for the velocity u. Then following the discussion in #55, I modified the error estimation code in go_mesh_seq.py to the following:

          # Loop over each timestep
            for j in range(len(sols[f]["forward"][i])):
                for f in self.fields:
                    # Update fields
                    transfer(sols[f][FWD][i][j], u[f])
                    transfer(sols[f][FWD_OLD][i][j], u_[f])
                    transfer(sols[f][ADJ][i][j], u_star[f])
                    transfer(sols[f][ADJ_NEXT][i][j], u_star_next[f])

                    # Combine adjoint solutions as appropriate
                    u_star[f].assign(0.5 * (u_star[f] + u_star_next[f]))
                    u_star_e[f].assign(
                        0.5 * (sols_e[f][ADJ][i][j] + sols_e[f][ADJ_NEXT][i][j])
                    )
                    u_star_e[f] -= u_star[f]

                for f in self.fields:
                    # Stephan's comment. In get_solver we first solve for h, then u, so
                    # the h solver uses u at the time step previous to the regularly exported one
                    if f == 'h':
                        transfer(self.u_prev[i][j], u['u'])
                    # Evaluate error indicator
                    indi_e = indicator_fn(forms[f], u_star_e[f])

                    # Project back to the base space
                    indi = project(indi_e, P0_spaces[i])
                    indi.interpolate(abs(indi))
                    indicators[f][i][j].interpolate(ufl.max_value(indi, 1.0e-16))

where self.u_prev is a list of velocity fields (deep copied) at timestep=export tilmestep-1, which is my dirty/temporary way of addressing @stephankramer's comment in #55. I get thickness and velocity fields, and their indicators, which look like this at the first tilmestep (similarly for other timesteps):

This all looks good (highest error around the grounding line). And then I compute isotropic and anisotropic metrics for each, but I get quite a different mesh with anisotropic metric for velocity, which doesn't look so good and doesn't improve with further fixed point iterations. The one for thickness doesn't look amazing either for that matter, but still better.

I am trying to figure out what is wrong. I guess it is either the computing of error indicators code above (which I don't think so because the indicator looks fine), or my construction of the metric, which I obtain like this:

mp = {
    "dm_plex_metric_target_complexity": target_complexity,
    "dm_plex_metric_h_min": 1.,
    "dm_plex_metric_h_max": 50e3,
    "dm_plex_metric_p": 1.0,
    "dm_plex_metric_a_max": 1e5,
    }

P1_ten = TensorFunctionSpace(msq[0], "CG", 1)
metric = RiemannianMetric(P1_ten)
metric.set_parameters(mp)

u_ind = indicators['u'][0][0]
u_sols = solutions["u"]["forward"][0][0]
u_mag = Function(msq.function_spaces["h"][0])
u_mag.interpolate(sqrt(inner(u_sols, u_sols)))

hessian = RiemannianMetric(P1_ten)
hessian.compute_hessian(u_mag)
metric.compute_anisotropic_dwr_metric(u_ind, hessian)
metric.normalise()

But I'm not sure I see an issue here either. I construct the thickness metric in the same way. Here indicators and solutions are the output of go_mesh_seq.indicate_errors, so it does not touch my touch msq.u_prev list.

Does someone have an idea what it could be please? :)

[ddundo]

And just for completeness, here is Hessian-based adaptation using velocity (3142 cells):

[jwallwork23]

Sorry, I only just saw this. Some thoughts:

• Have you tried using an updated Goalie? We just made some changes to the error indicator code that might have an effect? I just pushed (most of) your changes to branch 55_multiple_var_ee, which is on top of those changes, so maybe you could checkout a local copy and make that one edit for the transfer.
• I wouldn't recommend computing the Hessian of the speed because it's non-differentiable at 0. I prefer to compute the Hessians of each velocity component and then combine them.
• More generally, the anisotropic_dwr formulation is designed to take a Hessian which is representative of the entire solution. That is, I think it would be better if you computed Hessians of each scalar solution field and combined them all together, rather than just considering different components.

[ddundo]

Thanks a lot Joe! I was indeed using the updated Goalie, but computing the Hessian of each velocity component and intersecting them, as you suggested, gave me a better looking mesh (3241 cells):

and then intersecting the thickness field as well (3211 cells):

So I'd say I'm happy with this now! Thanks a lot!

But could you please points me towards some paper or explain briefly how " the anisotropic_dwr formulation is designed to take a Hessian which is representative of the entire solution"? I couldn't find this explained in the Carpio et al. paper that you reference.

[jwallwork23]

Okay great, in that case I'll close this issue. Oh yeah I suppose the paper only mentions the scalar case. I guess the formulation is built around having a Hessian that represents the PDE solution so for the case with multiple scalar components to the equation, you want a Hessian that captures the solution for each component? This could probably be proved rigorously.