Efficient and Accurate Coordinate Transform

Hello everyone,

I've been working a adaptive finite difference solving that can handle non-orthogonal meshes for arbitrary dimensions. This is being used to calculate the derivative terms of a set of vector fields for some fluid flow in a steady domain. Because the domain is steady and I'm typically looking at tens or hundreds vector fields. I'm trying to generate a set of coefficients upfront based on the characteristics of the mesh so I create a set of shifted vector fields and sum them together to get derivative of the vector field very cheaply.

Where I'm looking for suggestions is for a coordinate transform. Originally I created this set of functions assuming a 2D vector field. There I explicitly calculated the jacobain inverse upfront since it was only a 2x2 matrix inverse then multiplied it to each of vector field derivative calculated on a 'computational' domain i.e. just a set of index. This worked well because I could calculate the inverse once for a mesh then apply it to all of my flow snapshots. So now extending this now to 3D is causing problem with running into singularities and ill-conditioned matrices. So I know I could always produce the transform by using backslash and avoid essentially all of these problems, but then to my knowledge I would need to do this with very vector field. So is there a more stable means of calculating a 3x3 inverse than inv. While the end goal is to perform the transform, getting the inverse jacobain explicitly allows a massive speedup for my application.

Answer 1

[Matt Cohen]

Hi John,

I understand that you are interested in finding a stable method for computing a 3x3 inverse Jacobian.

If the matrix you are trying to find the inverse of is nearly singular and ill-conditioned, then you might want to try computing the Moore-Penrose pseudoinverse of the matrix using the function "pinv". The computation for the pseudoinverse is based on the matrix's singular value decomposition (SVD); it determines the singular values that would lead to singularity and removes these from the process to produce a more numerically stable result.

This is probably a decent option for your problem since the matrix is a small size. Using the backslash operator to solve the system for this size of a problem is likely very efficient, but it sounds like it is less convenient for the process. There are also other factorizations that can be helpful for this kind of problem, like LU factorization, but it would follow a similar workflow to using the backslash operator.

I hope this helps.

Matt Cohen