Hi there,

I have a feeling this may have been asked before but I can’t seem to find an answer online. The general question is: how can I compute the gradient of a function where the function depends on the solution of a PDE and the derivatives are taken wrt to the shape of the boundaries of the PDE? The problem, as far as I can tell, is that mesh objects are not differentiable, so I can’t use automatic differentiation to just magic my way through this. And the naive solution I came up with seems like overkill. So I was hoping that maybe someone here could educate me

For concreteness, imagine I wanted to optimize the shape of an object that was falling through a tube of water. For any particular object boundary \vec{\partial\Omega}(\{a_\alpha\}), where \{a_\alpha\} is a set of parameters that define the shape, I should be able to: (1) generate a suitable mesh of the 3D solution volume (still not totally sure how one does this, but let’s avoid the tangent); (2) specify that the velocity evaluated on the surface of the object is a constant vector \vec{v}(\vec{x_i})=\vec{v_o} where \vec{x_i}\, \epsilon\,\vec{\partial\Omega} are the vertices of the surface mesh or something similar; (3) similarly specify that the velocity at the tube walls is exactly zero; (4) solve the hydrodynamics problem using your favorite solver; and (5) finally compute the drag force on the object (the square of which, say, is our objective function).

Naively, I should be able to leave the mesh as is and directly differentiate my PDE wrt to my shape parameters and solve for \partial \vec{v}(\vec{x}, a_\alpha)/\partial a_\alpha by running my PDE solver over and over again for every shape variable. Thus computing the derivatives of the shear stresses and so on. Again, for concreteness, this would mean solving

\partial \vec{v}(\vec{x}, a)/\partial a\cdot\nabla\vec{v}(\vec{x})+\vec{v}(\vec{x})\cdot\nabla\partial\vec{v}(\vec{x}, a)/\partial a - \nu\nabla^2\partial \vec{v}(\vec{x}, a)/\partial a=-\frac{1}{\rho}\nabla\partial p(\vec{x}, a)/\partial a

with boundary condition

\partial \vec{x_i}/\partial a\cdot \nabla \vec{v}(\vec{x_i}) + \partial \vec{v}(\vec{x_i}, a)/\partial a=0

over and over again for every shape parameter using the \nabla\vec{v}(\vec{x}) and \vec{v}(\vec{x}) obtained from the initial solution.

Mercifully at least these PDEs are all linear. But there could be dozens or hundreds of them (squared and no longer linear if you wanted Hessians). And if I wanted to use this procedure to optimize the shape parameters I would then also have to automate and wrap the meshing operation inside of a function that computed the final objective function and gradients (which also assumes that discrete changes in the mesh from one step to the next don’t introduce non-differentiable noise into the gradient). So, first, is this correct/feasible and, second, is there no better way to do this?