What is the recommended order for enforce_spd, normalize and intersect with multiple fields?

[stephankramer]

In #224, @Patol75 wrote:

It took me a while, but I have something that seems to be working. However, given that it is the first time I get to that level of complexity with animate, I would not mind if you could give my workflow a quick look. In particular, I am not exactly sure when and in which order normalise, enforce_spd, and compute_hessian should be called. The function I use if the following:

def adapt_mesh(
    mesh: fd.MeshGeometry, mesh_fields: dict[str, dict[str, fd.Function | bool]]
) -> fd.MeshGeometry:
    for _ in range(prms.adapt_calls):
        M = fd.TensorFunctionSpace(mesh, "CG", 1)

        metric_fields = []
        for field, field_specs in mesh_fields.items():
            if isinstance(field_specs["add_to_metric"], list):
                for dim, add_to_metric in enumerate(field_specs["add_to_metric"]):
                    if add_to_metric:
                        metric_fields.append(field.subfunctions[dim])
            elif field_specs["add_to_metric"]:
                metric_fields.append(field)

        for field in metric_fields:
            if isinstance(field.ufl_element(), fd.VectorElement):
                metric_fields.remove(field)
                for dim in range(field.ufl_shape[0]):
                    metric_fields.append(field.sub(dim))

        metrics = []
        for field in metric_fields:
            # Firedrake function for a metric over a mesh where a field lives
            metric = RiemannianMetric(M, name=f"Metric ({field.name()})")
            metric.set_parameters(prms.metric_parameters)  # Set metric parameters
            metric.compute_hessian(field)  # Field Hessian
            metric.normalise()  # Rescale metric to achieve the desired target complexity
            metric.enforce_spd()  # Ensure boundedness (symmetric positive-definite)
            metrics.append(metric)

        overall_metric = metrics[0].copy(deepcopy=True)  # Overall metric
        overall_metric.rename("Metric (overall)")
        # Minimum in all directions across all metrics
        overall_metric.intersect(*metrics[1:])
        overall_metric.enforce_spd()  # Ensure boundedness (symmetric positive-definite)

        mesh = adapt(mesh, overall_metric)  # Generate new mesh based on overall metric
        # Ensure boundary coordinates are not exceeded
        mesh.coordinates.dat.data[:, 0].clip(
            0.0, prms.domain_dims[0], out=mesh.coordinates.dat.data[:, 0]
        )
        mesh.coordinates.dat.data[:, 1].clip(
            0.0, prms.domain_dims[1], out=mesh.coordinates.dat.data[:, 1]
        )
        mesh.cartesian = True  # Geometry tag informing other G-ADOPT objects

        fields = list(mesh_fields)
        for field in fields:
            field_specs = mesh_fields.pop(field)

            if isinstance(
                field.function_space().topological,
                fd.functionspaceimpl.MixedFunctionSpace,
            ):
                spaces = []
                for dim, element in enumerate(field.ufl_element().sub_elements):
                    spaces.append(fd.FunctionSpace(mesh, element))

                mixed_space = fd.MixedFunctionSpace(spaces)
                new_field = fd.Function(mixed_space, name=field.name())
                for sub_field, new_sub_field in zip(
                    field.subfunctions, new_field.subfunctions
                ):
                    new_sub_field.rename(sub_field.name())
                    new_sub_field.interpolate(sub_field)
            else:
                space = fd.FunctionSpace(mesh, field.ufl_element())
                new_field = fd.Function(space, name=field.name())
                new_field.interpolate(field)

            mesh_fields[new_field] = field_specs

    return mesh

mesh_fields is defined as follows:

mesh_fields = {
    stokes: {"add_to_metric": [True, False]},
    T: {"add_to_metric": True},
    psi: {"add_to_metric": True},
}

Thank you for your time and help.

[stephankramer]

Yes, this is a confusing topic (it is for me at least) and I don't think there's an easy answer that's correct for all purposes.

