DiscretizationCollection._base_to_geoderiv_connection is broken when using overintegration in lazy-mode

[thomasgibson]

The relatively recent change for storing affine geometric terms as constants for the pytato array context is broken when using overintegration. The code currently assumes that the from_discr is interpolatory, which is just not true when you have terms defined on a quadrature grid.

I have attempted to circumvent this issue in #172 but I would really like to settle this once and for all. And by this, I mean actually taking advantage of affineness even when overintegrating. Especially since having working overintegration with lazy evaluation is a high priority right now.

@inducer --- I don't know if I'm being dense here, but I don't see this as a simple 2-line change in the method. The problem is getting meshmode to do the right thing by transferring a DOFArray defined on a quadrature grid to the appropriate shape for the geo_deriv_discr. Is there something we can do in grudge to spoon-feed meshmode something it can work with?

[alexfikl]

With absolutely no actual knowledge of the problem..

Don't all geometry derivatives need to be computed on the base discretization and then interpolated to the quadrature grid from there? How is _base_to_geoderiv_connection getting a quad_dd?

[thomasgibson]

Don't all geometry derivatives need to be computed on the base discretization and then interpolated to the quadrature grid from there? How is _base_to_geoderiv_connection getting a quad_dd?

The derivatives are interpolated to the quadrature grid here: https://github.com/inducer/grudge/blob/main/grudge/geometry/metrics.py#L142-L143

The reason _base_to_geoderiv_connection is getting a quad_dd is in the operator application routines: https://github.com/inducer/grudge/blob/main/grudge/op.py#L427-L429 (and similar places for the weak derivative operator functions). This makes its way down to this part of the metric computation: https://github.com/inducer/grudge/blob/main/grudge/geometry/metrics.py#L423-L424, which is ultimately what is triggering the error.

The incoming dof descriptor must be on the quadrature grid when using overintegration (dd_in) and for the kernels to be correct, you need to have your geometry terms at each quadrature node (if non-affine).

When using affine, we just need to broadcast a constant value to each quad node.

[inducer]

Do you have a branch that exhibits the failure? Based on pure speculation, I suspect that the issue might be that you're asking the metrics to return a value on the quadrature grid. I don't think they're currently able to do that, but we should be able to teach them (simply compute the derivative on the base and then upsample that). (And that in turn would remove a bunch of repeated computation in the overintegrated case!)

[alexfikl]

That makes sense. Thank you for the detailed explanation!

Maybe an easy fix would be to just disable the interpolation to the quadrature grid (in the metrics) when using the base_to_geoderiv_connection? It's going to reduce everything to 1-dof anyway.

[inducer]

disable the interpolation to the quadrature grid (in the metrics)

Can you point to code that does that? I don't think the metrics attempt to project to a quadrature grid (but I'm suggesting we teach them how to do that).

[alexfikl]

Can you point to code that does that? I don't think the metrics attempt to project to a quadrature grid (but I'm suggesting we teach them how to do that).

From the code @thomasgibson linked in here? It seems to remove the quadrature tag, do the derivatives, and then interpolate back, right?

[inducer]

Ah, yep! Do all the metrics do that consistently? (I searched for project and hence came up empty.)

[alexfikl]

Ah, yep! Do all the metrics do that consistently? (I searched for project and hence came up empty.)

They should, yeah. All the other metrics just use the result from forward_metric_nth_derivative in various ways.

[thomasgibson]

Do you have a branch that exhibits the failure?

@inducer You can hit the exact error this issue is describing by running examples/wave/wave-op-var-velocity.py in lazy mode (python wave-op-var-velocity.py --lazy) --- use this branch: https://github.com/inducer/grudge/tree/lazy-wave-op-oi

[thomasgibson]

I believe I fixed this problem in this commit: 2c9c90c

In summary, this commit completely removes the _use_geoderiv_connection arg and instead metric routines will call the geometry discretization connection if the underlying actx (which is always passed through to the functions) supports nonscalar broadcasting and if the discretization has affine element groups (this information is provided directly by the DiscretizationCollection object).

Interpolation to the quadrature grid happens at the very last moment (just before memoization), and it will only happen if the actx doesn't support broadcasting or the mesh contains non-affine element groups.

Thoughts on this approach?

[inducer]

Can you explain the context of the approach a bit better? What does the removal of the _use_geoderiv_connection have to do with the issue you encountered?

[thomasgibson]

Just summarizing what @inducer and I discussed: The approach taken in #172 fixes this issue by moving the interpolation step (to a quadrature grid) at the very last moment (near the user-facing level in the relevant geometry routines; area_element, inverse_surface_metric_derivative_mat, mv_normal, and shape_operator specifically).

No major removals/additions were necessary 😄