How to modify the coordinates in the refinement hierarchy

[jiangzhongshi]

Hi,

In the geometric multigrid experiment, I wish my final level to be a better approximation towards the target shape. In practice, what I want is something like this. Here I only globally enlarge the mesh by 1.1, but in practice, I may want to perturb the vertices based on some other scheme.

mesh = UnitSquareMesh(8, 8)
hierarchy = MeshHierarchy(mesh, 4)
mesh = hierarchy[-1]
mesh.coordinates.dat.data[:, 1] *= 1.1

Here, the multigrid solver {"ksp_type": "cg", "pc_type": "mg"}gives the following message

Could not identify cell in transfer operator. Point: 1.56250000000000e-02 8.59375000000000e-01
Number of candidates: 4. Best distance located:   1.000000e+01[1]    86754 abort (core dumped)

My guess is that the coordinates are no longer nested in the physical space, but if I understand correctly, the transfer should be done in the parametric space, or is there something obvious that I missed?

Thanks! Zhongshi

[wence-]

The transfer has to happen in physical space. We need to know where the points are so as to know what the value of the coarse correction is on the fine grid.

We have some tolerance to allow for extrapolation near boundaries but it looks like your perturbation was too much in this case.

If what you want is that your mesh boundary conforms to some geometry (eg CAD of some kind), then you need to do this "as best as possible" on all levels to avoid the issue you observed in your test. The OpencascadeMeshHierarchy may be useful here. It builds a multigrid hierarchy with refinement, but takes a .step file and pushes the boundary nodes out to that geometry on every level.

[jiangzhongshi]

Hi Lawrence,

Thanks for the clarification.

I am quite new to the field, is there a reference for the interpolation/restriction scheme used in Firedrake? And is it possible to have parametric space-based transfer, and do you have some insight on how it affects convergence or performance?

The OpencascadeMeshHierarchy may be useful here.

Yes, that sounds more like my use case, is the transfer there also done in the physical space, and if so, would it fail if the base geometry gets coarsened too much?

After some debugging, the program complained in the function def inject_kernel(Vf, Vc):. (I am only starting to learn multigrid so please correct me if I am saying something stupid.) It seems to me that this function only queries the coordinate to determine the cell that it resides in, but no barycentric coordinates are computed (which I guess would be useful for the interpolation in physical space). Is this the place where the physical space interpolation is employed?

[wence-]

I am quite new to the field, is there a reference for the interpolation/restriction scheme used in Firedrake? And is it possible to have parametric space-based transfer, and do you have some insight on how it affects convergence or performance?

We do, effectively, the "standard" finite element approach. Let's consider prolongation/interpolation (restriction is dual to this, in the linear algebra sense its the transpose operator).

For prolongation we have a correct uc on a coarse grid and we want to represent it on the fine grid. uc can be evaluated at any point on the coarse mesh. So really what we need to to in the prolongation is figure out for each fine grid basis function, what the corresponding dof value is. This can be done by acting each fine grid dual basis on the coarse grid function in turn. Right now we only support dual bases that are single "point evaluations". So what this means is that for each fine grid basis function we do an evaluation of the coarse function uc at a point.

If the mesh hierarchy were created by regular refinement with midpoint bisection (so that the nesting is perfect both topologically and geometrically) then this evaluation can be encoded by a little matrix that doesn't depend on the geometry (or at least it does, but it's the same everywhere). But imagine the case where some of the newly introduced fine grid vertices are not at edge midpoints (for example you moved them at the boundary). Now even though the topologies still nest, the geometry is not quite as regular: if your fine grid vertex is 1/3rd the way along a coarse edge, then you can't use the same interpolation as the case where it's halfway along.

Our solution to this is to do the point evaluation by finding the location of each basis function in physical space, and then finding the corresponding coarse grid cell, mapping into reference space on the coarse grid cell and evaluating uc there.

Naturally, this evaluation is only technically valid if we find a point that is inside a coarse mesh cell. However, as long as we're "not too far away" then extrapolation doesn't really affect the performance of multigrid: the abstract theory of non-nested multigrid from Bramble Pasciak and Xu provides bounds (although they are not that tight).

tl;dr: the prolongation attempts to find, for each fine grid basis function, the coarse cell it is located in, and then evaluates the coarse function uc to get the value of the dof. If we don't find a cell, we look for the "nearest" one, and do extrapolation (just evaluate uc "in" that closest cell).

If we could always guarantee that our mesh hierarchies were completely regular, there could be a "fast path" that did a simpler thing.

Yes, that sounds more like my use case, is the transfer there also done in the physical space, and if so, would it fail if the base geometry gets coarsened too much?

The failure occurs if the locations at which we are evaluating fall too far away from any cell. So I think it could still fail in the same way, but is less likely too.

After some debugging, the program complained in the function def inject_kernel(Vf, Vc):. (I am only starting to learn multigrid so please correct me if I am saying something stupid.) It seems to me that this function only queries the coordinate to determine the cell that it resides in, but no barycentric coordinates are computed (which I guess would be useful for the interpolation in physical space). Is this the place where the physical space interpolation is employed?

Everything ends up being done in (or at least via) physical space. But it's possible you passed "unsupported" functions to injection. Perhaps you can show a sketch of what you're attempting in code.

[jiangzhongshi]

Thanks so much for the thorough explanation! It really helps a lot, as cubes/squares in my textbook often lead me to consider only things happening in the parametric space with a canonical stencil.

if your fine grid vertex is 1/3rd the way along a coarse edge, then you can't use the same interpolation as the case where it's halfway along.

Yes, the explanation makes a lot of sense. But if the transfer happens in physical space, the function space of the transfer is no longer nested? Is Firedrake designed/suitable for non-nested multigrid?

On the other hand, for a case where everything is enlarged/perturbed a little bit between hierarchies, physical space transfer will result in extrapolation error outside the domain and parametric space transfer results in "amortized" error everywhere, do you have insight of reference if one case is better than the other? (apparently for large enough perturbation they both fail)