The way it's done above, where you normalize the metric for each field (component) separately - you are basically assembling a metric for each field individually first, which answers the question: for each field individually what is the required edge length everywhere (and in all directions) to minimize the error norm in that field and such that the resulting mesh would have roughly "target complexity" number of elements. So imagine separate meshes that are optimal for each field individually, all with roughly the same n/o elements. Then intersect will everywhere (and in all directions) take the minimum of the edge length in these individual meshes. So you end up with a combined mesh that in general will actually have more elements. Thus if you want to obtain a mesh with the desired target complexity, you should call normalize again on the final intersected mesh.

This is all perfectly valid and may give you a good mesh in many cases. The advantage is that you don't really need to think how to scale the Hessian to get a desired estimated error, or even how big that desired error should be - that's all taken care of by normalize. However it does mean that you assume that the number of elements you need for each field individually to obtain the desired error for that field, is the same for each field - and that in practice is often not the case. Think of your (@Patol75) typical use case: you have the components of velocity of Stokes flow driven by some density (related) scalar fields. Initially your velocity might be (near) zero and you just want to adapt to sharp jumps in your density. If you were to apply the procedure above, the metric derived from the velocity components would be more or less uniform, and have an edge length determined by what's needed to achieve the target complexity. Then when combining it with the metric for the scalar density field, that edge length would act like a maximum edge length, meaning that it maintains where the density-based metric wants to refine, but where the density-based metric wants to coarsen it just applied a blanket uniform resolution - and finally if we normalize the combined intersected metric, it would need to use larger edge lengths in the refined areas.

The only way to deal with that is to actually consider how important you consider the different errors to be relative to each other. So this could be based on the way of thinking about an "estimated interpolation error" - as you are familiar with in Fluidity. Basically you look at the 2nd order Taylor term: $$f(x+h) = f(x) + f'(x) h + \tfrac 12 h^T H h$$ if we were to use a metric $$M=\epsilon^{-1} H$$ which gives us a mesh with edge lengths/vectors $$h$$ s.t. $$h^T M h\approx 1$$ then we can get any desired value for the "interpolation error" $$h^T H h=\epsilon$$ by choosing $$\epsilon$$, i.e. rescaling the metric by $$\epsilon^{-1}$$, where $$\epsilon$$ is the target "estimated interpolation error". If you do this for each field (component) individually, you typically want to choose a different value for $$\epsilon$$ for each - for starters it will depend on the dimensionalisation of the field (e.g. velocity units v.s. density units). You shouldn't then normalize these before intersection, as that rescales it again by something completely different, but you can then still normalize after intersection, which means that you don't need to worry about the absolute value of each $$\epsilon$$ it only matters what the relative value is between the different fields.

[stephankramer]

So that's about the order of normalize vs. intersect. Then there's also the order of enforce_spd vs. normalize.

So enforce_spd can do a number of related things, all to do with eigenvalues (look at the arguments to enforce_spd()):

1. ensure all eigenvalues are positive
2. restrict the edge lengths are between $$h_{min}$$ and $$h_{max}$$ by ensuring $$\lambda_{min}=1/h_{max}^2 \leq \lambda \leq \lambda_{max}=1/h_{min}^2$$
3. restrict the anisotropy, which can be translated into a bound on the ratio $$\lambda_{max}/\lambda_{min}$$

You kind of need 1) for any of the other operations to be valid. So I would do that immediately after compute_Hessian. Then 2) and 3) are done automatically inside normalize, by default it does 3) before the normalize, and 2) afterwards. This means that the number of elements can be quite different from the target if the minimum or maximum bound is hit in a lot of places. So if you care more about controlling the n/o elements, you might want to do it (in a separate call to enforce_spd) before normalize(), but then you won't actually get the min. and max. edge lengths that you've imposed as they'll be rescaled afterwards.

[stephankramer]

@joewallwork : any further thoughts/considerations? This is just my thoughts - but I'm more used to the approach where you don't do any normalisation, you just choose $$\epsilon$$ and apply it as an infinity norm, which often makes it a bit simpler to reason about what happens when combining and restricting in different orders. But I do actually agree that in practice really what we want is the best mesh with a chosen n/o elements, and also that the "best" may be expressed better in a L1 or L2 norm than $$L_\infinity$$ - so hence my inclination to only normalize at the end...

[joewallwork]

This (and previous discussions with @stephankramer and others) make me think we should provide some demos that go into this in this level of detail. Let's use this Issue as a marker to do that.

