The underlying SVD for definition of sectional DVs can give inconsistent results.

[anilyil]

Description

SVD (U \Sigma V^T = A) is used to figure out a local coordinate frame. This is done by looking at the columns of the U matrix. The first column is a best fit line through the points that define that FFD section. The second column is the next best fit that is orthogonal to the first one; this ends up being the second vector that defines the section tangents (although not exactly). Then the third one is the last orthogonal vector that spans the 3 dimensional space, and that ends up being orthogonal to the section surface itself.

The SVD is only unique up to a sign change of these vectors in the U and V matrices. So while the U1 U2 and U3 vectors’ orientations are correct, their “sense” (what way the arrow is pointing to) can vary based on what machine its running on. I realized that the sectional DVs on the dlr f6 case ended up getting a wrong direction on some of the sections when I ran it on docker, vs natively on the same computer. So far, seems like the SVD has been giving consistent results for everyone, but it can totally flip direction if not careful. I think even if you provide an “orient0" direction, the other two directions can still flip. If you do not provide any orient directions, then all 3 vectors can flip either way.

The method to create sectional DVs: https://github.com/mdolab/pygeo/blob/master/pygeo/DVGeometry.py#L689

During its setup, it calls https://github.com/mdolab/pygeo/blob/master/pygeo/DVGeometry.py#L3407 to figure out the sectional coordinates.

Here is the SVD where we take the U matrix from: https://github.com/mdolab/pygeo/blob/master/pygeo/DVGeometry.py#L3505

The sectionFrame method can take in an optional argument orient0, which is a way to correctly orient the first basis vector. However, this does not fully fix the solution as the other two vectors can still flip among themselves. If orient0 is not provided, then all 3 vectors can flip direction.

Suggested Solution

My suggested solution is to make sure we align these vectors with some of the directions introduced by the FFD. So I am proposing computing the “finite difference” version of these vectors with the FFD, and then using the sign of the dot product of these FFD vectors with the SVD vectors to flip the SVD vectors if required.

At least for Ney’s F6 stuff, this fix flips the direction of the sectional DV changes. However, this fix does make it consistent between my native and docker runs.

My suggested solution will most likely break the runscripts that have been using this in the past.

Steps to reproduce issue

I am attaching a simple FFD box and a python test script. There are 16 total FFD control points we select and we define normal DVs using these. When run with the current default branch of pygeo (commit: f8176b01a2bd3a01b59078db3048fec261d3bf45) the control points on the upper side ends up flipping in direction when I run it natively on my mac vs docker running on my mac.

pygeo_sectional_test.zip

Current behavior

The SVD can result in vectors that change in sense between computers, and this can lead to "flipping" of design variable directions.

Expected behavior

The local coordinate system computed for FFD sections should be unique.

[anilyil]

Any suggestions @nbons? I believe you worked on this at some point in the past?

[nbons]

I think your solution sounds good. It should be noted that if you use orient2='ffd', it uses the finite difference approach to compute the normal vector (I believe Gaetan added this as an option). This should work okay in most standard wing cases, but might give the wrong normal if you have sections that are changing orientation rapidly. I think correcting the sign of the SVD with the FFD approximation is probably the best "automatic" way to determine the normal.

In terms of breaking old cases, I don't think we should worry about that. This is definitely the better way forward, and it should really only flip the signs of the final design variables.

[anilyil]

Right, the 'ffd' version is nice, but the SVD is more accurate as you might have guessed. I think if the user specifies both orient2='ffd' and orient0 to some value, then the coordinates are fully fixed because I think the code does an internal cross product of these vectors to get the third one.

If orient2 is not provided, then the orientation of the normal vector can flip, causing the reference frame to flip (this is true regardless of the orient0 option).

If orient0 is not provided, then the code takes the first two vectors from SVD as is, and then the order can even be wrong (u0 cross u1 may not be in the same direction as u2).

Your comment helps a lot because I think you know about this with much more detail. I will work on the fix and reference this issue when I open the PR. Thanks for the help @nbons!