symbolic differentiation of scalars wrt tensors

[liamoconnor9]

It'd be nice if we could differentiate dot products with respect to their (vector) arguments. This would be particularly useful for the adjoint-looping optimization work I'm doing.

Here's an example of the functionality I'm looking for:

    coords = d3.CartesianCoordinates('x', 'z')
    dist = d3.Distributor(coords, dtype=np.float64)
    xbasis = d3.RealFourier(coords['x'], size=64, bounds=(0, 2*np.pi), dealias=3/2)
    zbasis = d3.Chebyshev(coords['z'], size=64, bounds=(-1, 1), dealias=3/2)
    x, z = dist.local_grids(xbasis, zbasis)

    u = dist.VectorField(coords, name='u', bases=(xbasis, zbasis))
    u['g'][0] = np.sin(x)*z
    u['g'][1] = np.cos(x)*z**2

    usqrd = d3.dot(u, u)

    assert np.allclose(usqrd.sym_diff(u).evaluate()['g'], 2*u['g'])

Given a TensorField u, taking u.sym_diff(u) returns 1. If instead it returned an identity tensor of the proper rank, then the test case above would pass

[kburns]

The correct way to do this is to compute the Frechet differential, not to use sym_diff.

[kburns]

Whoops didn't mean to close

[liamoconnor9]

The correct way to do this is to compute the Frechet differential, not to use sym_diff.

Ok my mistake. I think we can recycle a fair amount of existing code to make this in Dedalus. Would that be a useful feature? Or should I add that feature to my own project outside of Dedalus?

[kburns]

This already works -- it's how the Jacobian is computed for the NLBVP: https://github.com/DedalusProject/dedalus/blob/579b32b774eac77db9abbd1db9f546f171ece743/dedalus/core/problems.py#L397

[liamoconnor9]

Thanks for pointing me in the right direction, I'll take a look at that.