A few additional comments:

• I added a PETSc option that's exposed in Animate via the RiemannianMetric.set_parameters interface: dm_plex_metric_restrict_anisotropy_first. This controls whether anisotropy is restricted before or after metric normalisation. It defaults to True. Perhaps play around with that to see what works best for you.
• A workflow that I've often followed is to first call enforce_spd on each metric to be combined without imposing sizes or anisotropy, combine them, normalise, and then enforce_spd imposing sizes and anisotropy. Again, I don't think there's necessarily one best way to do this, so experiment to see what works best.
• In lieu of Animate demos, you might find it useful to look at some of the Goalie demos, which use Animate. Goalie implements space-time normalisation but otherwise the methods are the same. In particular, see https://github.com/mesh-adaptation/goalie/blob/main/demos/burgers-hessian.py.

[Patol75]

Thank you, Stephan and Joe, for your help here. I have started to have a look at your suggestions. I have not had enough time yet to investigate properly, but I am seeing an improvement towards some aspects I would like to see in the mesh if I apply normalise only after intersect and rescale each field metric after enforce_spd. More from me tomorrow.

[Patol75]

As I said yesterday, I have so far investigated Stephan's suggestions as a first attempt. I have implemented the scaling idea and tried to downscale the velocity metric relative to others. I have also moved normalise and enforce_spd calls based on your advice and actual results. My code is now as follows:

Updated code

mesh_fields = {
    stokes: {"add_to_metric": [True, False], "scaling": [1e-6, None]},
    T: {"add_to_metric": True, "scaling": 1e-1},
    psi: {"add_to_metric": True, "scaling": 1e0},
}
if prms.free_surface:
    mesh_fields[stokes]["add_to_metric"].append(False)
    mesh_fields[stokes]["scaling"].append(None)

adapt_calls = 3
metric_parameters = {  # For further information: `set_parameters` in animate/metric.py
    "dm_plex_metric": {
        "target_complexity": 20_000,  # Metric complexity, analogous to vertex count
        "h_min": 3e3 / distance_scale,  # Minimum metric magnitude (i.e. cell size)
        "h_max": 5e5 / distance_scale,  # Maximum metric magnitude (i.e. cell size)
        "a_max": 10.0,  # Maximum metric anisotropy (cell aspect ratio)
        "p": 1,  # Metric normalisation order
        "gradation_factor": 1.3,  # Maximum variation in length between adjacent edges
    }
}

def interpolate_fields(
    mesh: fd.MeshGeometry,
    mesh_fields: dict[str, dict[str, fd.Function | bool]],
    new_fields: list[fd.Function] | None = None,
) -> None:
    if new_fields is not None:
        for field, new_field in zip(mesh_fields, new_fields):
            field_specs = mesh_fields.pop(field)
            new_field.interpolate(field)
            mesh_fields[new_field] = field_specs

        return

    fields = list(mesh_fields)
    for field in fields:
        field_specs = mesh_fields.pop(field)

        if isinstance(
            field.function_space().topological,
            fd.functionspaceimpl.MixedFunctionSpace,
        ):
            spaces = []
            for dim, element in enumerate(field.ufl_element().sub_elements):
                spaces.append(fd.FunctionSpace(mesh, element))

            mixed_space = fd.MixedFunctionSpace(spaces)
            new_field = fd.Function(mixed_space, name=field.name())
            for sub_field, new_sub_field in zip(
                field.subfunctions, new_field.subfunctions
            ):
                new_sub_field.rename(sub_field.name())
                new_sub_field.interpolate(sub_field)
        else:
            space = fd.FunctionSpace(mesh, field.ufl_element())
            new_field = fd.Function(space, name=field.name())
            new_field.interpolate(field)

        mesh_fields[new_field] = field_specs


def adapt_mesh(
    mesh: fd.MeshGeometry, mesh_fields: dict[str, dict[str, fd.Function | bool]]
) -> fd.MeshGeometry:
    output_file = fd.VTKFile("output_adapt_mesh.pvd", adaptive=True)
    output_fields = []
    for name in ["Stokes", "Temperature", "Level set"]:
        for field in mesh_fields:
            if field.name() != name:
                continue

            if isinstance(
                field.function_space().topological,
                fd.functionspaceimpl.MixedFunctionSpace,
            ):
                for sub_field in field.subfunctions:
                    output_fields.append(sub_field)
            else:
                output_fields.append(field)
    write_output(output_file, 0, *output_fields)

    for _ in range(prms.adapt_calls):
        M = fd.TensorFunctionSpace(mesh, "CG", 1)

        metric_fields = {}
        for field, field_specs in mesh_fields.items():
            if isinstance(field_specs["add_to_metric"], list):
                for dim, (add_to_metric, scaling) in enumerate(
                    zip(field_specs["add_to_metric"], field_specs["scaling"])
                ):
                    if add_to_metric:
                        metric_fields[field.subfunctions[dim]] = scaling
            elif field_specs["add_to_metric"]:
                metric_fields[field] = field_specs["scaling"]

        fields = list(metric_fields)
        for field in fields:
            if isinstance(field.ufl_element(), fd.VectorElement):
                scaling = metric_fields.pop(field)
                for dim in range(field.ufl_shape[0]):
                    metric_fields[field.sub(dim)] = scaling

        metrics = []
        for field, scaling in metric_fields.items():
            # Firedrake function for a metric over a mesh where a field lives
            metric = RiemannianMetric(M, name=f"Metric ({field.name()})")
            metric.set_parameters(prms.metric_parameters)  # Set metric parameters
            metric.compute_hessian(field)  # Field Hessian
            metric.enforce_spd()  # Ensure boundedness (symmetric positive-definite)
            metric.assign(metric * scaling)
            metrics.append(metric)

        overall_metric = metrics[0].copy(deepcopy=True)  # Overall metric
        overall_metric.rename("Metric (overall)")
        # Minimum in all directions across all metrics
        overall_metric.intersect(*metrics[1:])
        overall_metric.enforce_spd()  # Ensure boundedness (symmetric positive-definite)
        overall_metric.normalise()  # Rescale metric to achieve the desired target complexity

        mesh = adapt(mesh, overall_metric)  # Generate new mesh based on overall metric
        # Ensure boundary coordinates are not exceeded
        mesh.coordinates.dat.data[:, 0].clip(
            0.0, prms.domain_dims[0], out=mesh.coordinates.dat.data[:, 0]
        )
        mesh.coordinates.dat.data[:, 1].clip(
            0.0, prms.domain_dims[1], out=mesh.coordinates.dat.data[:, 1]
        )
        mesh.cartesian = True  # Geometry tag informing other G-ADOPT objects

        interpolate_fields(mesh, mesh_fields)

        output_fields = []
        for name in ["Stokes", "Temperature", "Level set"]:
            for field in mesh_fields:
                if field.name() != name:
                    continue

                if isinstance(
                    field.function_space().topological,
                    fd.functionspaceimpl.MixedFunctionSpace,
                ):
                    for sub_field in field.subfunctions:
                        output_fields.append(sub_field)
                else:
                    output_fields.append(field)
        write_output(output_file, _ + 1, *output_fields)

    for field in output_fields:
        assert mesh is field.ufl_domain()

    return mesh

Here is the generated mesh and level-set fields in the initial state:

For reference, to obtain the initial state, I generate a starting mesh from Gmsh, impose initial conditions on transported fields (Temperature and Level set), impose boundary conditions on other relevant fields (Velocity), solve the Stokes system, adapt the mesh, and finally reapply initial conditions on transported fields.

Essentially, what I do not understand is why I cannot manage to have the intersected mesh intentionally biased towards the level-set field. For comparison, here is the mesh when I adapt only on the level-set field:

[Patol75]

I will also note that despite the changes made so far to the adaptivity strategy, I still get the following mesh when commenting out the lines related to VTK output. I had already mentioned this to Stephan on a separate channel, but I just wanted to document it here. I will have to investigate.

adapt_mesh with VTK output commented out

def adapt_mesh(
    mesh: fd.MeshGeometry, mesh_fields: dict[str, dict[str, fd.Function | bool]]
) -> fd.MeshGeometry:
    # output_file = fd.VTKFile("output_adapt_mesh.pvd", adaptive=True)
    # output_fields = []
    # for name in ["Stokes", "Temperature", "Level set"]:
    #     for field in mesh_fields:
    #         if field.name() != name:
    #             continue

    #         if isinstance(
    #             field.function_space().topological,
    #             fd.functionspaceimpl.MixedFunctionSpace,
    #         ):
    #             for sub_field in field.subfunctions:
    #                 output_fields.append(sub_field)
    #         else:
    #             output_fields.append(field)
    # write_output(output_file, 0, *output_fields)

    for _ in range(prms.adapt_calls):
        M = fd.TensorFunctionSpace(mesh, "CG", 1)

        metric_fields = {}
        for field, field_specs in mesh_fields.items():
            if isinstance(field_specs["add_to_metric"], list):
                for dim, (add_to_metric, scaling) in enumerate(
                    zip(field_specs["add_to_metric"], field_specs["scaling"])
                ):
                    if add_to_metric:
                        metric_fields[field.subfunctions[dim]] = scaling
            elif field_specs["add_to_metric"]:
                metric_fields[field] = field_specs["scaling"]

        fields = list(metric_fields)
        for field in fields:
            if isinstance(field.ufl_element(), fd.VectorElement):
                scaling = metric_fields.pop(field)
                for dim in range(field.ufl_shape[0]):
                    metric_fields[field.sub(dim)] = scaling

        metrics = []
        for field, scaling in metric_fields.items():
            # Firedrake function for a metric over a mesh where a field lives
            metric = RiemannianMetric(M, name=f"Metric ({field.name()})")
            metric.set_parameters(prms.metric_parameters)  # Set metric parameters
            metric.compute_hessian(field)  # Field Hessian
            metric.enforce_spd()  # Ensure boundedness (symmetric positive-definite)
            metric.assign(metric * scaling)
            metrics.append(metric)

        overall_metric = metrics[0].copy(deepcopy=True)  # Overall metric
        overall_metric.rename("Metric (overall)")
        # Minimum in all directions across all metrics
        overall_metric.intersect(*metrics[1:])
        overall_metric.enforce_spd()  # Ensure boundedness (symmetric positive-definite)
        overall_metric.normalise()  # Rescale metric to achieve the desired target complexity

        mesh = adapt(mesh, overall_metric)  # Generate new mesh based on overall metric
        # Ensure boundary coordinates are not exceeded
        mesh.coordinates.dat.data[:, 0].clip(
            0.0, prms.domain_dims[0], out=mesh.coordinates.dat.data[:, 0]
        )
        mesh.coordinates.dat.data[:, 1].clip(
            0.0, prms.domain_dims[1], out=mesh.coordinates.dat.data[:, 1]
        )
        mesh.cartesian = True  # Geometry tag informing other G-ADOPT objects

        interpolate_fields(mesh, mesh_fields)

    #     output_fields = []
    #     for name in ["Stokes", "Temperature", "Level set"]:
    #         for field in mesh_fields:
    #             if field.name() != name:
    #                 continue

    #             if isinstance(
    #                 field.function_space().topological,
    #                 fd.functionspaceimpl.MixedFunctionSpace,
    #             ):
    #                 for sub_field in field.subfunctions:
    #                     output_fields.append(sub_field)
    #             else:
    #                 output_fields.append(field)
    #     write_output(output_file, _ + 1, *output_fields)

    # for field in output_fields:
    #     assert mesh is field.ufl_domain()

    return mesh

[Patol75]

Hi @stephankramer and @joewallwork,

I am sorry it took me an awful long time, but I am finally here with good news. I have managed to obtain very good mesh generation from animate (GitHub is not letting me share the animation, but I have put it on Zulip for now @stephankramer). In the last two months, quite a few bugs across Firedrake, PETSc, animate, and my code had to be fixed to allow mesh evolution to behave properly. In the end, I took the road suggested by @stephankramer: rescale each metric before normalising. My current code is available here (dismiss main.py and simulation.py). The parameters provided to animate are here, and the class methods used to adapt the mesh and perform cross-mesh interpolation are here. Thank you again for your support, patience, and ideas.